Running vacuum in QFT in FLRW spacetime: the dynamics of ρvac(H)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\textrm{vac}}(H)$$\end{document} from the quantized matter fields

Phenomenological work in the last few years has provided significant support to the idea that the vacuum energy density (VED) is a running quantity with the cosmological evolution and that this running helps to alleviate the cosmological tensions afflicting the Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}CDM. On the theoretical side, recent devoted studies have shown that the properly renormalized ρvac\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\textrm{vac}}$$\end{document} in QFT in FLRW spacetime adopts the ‘running vacuum model’ (RVM) form. While in three previous studies by two of us (CMP and JSP) such computations focused solely on scalar fields non-minimally coupled to gravity, in the present work we compute the spin-1/2 fermionic contributions and combine them both. The calculation is performed using a new version of the adiabatic renormalization procedure based on subtracting the UV divergences at an off-shell renormalization point M. The quantum scaling of ρvac\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\textrm{vac}}$$\end{document} with M turns into cosmic evolution with the Hubble rate, H. As a result the ‘cosmological constant’ Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document} appears in our framework as the nearly sustained value of 8πG(H)ρvac(H)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8\pi G(H)\rho _{\textrm{vac}}(H)$$\end{document} around (any) given epoch H, where G(H) is the gravitational coupling, which is also running, although very mildly (logarithmically). We find that the VED evolution at present reads δρvac(H)∼νeffmPl2H2-H02(|νeff|≪1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta \rho _\textrm{vac}(H)\sim \nu _{\textrm{eff}}\, m_{\textrm{Pl}}^2 \left( H^2-H_0^2 \right) \ (|\nu _{\textrm{eff}}|\ll 1)$$\end{document}. The coefficient νeff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _{\textrm{eff}}$$\end{document} receives contributions from all the quantized fields, bosons and fermions, which we compute here for an arbitrary number of matter fields. Remarkably, there also exist higher powers O(H6)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{O}(H^{6})$$\end{document} which can trigger inflation in the early universe. Finally, the equation of state (EoS) of the vacuum receives also quantum corrections from bosons and fermion fields, shifting its value from − 1. The striking consequence is that the EoS of the quantum vacuum may nowadays effectively appears as quintessence.


Introduction
Despite having coexisted for many decades, a completely successful theory of gravity that combines Quantum Field Theory (QFT) and General Relativity (GR) does not exist yet, unfortunately. However, a variety of different approaches and techniques are available in the literature which allow one to study the subject of quantum fields in the gravitational context, and more specifically the physics of the expanding universe and its current speeding up. Our aim is to understand such an acceleration on fundamental grounds. To be precise, in this work we will concentrate on the well-known semiclassical approach which goes under the name of QFT in curved spacetime [1][2][3]. This means that gravity is still a classical external (background) field, whereas the matter fields are quantum field operators obeying suitable commutation or anticommutation relations [4,5]. A further step in the path of understanding gravity in the QFT context is quantum gravity (QG), in which spacetime itself (the metric) is quantized and hence functional integration over metrics is mandatory, see e.g. [6][7][8], and the review [9] and references therein. At the same time a lot of exciting QG phenomenology is being investigated in the current multi-messenger era, characterized by an outburst of experimental data that are being obtained from the detection of the various cosmic messengers (photons, neutrinos, cosmic rays and gravitational waves) from numerous origins [10]. On the theoretical side, effective field theory methods and the possibility of quantum gravitational effects leading to quantum hair may provide useful hints of QG which have been explored recently [11][12][13]. However, while the QG option has, of course, to be kept in mind since it can be very important when QG can be (hopefully) formulated in a fully consistent way [9], QFT in curved spacetime may still be of great help to further describe the role of quantum fields in a gravitational context. In this work, we will continue dwelling upon these lines and shall focus exclusively on the, more modest, but effective, semiclassical approach. It goes without saying that the latter has had also its own problems and successes over the years, and still has [14]. However, new perspectives have recently been explored in this context concerning the vacuum energy and the cosmological constant (CC) [15][16][17] which may be of significance, and for this reason we wish to further pursue this line of approach here.
The agent responsible for the accelerated cosmic expansion is generically called Dark Energy (DE), an entity which constitutes a key piece in the cosmological puzzle, but whose fundamental nature is still undisclosed [18]. Although it might be caused by deviations from GR connected with extended theories of gravity [20][21][22][23][24], the canonical possibility is still that the DE is related to the CC in Einstein's equations, Λ, as done routinely in the standard (or 'concordance') model of cosmology, aka ΛCDM [25][26][27]. The model has been a rather successful paradigm for the phenomenological description of the universe for about three decades, but it became consolidated only in the late nineties [28,29] and especially after the consistent measurements of Λ made in the last twenty years using independent cosmological sources, in particular including distant type Ia supernovae (SnIa), baryon-acoustic oscillations (BAO), the data on large-scale structure formation and, of course, the anisotropies of the cosmic microwave background (CMB). All in all they have put the very experimental basis for the concordance ΛCDM model of cosmology [30][31][32][33][34][35]. The situation is far from being satisfactory, though. The problems with the ΛCDM are both of theoretical and observational nature. As for the theoretical problems, recall that the value of Λ is traditionally associated to a parameter called the vacuum energy density (VED) in the universe, which in the context of the ΛCDM is nothing but a name for the following quantity with dimensions of energy density: ρ vac = Λ/(8πG N ) (G N being Newton's constant). Its theoretical significance is not explained at all in the context of the standard cosmological model. If, however, we take quantum theory seriously, the most universal contribution to this vacuum energy density is the zero-point energy (ZPE) of the massive quantum fields in the standard model of particle physics, and in fact in any realistic QFT model. However, it is well-known that a naive calculation of this quantity leads to very large contributions proportional to the quartic power of the mass of the particles [36], ρ ZPE ∼ m 4 , which is in blatant discordance with the order of magnitude obtained for this quantity from cosmological observations: ρ obs vac ∼ 10 −47 GeV 4 (expressed in natural units, with = c = 1). Even taking, for instance, the electron field one finds a mismatch of 34 orders of magnitude: ρ obs vac /ρ ZPE ∼ 10 −34 . The huge discrepancy between a typical standard model contribution to the ZPE and the measured value of VED constitutes the so-called Cosmological Constant Problem (CCP) [37][38][39]. See also [40,41] for a recent account. Despite the enormous discrepancies between usual theory predictions and factual measurements, estimates on the value of Λ within the right order of magnitude have been attempted under certain assumptions in the context of QG in different approaches, see e.g. [42][43][44][45][46]. While the aforementioned measurements of ρ vac indicate that the vacuum can gravitate within an energy density order of magnitude of ∼ 10 −47 GeV 4 , what is difficult to understand theoretically is why the vacuum can only gravitate in that tiny range, given the fact that any typical quantum effect rockets its contribution to much larger values. This is of course a rephrasing of the same puzzle associated to the CCP, expressed in the QFT context. However, new avenues for a possible solution have been suggested recently. The renormalization approach presented in the present work and in the preceding studies [15][16][17] offers some hope to eschew part of these difficulties. First and foremost, the renormalized quantum effects found here endow the VED with a mild dynamical nature. The latter thus appears as a slowly varying function of the cosmic expansion, specifically of the Hubble rate H, see below. Second, the renormalized VED as reported here proves well behaved and can perfectly accommodate the measured value of Λ from observational cosmology without fine-tuning. Technically, this is because the "running" of Λ is proportional to the tiny values of the β-function coefficients for bosons and fermions, which are responsible for the renormalization group evolution of the VED. As a result, at any given epoch of the late universe Λ appears essentially as constant, but it is not strictly so. Finally, a third crucial ingredient of our approach is that, in the very early universe, the VED becomes, in contrast, very large and fast evolving. There it can take the capital role of bringing about inflation, as we shall see.
As previously mentioned, in addition to the traditional theoretical problems, other issues of more practical and mundane nature have been perturbing cosmologists in the last few years, which put the concordance ΛCDM against the wall. The practical problems are the presently irreconcilable observational differences between the concordance model predictions and a number of different kinds of cosmological observations. For example, those involving structure formation data (the so-called σ 8 tension), remain at a moderate level of 2 − 3σ (where σ 8 is the root mean square of fluctuations in matter density perturbations within spheres of radius 8h −1 Mpc); and, above all, the notorious conflict occurring at ∼ 5σ level between the local value of the Hubble parameter, H 0 , obtained from the traditional distant ladder techniques, and the value extracted from the early universe using CMB data. See e.g. [47,48] for a summary of these problems, and [49][50][51] for more comprehensive reviews and a long list of related references. For some people the disagreement in the latter case is sufficiently severe (and persistence over time) as to still be pretending that it can be attributed to a fluke, thus increasing the odds that its origin may come from physics beyond the ΛCDM. [52,53]. So, even if we remain agnostic, we cannot exclude that the prevailing model of cosmology might well be facing a crisis. Science, however, thrives on crisis since new ideas are then stimulated which could help to overcome it and maybe refine some aspects of the paradigm in force, or even originate a new one subsuming the old paradigm. Many proposals have indeed been made to alleviate these tensions, which include different forms of DE as well, despite that many of them are essentially ad hoc. As indicated above, clues to eventually substantiate the nature of the DE may come from a variety of cosmic and even astrophysical messengers [10]. For instance, the possibility of measuring the bending of light in the Solar System scale has been proposed [54]. But whatever the nature of the DE might be, we must provide an explanation for the role played by the vacuum energy in QFT. Indeed, in the absence of a correct understanding of the VED from first principles many DE proposals may look as an escape forward rather than a real alternative. This work, in contrast, intends to dwell further on the methods of QFT in curved spacetime so as to shed some useful light on these difficult problems. Recall that in QFT we treat the Λ term in the gravitational action as a formal quantity from which together with the ZPE (a computable quantity in curved spacetime) one can determine the renormalized VED, ρ vac , a fundamental concept in QFT. The physical cosmological term can then be defined from the VED as follows: Λ phys = 8πGρ vac . This is nevertheless not just the usual naive relation since G and ρ vac are now properly renormalized quantities in QFT in curved spacetime. In fact, in this work we report on progress made along the lines of the preceding comprehensive works [15][16][17], where a detailed account was made of the virtual contribution to the VED from quantized scalar matter fields. We found that these effects, when appropriately renormalized, translate into a (mild) dynamical evolution of the VED and G with the cosmic expansion, ρ vac = ρ vac (H) and G = G(ln H). This is not excluded by the cosmological principle, as it permits a homogeneous and isotropic dependence in time of physical quantities. More specifically, it was shown by explicit calculation and adequate renormalization that the VED behaves in the characteristic manner of the running vacuum model (RVM), see [41] for a recent comprehensive review (for a shorter summary, see e.g. [55]). Let us also note that a similar RVM interpretation of the vacuum energy density can be achieved in the string context, what has been called 'stringy-RVM' [56][57][58][59][60][61]. The fact that a QFT calculation and an effective string theory approach can lead to the same kind of RVM solution seems to indicate that the VED as a rigid concept is not very natural and that a dynamical evolution of the vacuum energy (density) should be more plausible. In effect, it has been recently shown that this can help significantly to improve the description of the overall cosmological data and in particular opens a viable solution to the well-known tensions afflicting the ΛCDM, see particularly the last phenomenological analysis [62] which was preceded by several other akin works, such as e.g. [63][64][65][66][67][68][69][70]. It is also interesting to remark that the RVM structure of the vacuum energy has been successfully tested against competing models (e.g. ghost models and holographic models of the DE) using cosmographical methods, which are essentially model-independent -see e.g. [71][72][73] for details. The model has indeed passed a battery of different tests [74,75] and the outcome is that the quality fit provided by the overall cosmological data is comparable, actually better, than that of the ΛCDM, if we attend to the verdict of the standard information criteria [62].
In this paper we continue the task of computing the dynamics of ρ vac (H) induced by the quantum effects of the quantized matter fields in Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime, which two of us (CMP and JSP) initiated in previous works [15][16][17], see also [41] for a comprehensive review. The result of the present, more complete, calculation (involving for the first time the spin-1/2 fermionic contributions) reconfirms that the combined dynamics of the vacuum adopts the RVM form indicated below. In these works, a new (off-shell) implementation of the adiabatic regularization prescription (ARP) was used to compute the renormalized ρ vac for a nonminimally coupled real scalar field. The method is based on a series expansion in the number of derivatives of the scale factor which introduces a hierarchy in some physical quantities evolving in a dynamical background [1][2][3]. Not only it was shown that the running of the VED was free from dangerous large contributions proportional to m 4 (quartic powers of the mass of the particles), but ρ vac was shown to be mildly evolving with the Hubble rate and hence with the cosmic expansion. As a result, if t 1 and t 2 are two particular values of the cosmic time, both close to the present, the corresponding values of ρ vac are connected through the approximate relation where H 1 ≡ H(t 1 ) and H 2 ≡ H(t 2 ) are the values of the Hubble function at times t 1 and t 2 , respectively and |ν eff | ≪ 1 is a small parameter. While the above relation is relevant for the (very mild) evolution of the VED in the current universe, the corresponding analysis of the early universe leads to a new mechanism of inflation called 'RVM-inflation', which relies on the existence of quantum effects of 6th adiabatic order, i.e. up to terms ∼ O(H 6 ) which have been first accounted for in the case of scalar fields in [16]. Despite of the fact that the family of Running Vacuum Models (RVM) has been in the literature for quite some time (cf. [40,41] and [76][77][78], and references therein), a full-fledged account based on QFT principles is much more recent [15][16][17]. This work is tightly related to the preceding studies, in which the adiabatic regularization was applied to the 'simple' case of one real scalar field. It was natural to perform the next step and check if spin-1/2 fermions do preserve the main conclusions derived from scalars, above all to verify if the corresponding vacuum fluctuations induce also a running of the VED independent of the quartic powers of their masses and hence remain also free from the traditional fine-tuning illness. So the main goal of this paper is to extend the computations done for the scalar field, by considering the quantization of spin-1/2 fermions in the FLRW background. The extension proves rather non-trivial since we are in curved spacetime and the computations with fermions are no less involved owing to the Fermi-Dirac statistics and the peculiarities inherent to the spinor calculus. It is however reassuring to find that despite the many new technicalities involved in the computation they do not alter the main conclusions derived from the calculation with scalars. The RVM form (1) at low energies is once more attained, but the contribution to the coefficient ν eff is, of course, different and involves non-trivial computational details. Similarly, we compute the fermionic contribution to the ∼ O(H 6 ) terms which are involved in the RVM mechanism of inflation occurring in the very early universe. The final results concerning the renormalized VED can be obtained by considering the contributions from an arbitrary number of quantized scalar and fermion fields. We combine the two types of effects and present a final formula which we refer to as the bosonic and fermionic contributions to the VED, with the understanding that additional effects from gauge fields and their interactions with matter would be necessary in our calculation in a more complete approach. It is nevertheless not necessary in our study since our fields are free, except for the non-minimal coupling of the quantized scalar fields with the external (non-quantized) gravitational field and the necessary spinorial affine connection of the fermion fields. All that said, the computation of the free field contribution from bosons and fermions in curved spacetime is already a rather formidable task. So, for the sake of a stepwise and clearer presentation, we will address the fermionic contributions here on equal footing to the presentation of the scalar part performed in the preceding studies [15,16].
This work is structured as follows: In Sec. 2, we consider a quantized scalar field φ nonminimally coupled to curvature and review the computation of its energy-momentum tensor (EMT) and corresponding vacuum expectation value (VEV) induced by the vacuum fluctuations of that field. The original study of this part was made in [15,16]. The 00th component of the VEV of the EMT constitutes the ZPE of φ. We define also its associated vacuum energy density (VED), ρ vac , and vacuum pressure, P vac . All these quantities are unrenormalized at this point and hence UV-divergent. In the same section we review the off-shell adiabatic renormalization of the VED introduced in the previous references, which involves as a distinctive feature the subtraction of the on-shell EMT at an off-shell renormalization point M up to 4th adiabatic order (in 4 spacetime dimensions). The main subject of this work is to address the corresponding calculation for spin-1/2 fermions and combine it with the scalar field case. In Sec. 3, we review the quantization of a Dirac fermion in a curved background, the corresponding Dirac equation and its spinor solutions obtained from adiabatic expansion of the field modes. The computation of the EMT and of its VEV for the case of a free quantized fermionic field in a spatially flat FLRW background is performed in Sec. 4. The off-shell adiabatic renormalization of the EMT for spin-1/2 fermions is addressed and we extract the renormalized ZPE, VED (ρ vac ) and vacuum pressure (P vac ) in this context. Additionally, some remarks on the trace anomaly and its role in our approach are discussed. Sec. 5 contains the combined results from all the quantized matter fields. Specifically, we compute the renormalized VED for a system made of an arbitrary number of quantized scalar fields non-minimally coupled to curvature (with different masses and non-minimal couplings) and an arbitrary number of quantized spin-1/2 free fermion fields. In the same section we report on the corresponding running of the gravitational coupling G = G(H), which goes hand in hand with the running of ρ vac (H) in order to preserve the Bianchi identity. We also discuss the mechanism of 'RVM-inflation' with the combined contribution from all these fields, and derive the equation of state (EoS) of the quantum vacuum for that system of quantized bosons and fermions fluctuating in it. Remarkably, the vacuum EoS is no longer equal to w vac = −1, the reason being that the vacuum pressure and the VED are not exactly related in the usual manner (viz. P vac = −ρ vac ) since P vac and ρ vac are independent functions of the Hubble rate H and its time derivatives owing to the quantum effects. In the current universe, there is still some remnant of these quantum effects which induce a small (but potentially measurable) departure making the quantum vacuum mimic quintessence. The conclusions are delivered in Sec. 6 together with some additional discussion. Finally, three appendices are included. In Appendix A, we define our conventions and some useful formulas. The last two appendices, B and C, are rather bulky since they collect a number of cumbersome expressions related to the adiabatic expansion of the VEV of the EMT and the Fourier modes of the fermionic field (computed up to 6th order for the first time in the literature).

Vacuum energy density of a non-minimally coupled scalar field
In this section, we summarize the results for the vacuum energy density and pressure associated with a quantized scalar field in FLRW spacetime as obtained in [15,16]. In passing we introduce some notation which will be useful also for the fermionic calculation that will be subsequently reported. The Einstein-Hilbert (EH) action for gravity plus matter reads The term ρ Λ has dimension of energy density and sometimes is called the vacuum energy density, but this is inaccurate in the formal QFT context since renormalization is necessary and the physical vacuum energy density, ρ vac , is not just that term. In fact, ρ Λ is at this point just a bare parameter of the action, as the gravitational coupling G itself. Varying the action with respect to the metric provides Einstein's equations 1 8πG with G µν = R µν − (1/2)g µν R the usual Einstein tensor and T m µν the EMT of matter 1 : The matter action S m may contain a variety of contributions, including those from incoherent matter, but it will be enough to focus on fundamental effects from quantized scalar and fermion fields. Here we shall compute the fermionic contribution to the VED. But let us summarize first the situation with the scalar field part. The latter was dealt with in great detail in the two previous studies [15,16], in which the VED calculation was addressed under the assumption that no effective potential was present. However, we admitted a non-minimal coupling of the scalar field to gravity. That calculation in curved spacetime is already sufficiently demanding and in addition it furnishes the universal part of the VED through the zero-point energy (ZPE) effects in the curved background, see next section. The classical action for a non-minimally coupled real scalar field is the following: where ξ is the non-minimal coupling with gravity. It is well known that this action enjoys (local) conformal symmetry in the massless case with ξ = 1/6. However, the value of ξ is not fixed in our computation and in general we do not assume the presence of such a symmetry. Varying the above action with respect to the scalar field leads to the Klein-Gordon (KG) equation with non-minimal coupling: . The corresponding EMT of φ follows from the metric variation of the action (5) according to the recipe (4), and yields As indicated, we perform the calculation in cosmological (FLRW) spacetime with flat threedimensional metric. For convenience we use the conformal frame ds 2 = a 2 (τ )η µν dx µ dx ν , with η µν = diag(−1, +1, +1, +1) the Minkowski metric in our conventions (cf. Appendix A), a(τ ) is the scale factor and τ the conformal time. Differentiation with respect to τ will be denoted with a prime, so for example H ≡ a ′ /a is the corresponding Hubble function in conformal time. We will perform the explicit calculations using the conformal metric but our final results will eventually be rendered in terms of the usual Hubble function H(t) =ȧ/a in cosmic time t (where a dot denotes differentiation with respect to t). Recall that dτ = dt/a and hence H = aH.
If we switch on the quantum fluctuations of the φ field it is natural to consider the following decomposition: in which the background φ b and the fluctuating part δφ are understood to be independent. In particular, this is also the case for the corresponding Fourier decomposition in frequency modes. The decomposition of the fluctuating part reads with the usual commutation relations for the creation and annihilation operators, A k and A † k : By using these relations, the KG-equation (6) in terms of the frequency modes can be put as where Ω 2 k ≡ k 2 + a 2 m 2 φ + (ξ − 1/6) R and we recall that () ′ ≡ d/dτ (). The above differential equation does not possess a close analytic solution for the entire cosmological evolution. However, it can be solved by means of what is called an adiabatic series expansion, which is essentially a WKB-type solution. First of all, it is necessary to introduce the following ansatz for the mode functions: Notice that the modes are normalized through the Wronskian condition which is essential to preserve the canonical commutation relations for the quantized field φ. By introducing the above ansatz into (11), the function W k (effective frequency) is the solution of the (WKB-type) non-linear differential equation

Zero-point energy and adiabatic expansion
For a slowly varying effective frequency Ω k (τ ) one can proceed to solve this equation perturbatively with the help of an asymptotic series which can be organized through adiabatic orders. This approach constitutes the basis for the aforementioned ARP [79][80][81][82][83][84][85][86][87], see also [88][89][90][91][92][93] and [94][95][96][97][98] for more recent applications and extensions, and the textbooks [1][2][3] for a more systematic presentation. The quantities k 2 and a are taken to be of adiabatic order 0. Of adiabatic order 1 are: a ′ and H. The quantities a ′′ , a ′2 , H ′ and H 2 and linear combinations are taken of adiabatic order 2. It can be summarized by saying that each time derivative increases one unit the adiabatic order. So, the ansatz solution for W k can be written as an adiabatic series expansion: in which the superscript indicates the adiabatic order. We note that only even orders are allowed, which is justified from the general covariance of the result since only tensors of even adiabatic order can be present in the effective action and the field equations. Explicit calculation indeed corroborates the absence of the odd adiabatic orders. The "seed" to initiate the adiabatic series (i.e. the zeroth order contribution to W k ) is where M is an arbitrary off-shell scale. Nothing enforces us to take the mass of the particle at this point, we only need to preserve the adiabaticity of the expansion. The floating quantity M will play the role of renormalization scale, as it will be seen. In fact, as it was shown in [16], this parameter can also be used as the renormalization scale in the DeWitt-Schwinger expansion [4] of the vacuum effective action W eff [1], to wit: the effective action obtained from integrating out the vacuum fluctuations of the quantized matter fields. From the explicit expression of W eff one can also derive the VEV of the EMT -denoted T δφ µν -by computing its metric functional derivative as follows: This formula is of course similar to Eq. (4), but for the vacuum effective action. This alternative method provides exactly the same result as the WKB expansion of the field modes, as outlined below, and it was illustrated in great detail in [16] for the case of the quantized scalar fields. Such a cross-check in the determination of the VEV of the EMT involves a significant amount of calculations and provides a nontrivial validation of the entire renormalization procedure. The same holds good for the case of fermions but we shall not present the details of the DeWitt-Schwinger method for fermions here.
Next we summarize the mode expansion for scalar fields. By introducing W (0) k given above in the r.h.s. of (14), the terms of adiabatic order 2 can be collected to find the next term in the series W (2) k , with the result Here ∆ 2 ≡ m 2 −M 2 is the difference between the quadratic mass of the field and that of the off-shell scale, and is of adiabatic order 2 because it is necessary for renormalization. Loosely speaking, since M 2 and ∆ 2 appear together in the expansion they need to be of different adiabatic order so as to obtain a consistent adiabatic expansion exploring the off-shell regime. Since M 2 is of adiabatic order 0, the next-to-leading order for ∆ 2 to be compatible with general covariance is precisely order 2. This fact is reconfirmed on using the aforementioned DeWitt-Schwinger expansion [16]. Introducing W (2) k on the r.h.s. of (14) we can obtain W (4) k , . . . and so on. The higher order adiabatic terms become progressively more and more cumbersome since the number of terms involved in the expansion becomes larger and larger. In our case we reach up to order 6, i.e. we compute the series up to W (6) k . This was done for the first time in the literature in [16] for scalar fields, and we will also be done here to order W (6) k for the first time for fermions, see Sec. 4. Notice that attaining the order W (6) k is indispensable in order to study RVM-inflation in the early universe [16]. Extensive use of Mathematica [99] has been made to handle these bulky calculations.
Once obtained the expansion of W k we can compute the mode functions h k and other physical quantities such as the EMT to the suitable order, in particular the EMT trace, which can be used to compute the vacuum pressure (see below). The technical details for the scalars are not going to be repeated here 2 , but it is important to remark that these quantities present divergent terms up to 4th adiabatic order (in 4-dimensional spacetime). The approach we adopt here is the same as that which was proposed and amply tested in [15,16], namely we define the renormalized vacuum expectation value (VEV) of the EMT by taking the on-shell value (at an arbitrary adiabatic order) and subtracting from it the divergent orders at an arbitrary scale, which we denote M : Here the superscript (0 − 4) refers to the UV-divergent subtracted orders, i.e. from 0th up to 4th adiabatic order, all of them being UV-divergent (the higher adiabatic orders being all finite in n = 4 spacetime). Notice that for M = m φ the above definition provides the natural generalization of the subtraction of divergent constants performed to obtain finite results on trivial backgrounds (such as Minkowski spacetime). However, in curved backgrounds the mode by mode subtraction implied in the above prescription is not just a constant term; and moreover for arbitrary M the corresponding renormalized result allows us to test the evolution of the VED with the scale M . As previously noted, this feature obviously offers a floating scale which is characteristic of the renormalization group (RG) analysis in cosmology [40,41]. Let us clarify, however, that we distinguish M from the 't Hooft's mass unit µ in dimensional regularization (DR), which will not be used in this work at any point, although it can be invoked as an intermediate regularization procedure (not at all for renormalization though) if one likes [15,16]. The parameter µ is unphysical and is used in the minimal subtraction scheme (MS) with DR to define the renormalization point [100]. We should emphasize that we do not use MS at all in the present work, although one could make (optional) use of DR in intermediate steps. In these cases, the quantity µ always cancels out and the final renormalized expressions depend on M only, as it is the case e.g. of the effective action of vacuum. But the full effective action (which involves the classical and quantum parts) is of course independent of M as well, as the running of the couplings exactly compensates for the explicit M -dependence of the quantum effects. This is, of course, the standard lore of the RG program, see [16] for detailed considerations along these lines and making use of the effective action.
We refrain from writing out the unrenormalized expression for the EMT in the case of scalar fields, see [15,16] for full details. Let us however quote the renormalized result emerging from the ARP prescription (19). Expressing the final result in terms of the cosmic time and the corresponding Hubble function H =ȧ/a, we find for the ZPE the following result: We used the notation O(H 6 /m 2 φ ) to collectively refer to the terms of adiabatic order 6 (consisting of 6 time derivatives of the scale factor). It may include terms such as H 6 /m 2 φ , but also many other combinations such asḢ 3 /m 2 φ , H 2 ... H/m 2 φ , . . . There are actually many terms of this sort and they have been reported explicitly in [16]. We refrain from writing them down again here and invite the reader to check the aforementioned paper for more details. We will report explicitly on the 6th-order adiabatic terms only in the case of fermions (cf. Sec. 4) since these are computed for the first time in this work.

Renormalized vacuum energy and vacuum pressure
We are now ready to compute the vacuum EMT, T vac µν , which will lead us to the VED, ρ vac , and vacuum's pressure, P vac . As in [16], we write the vacuum EMT as the sum of the renormalized parameter ρ Λ and the renormalized VEV of the EMT, which embodies the finite form of the adiabatically renormalized vacuum fluctuations: Since the vacuum is expected to be a most symmetric state free of any new parameter, this expression must take on the form of a perfect fluid: T vac µν = P vac g µν + (ρ vac + P vac ) u µ u ν , where u µ is the 4-velocity (u µ u µ = −1). In conformal coordinates, in the comoving (FLRW) frame, u µ = (1/a, 0, 0, 0) and u µ = (−a, 0, 0, 0). Taking the 00th and iith-component (any i = 1, 2, 3 is good owing to isotropy, so we take i = 1) one finds the precise form of the vacuum energy density and pressure [16]: where isotropy allows to express the result in terms of the trace T δφ of the fluctuating part, if desired. Notice that ρ Λ (M ) in the above expressions is the renormalized form of the corresponding bare parameter appearing in the EH action (2) and it has units of energy density. The VED, however, is not just this renormalized parameter but the renormalized sum (23). Although is tantalizing to call ρ Λ (M ) "the CC density", and in fact this has been common in the literature (especially when the discussion is strictly classical without considering quantum effects), this is not strictly correct since the physical CC is not simply 8πGρ Λ but 8πGρ vac , that is, the physical vacuum energy density is connected with the physical Λ through ρ vac = Λ/(8πG). The parameter Λ which is measured in the observations is indeed defined through this expression, which is precisely computable in QFT from Eq. (23). We shall show once more for fermions (as we did for scalar fields in the previous works [15,16]), that the adiabatically renormalized form of the running VED is free from the huge ∼ m 4 contributions that are usually attributed to the VED in other (inappropriate) renormalization schemes, and therefore the renormalized expression that we will obtain can be perfectly consistent with the measured Λ. To be sure, it is not our aim to predict this value but rather to show that the theoretical formula points naturally to a value as small (in natural units) as measured by the observations. A simple way to condense these ideas is to say that the VED is related with the ZPE and ρ Λ as follows: "VED = ρ Λ + ZPE", i.e. Eq. (23). Parameter ρ Λ is initially just a bare coupling in the effective action and it has no direct phenomenological interpretation, not even after renormalization. On the other hand, the ZPE embodies the quantum fluctuations of the massive fields and calls also for renormalization since it is originally UV-divergent. The physical VED in this context is then the renormalized sum of these two contributions, and it can not be split apart since the separate terms make no sense in an isolated way. Observations are sensitive only to the sum. Furthermore, as we shall see explicitly for the fermionic case, there is a crucial cancellation of the quartic mass terms when we consider the evolution of the sum ρ Λ + ZPE as a function of the renormalization point, which does not occur if the two terms are dealt with separately. This was already pinpointed for the case of scalar fields in [15][16][17] .
With these provisos, the expression for the VED of the scalar field can be obtained. Notwithstanding, the final renormalized result still requires a physical interpretation since it depends on the renormalization scale M . In fact, recall that the result depends on both the values of M and H (and corresponding time derivatives), which are independent arguments. The scale M can not be left arbitrary at this point since we wish to provide an estimate of the VED at a given expansion epoch. As previously indicated, the vacuum effective action W eff is explicitly dependent on M despite the full effective action is of course RG-independent. Thus, following , [15,16] an adequate choice of the renormalization point M is to select it equal to the value of H at the epoch under consideration. This corresponds to choose the RG scale around the characteristic energy scale of FLRW spacetime at any given moment, and hence it should have physical significance. In actual fact this is in analogy with the standard practice in ordinary gauge theories, where the choice of the renormalization group scale is made near the typical energy of the process. In what follows we derive the 'low energy' form of the VED along these lines. This is actually the form that applies for the current universe. Subsequently we will focus on the running gravitational coupling G(M ) and its relation with the running ρ vac (M ).
We should also point out that the reach of our considerations concerns the calculation of the evolution (or 'running') of the VED only, rather than predicting its current value. Given that value, however, we can predict how it evolves with H around our epoch, or any other epoch. Now in the absence of an observational input at some expansion epoch H(t) we cannot compute ρ vac (H) at other values of H (i.e. at other epochs of the cosmic evolution). To compute the value of the VED at present is out of the scope of the renormalization program since the latter is based on the RG flow, which requires a boundary condition. This is exactly the same situation as in any renormalization calculation, we need the input values of the parameters at one scale to predict some observable (e.g. a cross-section) at another scale. The truly relevant feature of our calculational approach, as it should be clear at this point, is that the RG-flow is completely smooth. It only depends on the evolution of H and is completely free from spurious effects associated with the large contributions from the quartic masses of the fields. These quartic terms are the traditional kind of undesirable effects which spoil the physical interpretation of the renormalization program concerning the CC and the VED. They are typically found in calculations of the VED whose renormalization is based on, say, the MS scheme. Most existing approaches to the CC problem in the literature exhibit this unwanted feature, which is already at the basis of the Minkowskian calculation and is of course unacceptable for a realistic description of the VED in curved spacetime [41]. A similar situation is found in Schwarzschild and in de Sitter backgrounds, see e.g. [101,102].
Bearing in mind the above considerations, the final result for the running VED at low energies (specifically the part triggered by the quantized scalar fields), can be best written in terms of the evolution between two expansion history times. It is natural that we choose the current epoch (characterized by the value H 0 of the Hubble parameter) and relate it with the value of the VED at a nearby epoch H of the cosmic evolution 3 . The approximate final result can be rendered in the following very compact form [16]: where ρ 0 vac ≡ ρ vac (H 0 ) being the current value of the VED (accessible by observations) and H 0 today's value of the Hubble function. It is necessary to remark that ν eff is an effective parameter expected to be small due to its proportionality to m 2 φ /m 2 Pl . Remarkably, the above dynamical form of the VED turns out to adopt the canonical RVM form, see [40,41] and references therein. Phenomenological studies based on fitting the above RVM formula to the overall cosmological data indeed provide an estimate for ν eff at the level of ν eff ∼ 10 −3 and positive [62], see also [63][64][65][66][67][68]. This phenomenological determination turns out to lie in the ballpark estimate of the theoretical expectations [103]. The order of magnitude is reasonable if we take into account that the masses involved here pertain of course to the scale of a typical Grand Unified Theory (GUT) where, in addition, a large factor must be included to account for the large multiplicity of heavy particles. In Sec. 5 we provide a more general formula where an arbitrary number of species of bosons and fermion fields are included. It is worth noticing that the order of magnitude of ν eff picked out in the mentioned study is perfectly compatible with the result recently obtained from the Big Bang nucleosynthesis (BBN) bound in [104], although in the latter case the bound was not sensitive to the sign of ν eff .
On the other hand, the computation of the pressure in an analogous way (we refrain from providing more details on the scalar field contribution, see once more [16] and [17] for a fullfledged account) enables us to write an explicit expression for the equation of state (EoS) of the vacuum [17]. The leading expression for the current universe is the following: For very low redshift z and in terms of the current cosmological parameters Ω 0 This result is especially worth emphasizing since it predicts a small departure from -1 which could perhaps be measured around the present time. Recall that the traditional value of the EoS of the Cosmological Constant is just −1. The above result implies that the quantum vacuum receives small quantum effects which trigger a departure of its EoS from −1. For instance, if we adopt the positive sign for ν eff , as obtained from the latest fitting analysis to a large set of different kinds of observational data [62], then Eq. (28) predicts that the vacuum energy behaves as quintessence around the current time.
As noted, the EoS formula (28) is valid only for small values of the redshift z, but one can show that the departure is even bigger in the past, adopting a kind of chameleonic behavior by which the EoS of the quantum vacuum tracks the EoS of matter at high redshifts, see [17] and Sec. 5.4 for more details. All in all, these results have been predicted from first principles, namely from explicit QFT calculations in the FLRW background. In particular, the fact that the quantum vacuum may currently mimic quintessence is remarkable since the result does not rely on ad-hoc fields or on any other phenomenological ansatz.
3 Quantization of a spin-1/2 fermion field in curved spacetime As pointed out in the introduction, the main goal of this work is to extend the QFT results for the VED obtained for quantized scalar fields, which we have summarized in the previous section, to the case of quantized spin-1/2 Dirac fermion fields and then combine the two types of contributions in closed form. The calculation of the renormalized VED for free spin-1/2 fermions is also nontrivial and rather cumbersome, and requires a devoted study, which we present here (see also the appendices provided at the end for bulky complementary details). While the QFT treatment is analogous to the case of scalars, the specific technicalities are quite different and no less intricate, but fortunately the final result proves to be in consonance with the one previously derived for the scalars, so it is perfectly possible to furnish a close form for the combined contribution to the VED involving an arbitrary number of non-interacting scalar and spin-1/2 fermion fields, cf. Sec. 5.
The study of the solutions of the Dirac equation in curved spacetime goes back to the works from many decades ago by Fock, Tetrode, Schrödinger, McVittie, Bargmann, Wheeler and others: see e.g. [105][106][107][108][109], where the relevant historical references are given and different aspects of spin-1/2 fermions in curved spacetime are studied, including a detailed account for the solutions in FLRW spacetime -see also the review [110], with a rather complete list of references. On the other hand, the subject of adiabatic regularization for fermions has been previously treated in the literature in different applications, see e.g. [84] as well as the more recent papers [90][91][92][93] where emphasis is made on exact solutions e.g. in de Sitter spacetime. The calculation of the renormalized VED in FLRW spacetime is, however, more complicated for it does not admit an exact solution. Our strategy to circumvent this problem is based on using an off-shell variant of the ARP framework [15,16] which leads to the RVM behavior of the vacuum energy [40,41]. The RVM framework has proven rather successful in mitigating the cosmological tensions [49][50][51], as shown in different phenomenological analyses, such as [62] and [69,70]. On the theoretical side, attempts at computing the VED with other procedures has led to the traditional calamity with the quartic powers of the masses. Here we will show that using the off-shell ARP to tackle the VED contribution from fermions generates a result which is free from these difficulties and fully along the lines of what has been obtained for the scalar fields in the previous sections and originally in [15,16]. Therefore, the combined contribution from fermions and scalar fields to the VED is compatible with a smooth running of the cosmological vacuum energy and is consistent with the aforementioned phenomenological analysis of the RVM as a possible solution to the cosmological tensions.
Since it will be necessary a considerable amount of formalism to treat fermions within the adiabatic approach, it is convenient to summarize first the necessary aspects of that formalism before we can put forward our main results concerning their contribution to the vacuum energy density. It will be useful to fix some notation as well. Once more we perform the calculations in FLRW spacetime with flat three-dimensional metric. Consider a free Dirac spin-1/2 field, described by the four-component spinor ψ. In our conventions, the Dirac action in curved spacetime is given by In the above expression, m ψ denotes the mass of the Dirac field andψ ≡ ψ † γ 0 the adjoint spinor.
Since we are in a curved background, the partial derivative of a spinor ∂ µ ψ has been replaced with the corresponding covariant derivative ∇ µ ψ, which is defined below. Moreover, gamma matrices in curved spacetime are also needed, they are sometimes indicated (as above) with an underline to distinguish them from the Minkowski space gamma matrices. The former are γ µ (x) (which are generally functions of the coordinates) whereas the latter are the constant matrices γ α in flat spacetime. As it is well-known, to obtain a representation for the curved spacetime gamma matrices in terms of the Minkowskian gamma matrices we need to introduce the local tetrad or vierbein field (in 4-dimensional spacetime) e µ α . It is defined in each tangent space of the spacetime manifold and relates the curved spacetime metric with the Minkowskian one in the usual way: where η αβ is the Lorentz metric in the local inertial frame specified by the normal coordinates at the given spacetime point. The general relation between the two sorts of gamma matrices is γ µ (x) = e µ α (x)γ α . Specifically, in a spatially flat FLRW spacetime the vierbein in conformal coordinates is e µ α = diag (1/a(τ ), 1/a(τ ), 1/a(τ ), 1/a(τ )) where a(τ ) is the scale factor as a function of the conformal time. Whence the gamma matrices in this background are time-dependent and related to the constant flat spacetime ones as follows: γ µ (τ ) = γ µ /a(τ ). This relation insures that they satisfy the following anti-commutation relations: provided, of course, the (constant) flat space gamma matrices satisfy γ α , γ β = −2η αβ I 4 . In order to obtain the equation of motion, i.e. the covariant Dirac equation in curved spacetime, one has to vary the covariant action (29) with respect to the spinor field, giving The covariant derivative is defined through the spin connection, ∇ µ ≡ ∂ µ − Γ µ . The spinorial affine connection Γ µ satisfies the equation [106] Γ where Γ µ νρ are the Christoffel symbols. The above equation is tantamount to require the vanishing of the covariant derivative of the curved space gamma matrices: ∇ ν γ µ (x) = 0 [2], i.e. the curved-space gamma matrices are defined to be covariantly constant over the spacetime manifold. Using the Christoffel symbols in the conformally flat FLRW metric as given in Appendix A, the explicit solution of Eq. (32) can be found, with the following result: This expression can then be inserted in Eq. (31).
In this way we have obtained an explicit form for the Dirac equation in FLRW spacetime with spatially flat metric. We are now in position to attempt a solution by expanding the quantized fermion field in mode functions: Here B k,λ and D † k,λ are creation and annihilation operators which satisfy the standard anticommutation relations, The momentum expansion of the mode functions u k,λ and their charge conjugates v k,λ can be conveniently written in terms of two 2-component spinors ξ λ ( k) and corresponding spinor modes h I k and h II k : Using this representation, Eq. (31) splits into a system of two coupled first order equations for each of the two types of spinor modes h I k and h II k : After straightforward calculation, these equations can be rewritten as two second order decoupled equations: where with The fact that (38) only depends on the modulus of the momentum, k, justifies the notation used for the modes h I k , h II k , with no arrows. Following the same prescription as in the case of scalar fields (cf. Sec. 2), we have introduced an off-shell scale M , which again will take the role of renormalization scale. Correspondingly, we have defined ∆ 2 ≡ m 2 ψ − M 2 and once more assigned adiabaticity order 2 to it. We did not change the notation ∆ as compared to the scalar case since the final formulas do not depend on ∆ but on M and the respective physical masses. The argument of ω k will be omitted from now on, unless it takes a different value from M . The normalization conditions for the mode functions involved in ψ are implemented through the Dirac scalar product: and similarly for As mentioned in the previous section, the number of time derivatives of the cosmological scale factor a(τ ) that appear in a term of the expansion is called adiabatic order of the term. In order to solve the differential equations (38) we may follow a recursive process which preserves the adiabatic hierarchy, just as we did with the scalar fields. Let us first redefine h I k and the time variable as follow Substituting these relations into the equation for h I k in (38) we find Since ǫ 2 includes two derivatives, it contains terms of second and higher adiabatic order. We can ignore it to find the leading order solution so that we get a first approximation Notice that h I k,1 formally satisfies a differential equation with the same form as (38) for h I k . So that, we can repeat the process: The corresponding differential equation for h I k,2 is Once again, ǫ 4 consists of terms of adiabatic order 4 and higher. We can approximate a solution of (48) by neglecting ǫ 4 : whereby the approximation to h I k can be further improved: By iterating the procedure, we can obtain a better and better approximation to h I k , and after ℓ > 1 steps we find where, for ℓ ≥ 1, Now that the general method has been set up, let's find the 0th order solution for h I k . From (46), the most generic solution for h I k is where the time independent function f (0) k (of adiabatic order 0) accounts for the integration 'constant' (strictly speaking, a function of the momentum but not of conformal time) in the exponential. As for h II k , by comparing both lines of (38) it is clear that it is possible to proceed in an analogous manner. So we obtain where g (0) k has the same paper as f (0) k . To find the zeroth adiabatic order it is just enough to expand this solution and keep zero order terms. However, some extra caution is needed when dealing with the integrand in the exponential of (53), which may be expanded up to 1st order as where The reason is that the integration of the second term in the exponential factor is: so it yields a real term of adiabatic order zero, meaning that the expansion of Ω k up to 1st order in the integral was mandatory. We have not included an explicit multiplicative factor related with the constant of integration 4 since it is already represented by f such that the above solution can be compatible with mode functions in Minkowskian spacetime, so we can write Next we move on to the solution at 1st adiabatic order. As we have mentioned, the quantity ǫ 2 defined in (44), contains terms of adiabatic order two and higher, so it is not necessary to find the first order solution. It is enough to find the first order term from the denominator of (53). So, Similarly for the second spinor mode h II k : where f (1) k and g (1) k come from integration constants, as mentioned in the footnote of the previous page. By imposing the normalization condition (42), which has to be satisfied at each adiabatic order, it is possible to see that these constants are purely imaginary, that is To continue, we deal with the 2nd adiabatic order of the mode functions, i.e. h I,II(2) k . At this time, we have to include Ω 2 k,1 = 1 + ǫ 2 in our considerations (this term contains 2nd order adiabatic terms and beyond). Starting from Eq. (50), we have where ǫ 2 can be computed to be With this result, it is immediate to obtain an approximation for Ω k,1 valid up to the third adiabatic order: On the other hand, an expansion of the product Ω k Ω k,1 is necessary to improve the approximation of h I,II k , as one can see from equation (50). As earlier, if we wish to present a second order approximation of the modes we have to expand that product up to 3rd adiabatic order in the exponential. The expansion can be presented as follows: where the dots represent the contributions of adiabatic order bigger than 3, and the indicated terms in the expansion read As noted before, ω (2) k are real. Again, when integrated inside the exponential of equation (50) the former two give a real contribution, whereas the latter two become part of the phase of the mode and play the role of an effective frequency: The last result holds good up to an arbitrary function of momentum (constant in conformal time) multiplying the whole result. We account for this arbitrary constant by introducing the functions f k , . . . at each order. An efficient strategy to compute the integrals involved in the above calculation (and many other ones of a similar sort, see Appendix B for a sample of them) is to set up an ansatz which respects the adiabaticity order of the calculation. The ansatz consists of a finite number of terms (in fact, a linear combination of them) taken each at the given adiabatic order and with coefficients (or 'form factors') which must be determined. The terms of the ansatz are constructed out of the derivatives of the scale factor and the parameter ∆ 2 (which we recall is of second adiabatic order).
For instance, in order to compute the integral of ω (3) k in Eq. (66), we know that the result must be of second adiabatic order. Hence as a suitable ansatz we use a linear combination of second order adiabatic terms: where again the term 'const.' at the end means that it does not depend on the integration variable, τ . By taking derivatives with respect to (conformal) time of the last expression and comparing with ω (65) and (62), the expansion of h I k up to 2nd order is In a similar way, The normalization condition fixes the following relations: So far, the expansion for the modes h I k and h II k up to 2nd order has been presented. One can continue with the procedure formerly described to reach higher orders, although of course the calculation becomes more and more involved. We should keep in mind, though, that the adiabatic expansion is an asymptotic expansion. While for renormalization purposes it is enough to stop the expansion at 4th adiabatic order (in 4-dimensional spacetime), it is nonetheless necessary to reach up to 6th order to meet the finite terms ∼ H 6 that are dominant in the early universe and capable of triggering inflation in this framework (cf. Sect. 5.3) 5 . We shall refrain from presenting these cumbersome formulas in the main text, see Appendix B.
It is worth noticing that there is some residual freedom in the previous calculations since, we can not determine entirely the set of integration constants that appear during the calculations f k , . . . Because of the normalization condition (42) of the mode functions, some restrictions such as (61) and (71) apply. Fortunately, as commented in more detail in Appendix B, the satisfaction of these restrictions is enough for the observables to be independent of this residual freedom. So that, is enough to set all of them to 0 to get, for instance, the desired values of the energy density and pressure.

ZPE and VED for a spin-1/2 field in FLRW spacetime
The computation of the Fourier modes for a quantized fermion field through adiabatic expansion as explained in the previous section is just the first step to compute the vacuum energy density (VED). The next step towards the VED is to obtain the ZPE associated with Dirac fermions in curved spacetime. As it well known, traditional computations of ZPE suffer from the wellknown headache of carrying highly unacceptable contributions proportional to the quartic powers of the masses, ∼ m 4 . This is so both for scalar and fermion fields, and it is already the case in flat, Minkowskian, spacetime, see e.g. [40,41] for a detailed discussion and more references. In curved spacetime we have the same situation, in principle, but in addition we encounter subleading, curvature dependent, contributions which do not exist in the flat case, as we shall see in a moment. To handle this issue, an adequate renormalization prescription is called for.
The calculation of the ZPE performed here for spin-1/2 fermions is closely related with the one previously put forward for scalar fields in [15,16] and summarized in Sec. 2. Once more the computation will be done through adiabatic expansion of the field modes and will be carried out up to 6th adiabatic order, since this is the first non-vanishing order on-shell, i.e. when fixing the renormalization scale M to the value of the mass of the fermion m ψ . However, the off-shell computation at 4th order is already very useful as a means to determine the RG running of the VED as a function of the scale M . This is actually one of the main new features of the off-shell ARP method proposed in [15,16], which leads to the cosmic evolution of the VED. Next we consider the actual calculation for spinor fields.
To find out the ZPE, we start from the definition of EMT in Eq. (4). In this case we have to evaluate the functional derivative applied to the fermion action (29). Upon a straightforward calculation we arrive at the following symmetric expression: in which the equation of motion (31) and its hermitian conjugate have been used. We treat this spinor field as a field operator and upon using its expansion in Fourier modes and utilizing the anticommuting algebra of the creation and annihilation operators, Eq. (34), we can compute the VEV of the various components, which reflect the contribution from the vacuum fluctuations of the quantized fermion fields. The method is the same as for the scalar fields [15,16], but details are of course different and shall be omitted here. After significant work, we find that the VEV of the 00th component of the EMT can be cast as follows: where ρ k is a function of the previously defined mode functions (which can be computed through adiabatic expansion): The explicit form of the adiabatic expansion of ρ k is rather cumbersome; the reader may find the final result for T δψ 00 in the Appendix B. Let us note that for off-shell renormalization at a point M it suffices to adiabatically expand the solution up to 4th order (as it was prescribed in Eq. (19)), see Eq. (79) below. However, we will provide the result up to 6th order so as to be sensitive to the on-shell result (occurring for M = m ψ ) and also because it is important for the inflationary mechanism in the early universe (cf. Sec. 5.3). Renormalization of the above expressions is indeed necessary since the VEV of the EMT is formally divergent. The UV-divergent contributions appear up to 4th adiabatic order (in n = 4 spacetime dimensions), so that one has to subtract terms up to this order to obtain a finite result.

Divergence balance between bosons and fermions in vacuum
The unrenormalized VEV of the EMT can be split into two different parts, divergent (in the UV sense) and non-divergent. Explicit calculation using the formulas of Appendix C) shows that the divergent part reads as follows: As it is easy to see, there are terms diverging quartically, quadratically and logarithmically. The non-divergent part contains the remaining terms, all of them being finite. The above ZPE is, as warned, an unrenormalized result at this point. However, before we proceed to renormalize that expression in the next section, it may be instructive to check if there is a chance for a cancellation between UV-divergent terms between fermions and bosons in the supersymmetric (SUSY) limit, if only for the leading divergences. In the on-shell case (M = m and hence ∆ 2 = 0) the above equation (76) simplifies to It coincides with the Minkowskian result for a = 1 (since H = 0). Now, in a SUSY theory, in which the number of boson and fermion degrees of freedom (d.o.f.) is perfectly balanced, we should expect that the leading (quartic) divergences cancel among the fermionic and bosonic contributions in the vacuum state [111,112] since in such case the scalar and fermionic fields have the same mass m. Thus the quartically divergent contribution from the first term of (77) should be minus four times the corresponding result for one real scalar field 6 . We can check it is indeed so using the above formulas, for in the on-shell limit and projecting the UV-divergent terms of the first two adiabatic orders only, we find that the contribution from one real scalar field in FLRW spacetime with spatially flat metric is [16] T δφ We confirm that the first term (the quartically divergent one) of this expression is of opposite sign to the first one in (77) (77) it is clear that for fermions we only have subleading divergences of logarithmic type, which also hinge on curvature effects since they are again proportional to H 2 and would also vanish in Minkowski space. Hence there is no possible cancellation of these subleading divergences between bosonic and fermionic d.o.f., in FLRW spacetime, even in the exact SUSY limit. Of course, our framework is not placed in the context of supersymmetry, but it serves as a consistency check of our calculations. See also the discussion in [113,114].
Although it is possible to introduce a cutoff for a preliminary treatment of the subleading divergences (and maybe to speculate on its possible meaning) it is not really necessary. One simply has to implement appropriate renormalization since renormalization is anyway necessary to deal meaningfully with the VED, as there is no way to cure the divergences from the combined contributions from bosons and fermions and it is not useful to be left with a "physical" cutoff. Dealing with a cutoff is always ambiguous as it is generally not a covariant quantity. Renormalization gets rid of cutoffs and one can preserve covariance, which is safer for a physical interpretation of the final results. The adiabatic renormalization is ideal in this sense since the adiabatic expansion generates automatically a covariant result.
It is well-known that the renormalization program in QFT requires the presence of a renormalization point, as well as a renormalization prescription. The renormalization point is a floating scale characteristic of the RG. As in the ordinary adiabatic procedure, to implement the renormalization of the EMT in 4 spacetime dimensions we perform a subtraction of the first four adiabatic orders, which are the only ones that can be UV-divergent [1][2][3]. However, in contrast to the usual recipe, in which the subtraction is performed on the mass shell value m of the quantum field, we perform it at an arbitrary scale M since this enables us to explore the RG evolution of the VED and ultimately connect it with its cosmic evolution. This is the specific feature of the adiabatic renormalization procedure (ARP) for the VED that was proposed in [15,16] -see also [41] for additional details and a comparison with other renormalization schemes. The resulting renormalized VED ensuing from this procedure is free from the usual troubles with the quartic powers of the masses and their inherent fine tuning problems.
Finally, let us note that dealing with the CCP in Minkowski spacetime using, for instance, the MS scheme and assigning some value to the 't Hooft's mass unit µ in DR (as discussed so many times in the literature), is entirely meaningless. It is not only devoid of meaning in that a non-vanishing cosmological constant cannot be defined in Minkowski space without manifestly violating Einstein's equations; it is meaningless also on account of the fact that there is no sense in associating the scale µ with a cosmological variable, say H, since, if Einstein's equations are invoked, the Λ term as such in these equations cannot exist in Minkowski spacetime unless the VED is exactly ρ Λ + ZPE = 0. So there is no cosmology whatsoever to do in flat spacetime, despite some stubborn attempts in the literature. Persisting in this attitude leads to the nonsense of having to cope with ∼ m 4 effects which must then be fine tuned among all the particles involved. This point has been driven home repeatedly e.g. in [40] and also recently in [41], see also [115]. A realistic approach to the VED within QFT in curved spacetime must get rid of Minkowski space pseudo-argumentations. The approach that we present here is fully formulated in curved spacetime and the vacuum energy density evolves with the participation of the curvature effects (powers of H) rather than with only powers of the masses, i.e. we pursue the successful renormalization program of [15,16]. Therefore, when the background curvature vanishes, we consistently predict that the non-trivial effects which are responsible for the value of the vacuum energy density and the cosmological constant disappear (and hence we are left with no Λ nor VED in the universe). Such is, of course, the situation in Minkowski space. In practice, however, we cannot reach that flat spacetime situation in our universe since there exists four-dimensional curvature at all times during the indefinite process of expansion. But by the same token such an impossibility evinces the fact that the VED and its dynamical nature is a direct consequence of the expansion process (and of the spacetime curvature inherent to it). The expected size of the VED and of Λ in our framework is indeed provided by the magnitude of the spacetime curvature, which is of the typical value of the measured Λ. It is therefore not caused by the quartic power of the masses of the fields (which is the very root of the CC problem in most approaches). These powers do not affect the running of the VED in our framework. To put it in a nutshell: the renormalized VED in our framework is like a small quantum 'ripple' imprinted on the existing (classical) background curvature owing to the vacuum fluctuations of the quantized matter fields. In the absence of the background curvature, the ripple would disappear too since it is proportional to it through the coefficient ν eff , which encodes the quantum effects from the quantized matter fields.
Following the same approach as for scalar fields, in the next section we compute the quantum effects contributing to the VED from the quantized spin-1/2 fields and express them in renormalized form using the same subtraction scheme devised in [15,16].

Renormalized ZPE for fermions
Thus, following the same prescription (19) as in the case of the scalar field, the renormalized form of the fermionic VEV of the EMT reads Since our aim is to study the ZPE we will focus into the 00th-component of the former equation for the moment. Alternatively, it is written as 7 T δψ 00 ren where we have used the calculational results for the unrenormalized components of the VEV of the EMT recorded in Appendix C and we have introduced the notation The last line of (80) contains all the non-divergent terms, which constitute a perfectly finite contribution and is made of finite parts from the 4th order expansion and of the entire 6th order term, which is fully finite but rather cumbersome. On the other hand, the first two lines in the last equality are a collection of terms that are individually divergent, but whose combination makes the integral convergent. In fact, by making use of simple algebraic manipulations at the level of the integrand one can show that explicitly. For instance, the rearrangement in the integrand shows that terms seemingly diverging as ∼ k 4 organize themselves to eventually converge as ∼ 1/k 2 . Needless to say, this is the consequence of the subtraction that has been operated. Similarly with the second integral in (80), whose individual terms are logarithmically divergent, but overall the integral is once more convergent thanks to the involved subtraction. The above renormalized result (80) would, of course, vanish for M = m ψ if we were to stop the calculation at 4th adiabatic order, so in case that one wishes to obtain the renormalized onshell result one has to either compute the exact unrenormalized EMT on-shell before subtracting the divergent adiabatic orders -which is possible but only in simpler cases such as in de Sitter space [91,92] -or one has to face the calculation of the adiabatic expansion up to 6th-order at least. In the last case the term T δψ 00 (6) (m ψ ) indicated at the end of Eq. (80) must be computed. This is what we have done here since an exact solution in the FLRW case is not possible. The necessary work to reach up to 6th adiabatic order for fermions is again significant, as it was previously for the scalar case [15,16]. The unrenormalized components of the EMT up to the desired order are explicitly collected in Appendix C. To subsequently obtain the renormalized EMT one has to implement the subtraction (79) and compute all the involved integrals. Despite the considerable amount of work involved, the final result to the desired order can nevertheless be presented through a rather compact formula, as follows 8 : The final equality corresponds to the expression in terms of the cosmic time (d()/dt ≡()) with H ≡ȧ/a. We point out that there is an explicit dependency on the Hubble function (and its derivatives) coming from G µν . This justifies the notation T δψ 00 ren (M, H), with two arguments, 8 We refer the reader to Appendix A.2 of [16] for the computation/regularization of the involved integrals (depending on whether they are convergent or divergent) with the help of the master DR formula quoted there. Use of DR can be convenient since in certain cases the needed rearrangement of terms in the integrand to verify that the overall integral is actually convergent can be complicated. Let us emphasize, however, that DR is only used as an auxiliary regularization tool for intermediate steps. The final result has no memory of this intermediate step, see e.g. Appendix B of [15] for an explicit nontrivial example. To be sure, no MS prescription is used for renormalization at any point of our calculation. The crucial difference between the ARP and the MS-like schemes is that the subtraction (79) involves not just the UV-divergences but also the finite parts.
where the dependence on the time derivatives of H is omitted for simplicity. We note that in the fermionic case there are no terms of O(H 4 ) in the evolution of the ZPE (and the VED, see next section), in stark contrast to the situation with scalars, see the last line of Eq. (20), where we can recognize terms of the form H 2Ḣ , HḦ andḢ 2 all of them of O(H 4 ). We also remark what has been previously anticipated: for M = m ψ (on-shell point) only the 6th-order terms remain, which are the ones in the last two lines of Eq. (82). These terms are relevant for the RVM mechanism of inflation in the very early universe (cf. Sec. 5.3). However, for the study of the renormalized theory at the point M (generally different from the on-shell mass point m ψ ) it is enough to consider the terms up to 4th adiabatic order, those in the first line of Eq. (82), see the next section.
So far, we have been able to provide the desired formula for the renormalized ZPE at the energy scale M up to 6th adiabatic order, as expressed by Eq. (82). This is, however, not the end of the story, since a proper expression for the VED needs to take into account also the renormalized parameter ρ Λ in the Einstein-Hilbert action (2), as this parameter is part of the unrenormalized vacuum action and after renormalization it also runs with the scale M , i.e. ρ Λ (M ). Both the renormalized ZPE and ρ Λ (M ) run with the scale and this will be crucial to study the properties of the renormalized VED. The running of the ZPE part between two different scales M and M 0 can be illustrated by considering the difference of the respective ZPE values at these scales. From (82) we find The finite parts, and in particular the 6th order terms cancel of course in the above difference, but the latter will be essential in the on-shell case since the result would be zero without these higher order effects 9 . We should notice that, in contradistinction to the case with scalar fields, there are no contributions of O(H 4 ) such as H 2Ḣ , HḦ orḢ 2 in the expression for the ZPE, as can be seen on comparing equations (20) and (82). For this reason it is unnecessary to use the higher derivative (HD) tensor (1) H µν (cf. Appendix A) as part of the renormalized Einstein's equations in the case of the fermion fields, again in contrast to the situation with the scalar fields -see [16] for details. Therefore, for fermions the subtraction at the two scales of the renormalized form of Einstein's equations can be done using the ordinary form of Einstein equations (3) without higher order curvature terms, and we find By comparison equations (83) and (84), and taking into account the tensorial structure of (84) and the explicit form of G µν in FLRW spacetime (cf. Appendix A) , we can perform the following identifications:

Renormalized VED
Once the renormalized ZPE has been obtained, the same consideration as for the scalar field case (see (22) and (23) ) can be repeated intact here, thus leading to the expression for the renormalized VED of the fermionic field: Now if the subtraction of scales is done, where in the last equality (85) was used. As expected, when written in terms of the ordinary Hubble function H in cosmic time, the evolution of the VED does not depend explicitly on the scale factor. For the sake of emphasizing the point, in the above equation we have explicitly indicated the cancellation of the terms carrying along the quartic powers of the masses, see the third equality in the above equation. As we can see, it is essential that the structure of the VED is obtained from the sum "VED = ρ Λ + ZPE", i.e. as in Eq. (23), since the referred cancellation occurs between the renormalized expressions of ρ Λ and ZPE upon being subtracted at the two arbitrary scales M and M 0 . This means that the two values of the VED at these scales are related in a very smooth manner: in fact, they differ only by a term proportional to H 2 , as it is obvious from (87). Even though Eq. (87) is formally correct, our job is not finished in the physical arena yet. Despite of the fact that such an equation describes the precise mathematical evolution of the VED with the renormalization scale, M , it is necessary to associate the latter with a suitable physical scale in order to extract useful phenomenological information out of it, exactly as in the companion studies of the VED for scalar fields previously presented in [15,16]. As pointed out in these references and also in Sec. 2.2 regarding the contribution from the scalar fields, the Hubble rate H is a characteristic energy scale (in natural units) of the expanding universe in the FLRW metric, and hence proves to be a natural candidate for a representative physical scale in this context. Whereby by following the same prescription used in the aforementioned references, we set the renormalization energy scale to M = H(t) (at the end of our calculations) in order to track the physical evolution of the VED. In other words, this prescription should allow us to explore the VED at different expansion history times H(t) in a physically meaningful way. In this way we obtain a well behaved evolution of the VED, which means that, given its value at one scale all other values at nearby scales are very close to it. The dynamics of the VED is slow and can be encoded into an effective contribution to the ν eff parameter, as we did for bosons in Eq. (26). The combined contribution from bosons and fermions to this parameter will be given in Sec. 5. Let us finally clarify the sense of this scale setting, Namely, the full effective action does not depend on M , of course, but the renormalized VED indeed does since the effective action of vacuum is only a part of the full effective action. The scale dependence on M from the other terms of the action, for example the terms carrying the running couplings of the RG-improved classical action, compensates for the M -dependence of the vacuum action. Put another way, only the full effective action (involving the classical part plus the nontrivial quantum vacuum effects) is scale-(i.e. RG-) independent. This is of course the standard lore of the renormalization group (RG), see also [41] for an expanded discussion. The choice of a particular scale helps of course in enhancing the physical significance of particular sectors of the full effective action. The procedure is of course akin to the usage of the RG in conventional gauge theories of strong and electroweak interactions, except that here one has to pick out an appropriate cosmological energy scale which is most adequate for the description of the universe's expansion. The distinguished scale H appears to be the natural choice if the universe where we live is indeed suitably described by the FLRW metric. In the next section we apply this approach to derive the important RG equation of the VED itself.

Renormalization group equation for the vacuum energy density
One can also compute the β function of the running vacuum generated by the fermionic quantum fluctuations. Only the adiabatic terms below 4th order carry M -dependence by definition since the higher orders are finite and hence are not subtracted in the renormalization procedure. As it was noted before, in contrast to the scalar case the terms of 4th adiabatic order are not present for fermions. The computation follows the same strategy as for scalars [16]. In this case we make use of equations (83) and (86), and we find The second equality holds immediately after computing the β-function of the parameter ρ Λ . From the first equation (85) we find that and hence contains a term proportional to the quartic power of the particle mass; what's more, there is an exact cancellation between the terms of the ZPE containing quartic powers of M and m ψ and the expression of β ρ Λ . The result (88) can also be consistently obtained directly from Eq. (87). Notice that neither the parameter ρ Λ nor the ZPE have physical meaning in an isolated way, only the sum makes physical sense and defines the VED in the present context. Let us compare the above results with those following from the contribution of one real scalar field φ [16]: and where we omit the O(H 4 ) terms in the scalar case (not present in the fermionic case) since it is enough to check the comparison at low energies. We can see that in both cases the β-function of the VED is proportional to β ρvac ∝ H 2 M 2 − m 2 , where m = m φ or m ψ , and therefore has a very smooth behavior thanks to the factor H 2 . In contrast, the β-function for the parameter ρ Λ in the gravitational action (which is often incorrectly identified as the VED in some explicit QFT calculations of the vacuum energy in the literature) behaves in both cases as β ρ Λ ∝ M 2 − m 2 2 and hence leads to undesired quartic contributions ∼ m 4 to the running. These are the problematic terms leading to fine tuning problems, but as can be seen these terms exactly cancel in β ρvac for the vacuum energy both for fermions and bosons in our renormalization scheme. Notice that there is a factor of 4 between equations (89) and (91) and have opposite sign. In a SUSY context, Eq. (91) should be multiplied by 4 to equalize bosonic and fermionic d.o.f in a given matter supermultiplet, all of whose members possess the same mass. Then β δφ , and the sum of the two coefficients will indeed vanish in a supersymmetric context: But this is, of course, not a cancellation of the β-function coefficients for the VED of bosons and fermions in the SUSY limit, but only the cancellation of the contributions to the β-function coefficient for the formal parameter ρ Λ in the EH action (2). This property is obviously connected with the discussion in Sec. 4.1 about the balance of UV-divergences between fermions and bosons. In a SUSY theory the quartic divergences cancel prior to any renormalization process, as we have noticed, and the resulting β-function for the parameter ρ Λ is zero. By the same token the running of the VED is freed from ∼ m 4 effects, which cancel among fermions and bosons in a SUSY context. The quartic powers are independent of the curvature of spacetime. However, the subleading divergences do depend on the background curvature and do not cancel at all, even in the exact SUSY limit 10 . The "residual" (finite) parts left in the renormalization process do not cancel either; they are actually proportional to the curvature of the FLRW background, R ∼ H 2 . This fact translates into a correction to the physical vacuum energy density of order ∼ m 2 H 2 both for bosons and fermions, which is far smaller than m 4 . So the finite, curvature dependent, terms that remain after ARP renormalization are de facto the most important ones for our purposes since they lead to the RVM form of the VED! The renormalization of the formal parameter ρ Λ , in contrast, has no physical imprint in the final result for the VED, except that the unwanted m 4 terms cancel against those involved in the ZPE, thus rendering the renormalized V ED = ρ Λ + ZPE free from quartic mass dependencies. From the above RG equations we may write down the total contribution to the β-function of the VED from the matter fields. Consider N f species of fermion fields with masses m ψ,ℓ for each species ℓ ∈ {1, 2, . . . , N f }, and similarly let N s be the number of scalar field species with masses m φ,j , j ∈ {1, 2, . . . , N s }. Some of these species may have the same mass, but this aspect is not relevant here, our formulas will include a summation over all contributions irrespective if some of them may be equal. The total β-function of the VED from an arbitrary number of quantized matter fields can now be cast as follows: 10 The SUSY considerations we have made here in passing only intend to clarify that in curved spacetime, irrespective of whether the quantized matter fields belong to a supersymmetric theory or not, the renormalization program is in any case mandatory to finally get rid of all the UV divergences. The calculations in this work, however, do not presume any SUSY context at all, not even a SUSY-broken theory. Our treatment of scalar and fermion fields is indeed completely general, in the sense that we are dealing with an arbitrary number of matter fields of both species without enforcing any balance between bosonic and fermionic d.o.f. -see Sec. 5 for more details.
The net outcome, therefore, is that the β-function of the vacuum energy density is free from undesirable contributions proportional to quartic mass powers of the quantized fields, ∼ m 4 , and hence these contributions do not appear in the renormalized theory. This is of course an extremely welcome feature of our renormalization framework, which is, on inspection of the above equation, fully shared by both scalar and fermion fields. Indeed, up to numerical factors fermions and scalar fields provide the same kind of leading contribution to the time evolution of the cosmological vacuum energy. Overall we find that the running of ρ vac depends only on quadratic terms in the fermion mass, namely ∼ m 2 ψ ℓ H 2 , which are of the same type as in the case of bosons, namely ∼ m 2 φ j H 2 , as discussed in Sec. 2.2 and previously demonstrated in great detail in [15,16]. These terms are actually very smooth owing to the presence of the H 2 factor. Integrating the RG equation corresponding to the β-function (93) one finds the expression for the evolution of the VED as a function of the renormalization scale M in the presence of any number of matter fields, see Sec. 5. In particular, integrating (88) for the case of one single fermion it is easy to verify that it leads exactly to (87).
The kind of much tempered behavior of the VED evolution that we have found here within our ARP renormalization program is of the sort that was expected on the basis of semi-qualitative RG arguments and constitutes the characteristic running law of the so-called Running Vacuum Models (RVM), see [40,41] and references therein. Thus, there is no need for fine-tuning in this scenario, since in such a renormalization procedure we have already gotten rid of the ugly contributions carried along by the quartic powers of the masses. In other words, the 'problem' with the quartic powers of the masses does not appear in the physically renormalized theory. While the running of the formal parameter ρ Λ with M indeed carries ∼ m 4 contributions, as it is obvious from the formulas above, this fact has no physical implication since ρ Λ is not itself a physical parameter (if taken in isolation) and the unwanted terms carried by it exactly cancel out in the β-function for the VED, as we have just proven. As a result, the running of the VED is much softer, the 'slope' is ∼ m 2 H 2 rather than ∼ m 4 . At variance with this result, in the context of the MS renormalization approach, in which ρ Λ runs with the unphysical mass unit µ coming from dimensional regularization, one is enforced to fine tune ρ Λ (µ) against the large contribution proportional to ∼ m 4 terms [41].

Renormalization of the fermionic vacuum pressure
Taking into account the perfect fluid form of the EMT associated with the vacuum, the corresponding pressure is defined through the iith-components. Any of them can be utilized owing to the assumed homogeneity and isotropy. So, it is just enough to compute the VEV of the 11thcomponent 11 : From (73) and using once more the expansion of the spin-1/2 fermion fields in Fourier modes (cf. Appendix B and Appendix C) the result can be expressed in the following way: 11 One can either compute the VEV of the T11 component, as we do here, or use the formula (24), which allows to compute the vacuum pressure from the 00th component and the trace of the EMT. The result is the same, of course, owing to the isotropy of vacuum. In Ref. [16], for instance, we presented the computation of the pressure for the scalar fields using this second method.
and where the explicit expressions (in WKB-expanded form) for the fermion modes h I k and h II k can be found in the appendices. Notice that there is a relation between ρ k and P k , which follows from (75) using the mode equations (37). This relation can be used as an alternative way to calculate T δψ 11 from T δψ 00 : For the sake of simplicity, the remaining discussions of this section will be restricted to the case of one single field. We shall retake the multifield case in Sec. 5. Following the same steps and considerations made in the previous sections for the 00th-component of the EMT, we reach the following expression for the renormalized value of the VEV of the 11th-component of the EMT up to 6th adiabatic order: We may now proceed to compute the vacuum EoS for the fermion fields up to the sixth adiabatic order. The best strategy is to compute first the pressure through Eq. (99), which can be inserted into the relation (94). Using next the VED expression (86) for fermions -with T δψ 00 given by (82) -the vacuum pressure can be seen to be equal to minus the VED plus some additional terms: The additional terms represent a small (but worth noticing) deviation from the classical vacuum EoS relation P vac = −ρ vac . The dominant vacuum EoS is still the classical one up to a leading correction of O(Ḣ) (the second term on the r.h.s of the above equation) and several sorts of higher order corrections of O(H 6 ) indicated in the last two lines. The ∼Ḣ correction in the first line of Eq. (100) can obviously be relevant for the present universe, and in particular it can modify the equation of state of the vacuum it to depart from −1 at present (cf. Sec. 5.4). The higher order terms in the last two lines, in contrast, might be relevant only for the very early universe, in principle. However, these terms involve time derivatives and hence vanish for H =const. This fact will have implications for our discussion of RVM-inflation in Sec. 5.3, since inflation can be shown to exist in this framework for H =const., cf. Sec. 5.3. So at the end of the day, the higher order terms in the last two lines of Eq. (100) become irrelevant both at low and high energies in this framework. The consequence is that the EoS of the quantum vacuum stays very close to −1 during inflation, in contrast to the vacuum EoS in subsequent eras of the cosmic evolution (cf. Sec. 5.4).

Trace Anomaly
It is a very well known result that if a field theory has a classical action which is conformally invariant, then the trace of the classical EMT vanishes exactly. For this it is necessary to work with a massless field, otherwise the presence of a mass breaks the symmetry since it introduces a fixed length scale. For instance, for a massless scalar field, This follows immediately from (6) and (7). However, it is also true that this simple result does not hold when one takes into account the quantum effects from the scalar field and constitutes the scalar part of the conformal anomaly, [1]. This follows after a careful study of the diverging part of the vacuum effective action, W Div eff , in which W eff was defined in Eq. (17). The part W Div eff is not conformally invariant for an arbitrary number of spacetime dimensions n (although W eff is so in the massless limit), except for the case n = 4. As a consequence, W Div eff receives a finite payoff for n → 4 owing to the existing pole 1/(n − 4) in it. Correspondingly, the VEV of the on-shell EMT receives a nontrivial contribution in the massless limit coming from the divergent part of the effective action, even in the case ξ = 1/6: The term δφ 2 contains some elements of 4th adiabatic order proportional to 1/m 2 φ , so that the corresponding limit results in a finite contribution. The same idea applies in the fermionic case, Here the term ψ ψ contains 4th adiabatic order terms that are proportional to 1/m ψ which make the limit non-trivial. Technically speaking (102) and (103) are not yet what we call the trace anomaly or conformal anomaly . This is due to the fact that the total effective action is conformally invariant and the corresponding EMT is traceless, so the part of the trace associated with the finite and divergent parts should be equal but with opposite sign in the conformal limit [1]. The anomaly is generated from the finite part, so its actual value for the scalar field case is where the conversion of the anomaly result into an invariant expression in the last step can be performed using the formulae of Appendix A. This result was explicitly verified in the calculation of [16]. We remark that for an arbitrary curved background the expression for the conformal anomaly is more involved [1]. However, since the spatially flat FLRW spacetime is conformally flat (i.e. conformal to Minkowski space) the contribution from the Weyl tensor vanishes identically and hence also its square (entering the anomaly). Additional terms beyond 4th adiabatic order decouple when m φ → ∞, satisfying the Appelquist-Carazzone decoupling theorem [116]. These terms are not finite in the massless limit, and hence do not take part of the anomaly. In practice we have derived the anomaly (104) from the unrenormalized trace of the vacuum EMT for scalar fields, T δφ , which is given in full detail in [16]. The corresponding conformal anomaly for fermions can be similarly extracted from the unrenormalized T δψ and it is a bit cumbersome as well, so we shall spare details here. We limit ourselves to provide the final result. Once more we can recognize the expression of the anomaly as a linear combination of finite terms of adiabatic order 4 which are independent of the mass scale. We find One natural question is related with the physical consequences of the conformal anomaly. It is well-known that it is a valuable theoretical concept encoding essential information on the VEV of the renormalized EMT [1], although it need not be itself part of the observable quantities of the renormalized theory. There are some attempts in the literature to remove the anomaly by particular prescriptions or definitions of the renormalized EMT [117]. This is also the case of the renormalization procedure employed in this work, as defined in (19) and (79), where the anomaly has no physical effects. The reason is that the on-mass-shell VEV of the EMT is subtracted at an arbitrary scale, M , up to 4th adiabatic order. Since the anomaly is of 4th adiabatic order and it is independent of the mass of the fields and, of course, also from the arbitrary renormalization point, it gets cancelled exactly in our ARP renormalization procedure. Alternatively, one can think in terms of the effective action. Indeed, in [16], we defined the renormalized effective lagrangian density off-shell at an arbitrary scale M , and it was shown by expanding it through the DeWitt-Schwinger series that it eventually leads exactly to the same renormalized EMT defined by (19). This result was obtained explicitly for a scalar field φ and can be repeated for fermions, although we shall not provide details here. Now the anomaly is related with the divergent part of the effective Lagrangian, corresponding to the lowest adiabatic orders (up to 4h order). As a consequence it gets once more exactly cancelled in (106) analogously to the subtraction of the EMT. As previously indicated, the anomaly part of the trace is contained in the unrenormalized trace of the EMT (even though the anomaly itself is a finite part of it). In our framework, however, the anomaly cancels since the anomaly is independent of the mass scale and our renormalized VEV of the EMT is defined through a subtraction of its value at two different scales, see equations (19) and (79). Thus the conformal anomaly is not involved in the renormalized expressions for the vacuum energy density and pressure in our framework. Despite it not having physical consequences in our approach, the explicit calculation of the anomaly is certainly useful as a nontrivial cross-check of our intermediate results.

Combined fermionic and bosonic contributions
Let us now determine the combined vacuum contributions from a multiplicity of non-interacting fermionic and bosonic degrees of freedom. As defined before (cf. Sec. 4.4 ), we consider N f species of fermion fields with masses m ψ,ℓ (ℓ ∈ {1, 2, . . . , N f }), and N s scalar field species with masses m φ,j (j ∈ {1, 2, . . . , N s }).

Running vacuum from an arbitrary number of quantized matter fields
The renormalized expression of the vacuum energy density is, in that case, Notice the appearance of the 00th component of (1) H µν , which is a HD tensor of O(H 4 ), hence of adiabatic order 4, see Appendix A. Its presence in the generalized Einstein's GR equations is indispensable for renormalization purposes and constitutes a UV completion of the field equations. No additional HD tensors are needed for conformally flat spacetimes [1]. In our case, (1) H µν is necessary for the renormalization of the short-distance effects produced by the quantum fluctuations of the scalar fields, as these involve O(H 4 ) corrections. However, as previously indicated in Sec. 4.2, the renormalized EMT for fermions does not contain O(H 4 ) terms. By using the expression of (1) H 00 in Appendix A we can recognize the tensorial structure of the various terms, and from it we can pin down immediately the running of the couplings: From the above formulas we can now use Eq. (107) to find out the difference between the values of the VED at two different scales: In the last line, the dots collectively represent all the terms of adiabatic 8 or beyond, which are not considered in our analysis. Notice that in the previous expression we have used the important relation (109), which is essential to cancel the quartic mass contributions from the matter fields. Following the same prescription that we used before to derive equations (25) and (26) for a single scalar field, we may implement now the scale settings M = H and M 0 = H 0 in order to compare the evolution of the VED between these two points, in the present case involving the full contributions from all the matter fields. For simplicity, let us call ρ vac (H) ≡ ρ vac (H, H) and ρ vac (H 0 ) ≡ ρ vac (H 0 , H 0 ) when using the above expression (112). The expansion history times H and H 0 can be arbitrary, of course, but for obvious reasons we choose H 0 = H(t 0 ) to be the value of the Hubble function at the present time, t 0 , and H = H(t) a value at a point in our past (t < t 0 ). Therefore, the running of the VED between these two points can be expressed as follows: Obviously, if the point H is in the nearby past we can neglect all the O(H 4 ) terms generated in the above expression since they are much smaller than the O(H 2 ) contributions. We will do this in the next section, where we study in more detail the low energy regime, in particular the late time universe where we live. Let us however clarify that the O(H 2 ) terms are dominant not only for the late time universe around our time, but in actual fact for the entire post-inflationary regime. Finally, we can extract the running of the gravitational constant from eq. (110), with the following result: (114) We discuss the running of ρ vac and G in terms of H in the next section.

The low energy regime: evolution of ρ vac and G in the present universe
Of paramount importance is the evolution of the VED and of the gravitational coupling G in the low energy regime, especially around our time. Therefore, following our prescription, we evaluate (113) for the late universe, when the dominant powers of H are the H 2 ones. Such an expression then boils down to where the function ν eff (H) is defined as follows: It is indeed a function evolving with the Hubble rate, but is almost constant since the dependence on H is very mild, as we shall make manifest in a moment. Let us emphasize that the O(H 4 ) terms correcting the r.h.s. of Eq. (115) are completely irrelevant for the current universe, and hence they can be safely ignored for the FLRW regime, that is to say, during the entire period following the inflationary stage (cf. next section). Therefore, equation (115) should actually be relevant for the full cosmological evolution that is accessible (directly or indirectly) to our physical measurements and observations. It is convenient to define the parameter This parameter is connected to the β-function (93) at low energies. Indeed, when we consider M = H in the low energy regime, it is obvious that H 2 ≪ m 2 for any particle mass, and hence Eq. (93) reduces to within a very good approximation. Quite obviously we can see that ǫ plays the role of coefficient of the low-energy β-function of the VED. However, the eventual coefficient that effectively controls the final evolution of the VED is actually enhanced with respect to ǫ by a big logarithmic factor. To see this, let us take the current limit (H → H 0 ) of the function (116): A simple rearrangement now shows that we can rephrase (116) in terms of ǫ and ν 0 eff : This formula is exact, but in practice some simplifications are perfectly possible. For example, consider the big logarithms ln m 2 i /H 2 0 (with m i any particle mass, boson or fermion) involved in ν eff but not in ǫ. For any known massive particle, we have ln m 2 i /H 2 0 ≫ 1, this being true even for the lightest neutrinos (recall that H 0 ∼ 10 −42 GeV= 10 −30 meV). Typically ln m 2 i /H 2 0 = O(100) in all cases. But as a matter of fact the only relevant contributions to ν eff (H) come from the heavy massive particles that belong to a GUT at a characteristic scale M X ∼ 10 16 GeV. For these particles (whether bosons or fermions, with masses m i ∼ M X ) we have m 2 i /m 2 Pl ∼ M 2 X /m 2 Pl and this number is not so small since it may thrust the value of ν eff up to ν eff ∼ 10 −3 , if one takes into account the large multiplicities of heavy fields existing in a typical GUT. This was first estimated long ago in [103]. Thus, it is natural to expect ν 0 eff ≫ |ǫ|. It follows that we can safely neglect the term proportional to ǫ in (120). It also means that we can neglect the very mild time-dependence of ν eff (H) and replace it with the constant coefficient ν 0 eff in which H is evaluated at the current time, H = H 0 . By the same token we can also ignore the numerical additive terms accompanying the big logarithms in (119). All in all, in very good approximation the evolution of the vacuum energy density can be described through the formula with an effective ν eff ≃ ν 0 eff given by As it turns out, in practice it all amounts to replace ν eff (H) → ν eff in (115) since the result still retains a great degree of accuracy. From the previous two equations, coefficient ν eff is seen to play the role of β-function for the running vacuum directly as a function of H. If we compare (122) with (117) we can see that ν eff and ǫ are roughly 'proportional' through a big log: Despite there being a summation over different masses, and hence such a proportionality not being strict, the above relation is nevertheless true in order of magnitude. The presence of the big log factor in ν eff makes the running of the VED faster than the tiny value of ǫ might convey at first sight. On the other hand, as we shall see below, it is ǫ alone that controls the (much softer) running of the gravitational coupling G, which does not receive any enhancement from big log factors. Equations (121) and (122) suffice to study the behavior of the VED near our time 12 . These simplified formulas have been previously used de facto to fit the value of ν eff from the latest cosmological data, see e.g. [62] and references therein. Here, however, we provide for the first time the full theoretical structure behind this parameter in the QFT context from the quantum effects induced by an arbitrary number of quantized matter fields. The typical fitting value obtained in the mentioned reference is ν eff ∼ 10 −3 and positive, which is well within the said expectations. This phenomenological determination picks up of course the net outcome from the various quantum matter fields involved in (122), which at this point cannot be discriminated in an individual way.
Finally, insofar as the running gravitational constant is concerned, it can be written using the same renormalization scale as follows: The former expression can be derived straightforwardly from eq. (114) by setting M = H and M 0 = H 0 , with H 0 being the current value of the Hubble function, and we have defined G N ≡ G(M 0 ) (the current value of the gravitational coupling). We follow exactly the same recipe as for the VED. In the low energy regime, where H 2 ≪ m 2 Pl , we can approximate with high accuracy the former expression by just where ǫ receives contributions from all the matter fields, see Eq. (117). Recall from (123) that |ǫ| ≪ |ν eff |, and also that the running of G(H) is logarithmic, in contrast to the running of ρ vac (H) which is quadratic in H 2 at low energies. Therefore, the running of G is much lesser than that of the VED. It may, however, be interesting to note that when H approaches m Pl the term ∼ H 2 /m 2 Pl in the denominator of the more accurate formula Eq. (124) could be dominant over the logarithmic one. If the multiplicity of matter fields is large enough, such a term could make the gravitational coupling to evolve asymptotically free at very large energies when we approach the Planck scale. Until that point, G increases at high energies for ǫ > 0. Beyond that point, decreases. As an additional cross-check, we can see that the running of the vacuum energy density and of the gravitational coupling are compatible through the Bianchi identity. This can be translated into a local energy exchange between the vacuum fluid and the background gravitational field due to the quantum fluctuations. A deeper insight on the local (covariant) energy conservation and the Bianchi identity can be found in the previous works [16,17], where the reader may find a detailed derivation of the logarithmic evolution law, that is to say, equation (125), in the simpler scenario of one real scalar field. In fact, one finds that the β-function (90) for the VED running is crucially involved also in the local conservation law of the VED, which can be written in two alternative ways [16]:ρ The first equality expresses the fact that the non-conservation of the VED is due to both the running of ρ vac with M (i.e. the fact that β ρvac = 0) and also to the cosmic time dependence of M (viz.Ṁ = 0), whereas the second equality is a direct reflex of the Bianchi identity in Einstein's equations with variable ρ vac and G, and hence provides a link between the time variation of the VED and that of the gravitational coupling G. The former equation does not depend on the number or nature of the fields involved, and holds as long as the matter components are covariantly conserved on their own. The non-conservation of ρ vac , however, preserves the Bianchi identity thanks to the corresponding running of the gravitational coupling. This does not preclude, however, that one can still formulate scenarios where matter can exchange energy with ρ vac , but we do not address this situation here. Taking the leading terms of β ρvac from the r.h.s. of (93) for the present universe, and setting M = H, we immediately obtain the following differential equation where ǫ is the full expression (117) involving the contributions from all the matter fields. Its solution is precisely Eq. (125), as can be readily checked. It is also interesting to note from the above formulas that this framework predicts a (mild) cosmic time variation of the "fundamental constants", such as the gravitational coupling G and ρ vac , as a function of H(t). Whence, it implies a small evolution of these 'constants' with the cosmological expansion. The possibility for such a variation is not new. It has long been discussed in the literature [118] and is still a hot matter of debate and of intensive test by different groups, see also [119] and the ample bibliography provided in them. Specific theoretical models accounting for such a possible variation are manifold, and in some cases they imply a time-dependence of the running couplings and masses in the particle and nuclear physics world, see e.g. [120,121]. While most of the proposals are based on strict particle physics scenarios, in particular on GUT's, testing the evolution of the VED in curved spacetime is a novel feature suggested in our framework, which was actually put forward on more phenomenological grounds sometime ago in [122][123][124]. The QFT calculations presented in the current work provide indeed a solid theoretical support to these same ideas but from first principles. Recall the golden rule in this arena: when one "fundamental constant" varies, then all of them vary! The formulas discussed above concern important epochs of the cosmological expansion such as the radiation-dominated epoch, matter-dominated epoch and the current epoch in which the vacuum energy resurfaces and became finally dominant over matter. During the entire FLRW regime the dominant power of the Hubble rate in the VED is H 2 orḢ (which are of the same adiabatic order). The terms with powers of H (or of equal adiabatic order) higher than H 2 (indicated by O(H 4 ) in (115)) acquire real relevance much early on in the expansion history since only during the most primitive stages of the universe we encounter a truly high energy scenario. As previously anticipated, an interesting feature regarding these higher powers is that they can provide us with a possible mechanism for inflation. Namely, if the early cosmic era possesses a short period where H remains approximately constant and very large (typically near a GUT scale), the universe may go through a phase of exponential expansion in which the VED starts from a huge value which then quickly decays into radiation and triggers the ordinary FLRW regime. This situation is possible also in the RVM framework, and is called 'RVM-inflation' [16,77]. See also [56] for a 'stringy' version. We reassess RVM-inflation in the next section in the extended QFT context of the present considerations, where we now have both scalar and fermion fields.

Inflation from running vacuum
It was noted in [16] that the quantum effects computed from the adiabatic expansion lead to higher powers of the Hubble rate and its derivatives, which are irrelevant for the current universe but capable to bring about inflation in the very early universe. They are characterized by a short period where H=const., provided this constant value is, of course, very large, namely around a characteristic GUT scale. The regime H=const. in our case is totally unrelated to the ground state of a scalar field potential and therefore this new mechanism does not require any ad hoc inflaton field. As said, it is called 'RVM-inflation'. Here we consider the contribution from the fermions fields and provide a formula for the dominant term of the energy density receiving contributions from an arbitrary number of non-minimally coupled scalar fields and also an arbitrary number of fermions fields. The payoff from the latter stems from setting H =const in Eq. (82), where we can see that all the time derivatives of the Hubble rate vanish except for a single term which is proportional to H 6 /m 2 ψ . The contribution from a non-minimally coupled scalar field was computed in [16] and here we just combine it with that of fermions assuming any number of both species. Overall, we find that the total VED involving the contributions from bosons and fermions at very high energies (hence relevant for triggering RVM-inflation in the very early universe) can be put in the following fashion: where The terms collected in the function F (Ḣ,Ḧ, ... H...) depend on different combinations of powers of H involving derivatives of H in all cases, and hence they all vanish for H =const. Thus F = 0 for H =const. Overall we see that the dominant contribution is of the form ρ inf vac ∝ H 6 with a complicated coefficient C inf which depends on the number of scalar and fermions fields, their masses, multiplicities and also on the non-minimal couplings of the different scalars. In the case of fermions this contribution is seen to be negative-definite, whereas in the case of the scalars it can be positive. Let us note that during the inflationary period the EoS of the quantum vacuum is essentially −1, with very tiny deviations caused by terms which depend on the various time derivatives of the Hubble rate. To the extent that the condition H=const. is fulfilled these deviations are extremely small, see Eq. (100). In the next section we shall see that, in contrast, the EoS of the quantum vacuum in the present time can deviate from −1 by a small amount which is not as negligible as in the very early universe and therefore could be detected and even mimic quintessence behavior.
The solution of the cosmological equations proceeds along the same lines as in [16], except that now the fermionic contribution is also included but it only modifies the specific coefficient of H 6 . Therefore, one finds again that a short period of inflation can be generated with H ≈ const. and subsequently the vacuum decays quickly into radiation [16]: in whichâ ≡ a/a * , with a * is the transition point from the regime of vacuum dominance into one of radiation dominance, which can be estimated to be around a * ∼ 10 −30 (see Eq. (137) below). Moreover, H I and ρ I are the value of H and ρ vac , respectively, at the beginning of inflation, with ρ I = C inf H 6 I . Applying the Friedmann equation, we find where G I ≡ G(H I ) is the gravitational coupling at H = H I , the latter being value of the Hubble parameter at the inflationary era. Needless to say, the difference between G I and the usual G N is not very important here since the running of G is logarithmic, and hence the effect is very small as compared to the fast evolution of the H 6 term, so in practice we can neglect the running of G for these considerations. To trigger inflation in an effective way, we must have a positive coefficient C inf > 0. In the light of Eq. (129) we can see that this is perfectly possible since the couplings ξ j and masses of the fields can take a variety of values that make this possible, as can be shown in a devoted study that will be presented elsewhere. The masses of the relevant fields involved must be very large, say around a typical GUT scale, M X ∼ 10 16 GeV. This may not be obvious at first sight. A naive interpretation of the higher order terms of the VED, which are related with the 6th adiabatic order of the ZPE (see equation (82)), may give the erroneous impression that the relevant masses are to be the lightest possible ones, but this is by no means true since in such a situation the adiabatic expansion would break down. On the other hand, the analysis through the Friedmann equations reveals the correct dependency of the VED and of the Hubble function on the masses during the inflationary regime. From equation (129) it is obvious that C −1/2 inf ∝ m φ,ψ , where the notation stands for a linear combination of the typical masses of the matter fields. The inflationary parameters (132)-(133), therefore, depend on a positive power of m φ,ψ , and as a result the process of RVM-inflation is actually dominated by the heaviest masses, in contrast to naive expectations; namely, masses m φ,ψ ∼ M X ∼ 10 16 GeV of order of a typical GUT, as mentioned above. It follows that the same heavy masses which may generate a mild (but non-negligible) quadratic running ∼ H 2 of the VED (with a coefficient ν eff ∼ 10 −3 ) in the late universe can also be responsible for driving fast inflation in the primeval stages of the cosmological evolution. To see this feature more explicitly, let us recall that the differential equation driving the Hubble function in the presence of a high power H 6 in the VED (128) reads [125][126][127] where ω m = 1/3 is the EoS of matter in the relativistic epoch, and H I is given in (132). We have neglected the influence of the term H 2 and also of the constant term in the very early universe. It is obvious from (134) (132) and (133), we find that the order of magnitude of the physical scales involved in RVM-inflation is the following: up to numerical coefficients and multiplicity factors, of course. Thus, if the masses of the relevant matter fields lie in the expected range for a GUT, the right order of magnitude for the relevant physical parameters at the inflationary epoch can be obtained. Keeping in mind that the mechanism of RVM-inflation can also be motivated in 'stringy' scenarios [56][57][58][59][60][61], it should be natural to expect RVM-infation in the range between the GUT scale and the Planck scale. This is exactly what the above estimates suggest in order of magnitude. One more observation is in order. One can easily check from (131) that forâ ≫ 1 (i.e. a ≫ a * ) we retrieve the standard decaying behavior of radiation, ρ r (a) ∼ a −4 . This condition enforces the following relation between ρ 0 r (the current value of radiation energy density) and ρ I (the energy density at the inflationary time): where ρ 0 c ∼ 10 −47 GeV 4 is the current critical density. In the meantime the vacuum energy becomes negligible and does not disturb primordial BBN, see [16,17] for more details. See also [77,[125][126][127] for interesting phenomenological applications prior to the QFT treatment of RVM-inflation, first presented in [16].
Finally, we should stress that RVM-inflation is genuinely different from e.g. Starobinsky's inflation [128], as explained in detail in [77]. While it may be natural to conceive that a consistent inflationary model of the very early Universe should be a good candidate for an effective theory of quantum gravity, at least at energies much less than the Planck scale, RVM-inflation reveals itself as one such possible candidate, all the more if we take into account that a low-energy 'stringy' version of RVM-inflation has been also identified and sharing most of the virtues of the current QFT formulation [56,57].

Equation of state of the quantum vacuum
The quantum effects of the fields have an imprint on the vacuum equation of state, which is not exactly the traditional one P vac = −ρ vac . From the expressions of the renormalized energy density and pressure of the vacuum that have been obtained in the previous sections and considering their generalization to an arbitrary number of fermions and scalars, we arrive at the following expression for the EoS of the quantum vacuum: Here O H 6 stands for the terms of adiabatic order 6 or higher, such as H 6 /m 2 ,Ḧ 2 /m 2 , . . . (m = m ψ j , m φ ℓ ). All these higher order terms can be neglected during the postinflationary era.
The above relation shows that the quantum vacuum EoS is of dynamical nature. As we can see there is a deviation from the rigid value −1, which is the traditional EoS ascribed to the cosmological constant in the ΛCDM framework. The correction terms due to both bosonic and fermionic fields are small in the present era, in comparison with the constant term −1, but need not be negligible since the particle masses involved can be from a typical GUT, and hence one can estimate that the effective parameter ν eff -the parameter defined in Eq. (122) -could reach up to 10 −3 [103]. Furthermore, if we focus only on the O(Ḣ) ∼ O(H 2 ) terms relevant for the current universe and the radiation epoch we may also neglect the higher order adiabatic terms of O(H 4 ) in the last lines of Eq. (138). Following the steps of [17], the EoS can finally be written in a rather compact form as a function of the cosmological redshift: where ν eff contains the combined effects from fermions and bosons, see Eq. (119), and we have defined the normalized Hubble rate with respect to the present time (H 0 ): Here Ω 0 vac = ρ 0 vac /ρ 0 c ≈ 0.7, Ω 0 m = ρ 0 m /ρ 0 c ≈ 0.3 and Ω 0 r = ρ 0 r /ρ 0 c ≈ 10 −4 are the current fractions of vacuum energy, dust-like matter and radiation, respectively. The EoS formula may be further simplified if we neglect the effect of the small coefficient ǫ in Eq. (120). This is justified since ν eff ≫ ǫ owing to the logarithmic extra terms ln m 2 /H 2 0 contained in ν 0 eff , which can typically be of O(100), see (123). Thus, to within a very good approximation, we can write Notice that the term proportional to ν eff in the denominator cannot be neglected at large z since it becomes dominant. In this case, the EoS takes on the form where ν eff has cancelled. For example, for z large enough but within the matter-dominated epoch the dominant term in equation (140) is the ∼ (1 + z) 3 one, and we can see from (142) that the vacuum EoS then mimics matter since w vac ≃ 0. Similarly, at much larger values of z already in the radiation-dominated epoch, where ∼ (1 + z) 4 is the dominant term, then w vac ≃ 1/3 and the vacuum imitates radiation. Such a 'chameleonic' behavior of the quantum vacuum was first noticed in [17], and in fact the formula for the vacuum EoS that we have found here is a generalization for an arbitrary number of fermion and boson fields of the expression previously found in [17]. Last but not least, the evolution of the vacuum EoS in the late universe is no less remarkable and striking. From (141) we find In this approximation we recover once more the form (28), but in this case with ν eff involving the contributions from all the quantized matter fields. Taking into account that the last fits of the RVM to the overall cosmological data favor a positive value of ν eff > 0 and of order 10 −3 [62] (see also the previous phenomenological studies on the RVM reported in recent years [63][64][65][66][67][68][69][70]), we learn that the deviation of the EoS from −1 in the present universe is not completely negligible. The vacuum appears disguised as quintessence since w vac −1.
To summarize, the running vacuum mimics the EoS of the dominant component at a given time of the cosmic evolution, and at present is dynamical. This was first confirmed for scalar fields in [17]. In the present work, we have found that the formal structure of the EoS does not change when we evaluate the fermionic contribution. Therefore, we find that w vac is −1 during inflation; stays close to 1/3 in the radiation dominated epoch, and then close to 0 in the matter dominated epoch. Furthermore, at present it mimics quintessence since ν eff is found to be positive in the phenomenological tests of the RVM and hence w vac is slightly above −1. Therefore, the quantum EoS deviates from the traditional value w vac = −1 of the classic vacuum. The remarkable consequence is that we can have at present an effective quintessence behavior of the quantum vacuum without need of invoking ad hoc scalar fields. Finally, the vacuum EoS in the remote future will be as that in the very early (inflationary) universe, i.e. w vac → −1, since the cosmic expansion asymptotes towards a new inflationary period.

Conclusions
In this paper, we have evaluated the contributions to the vacuum energy density (VED) from the quantized matter fields in a semiclassical gravity approach. By using a regularization technique of Quantum Field Theory (QFT) in curved spacetime called adiabatic regularization and making use of a specific (off-shell) subtraction prescription amply tested in previous works [15][16][17], we have been able to calculate the mode functions and the renormalized zero-point energy (ZPE) from spin-1/2 quantum fields in a FLRW background up to sixth adiabatic order. Together with the contribution from the ρ Λ term in the Einstein-Hilbert action, we have obtained the properly renormalized VED. Since the corresponding computation for scalar fields had already been accounted for in the aforementioned works, we have put forward here the combined contribution to the VED from an arbitrary number of quantized matter fields. We did not consider interactions among them, however, as the free field calculation in curved spacetime is already rather cumbersome in itself. One interesting difference between the expression of the ZPE of these two types of fields is that in the fermionic case the terms of fourth adiabatic order (viz. those involving four time derivatives) are not present. The final result is that the overall VED of the quantized matter fields upon adiabatic renormalization appears to be a soft dynamical quantity with the cosmological evolution. This is a most remarkable outcome of the present study. More specifically, the VED shows up in the form of an expansion in powers of the Hubble rate H and its time derivatives, all these powers being of even adiabatic order, a property which is fully consistent (and expected) from the general covariance of the theory. Such a series expansion appears to take the canonical form of the running vacuum model (RVM), see [40,41] and references therein. This means, in particular, that the leading quantum effects obtained for the late universe are of second adiabatic order, thus ∼ H 2 and ∼Ḣ. Obviously, this may have consequences for the present universe, and these consequences have been tested in previous phenomenological works. What is more, these quantum effects turn out to have a positive impact on a possible solution to the ΛCDM tensions, see e.g. [62][63][64][65][66][67][68][69][70] 13 In this comprehensive study, we have also discussed some of the theoretical difficulties in trying to renormalize the cosmological term, Λ, and its relation with the VED. To start with, it should be emphasized that these are two different concepts that can only be properly related in non-flat spacetime. If Λ is taken to be the physically measured value of the cosmological term at present, Λ phys , then its relation with the current VED is ρ 0 vac = Λ phys /(8πG N ). However, at a formal QFT level these quantities have to be derived from a gravitational action in curved spacetime and a lot more of care needs to be exercised. Leaving for the moment quantum gravity considerations for a better future (viz. for when, hopefully, the quantum treatment of the gravitational field becomes possible), the more pedestrian renormalization of ρ vac in QFT in curved spacetime proves to be already quite helpful at present [15,16]. It shows, for example -and the present works attest once more for this fact -that the renormalized VED in the FLRW background is a mild dynamical quantity evolving with the cosmic expansion, and hence that ρ 0 vac is just its value at present. There is in fact no such thing as a rigid, everlasting, cosmological constant in the context of QFT in the FLRW background. In general, ρ vac = ρ vac (H) is a function of the Hubble rate and its time derivatives. What we call the 'cosmological constant' Λ appears in our framework as the nearly sustained value of the renormalized quantity 8πG(H)ρ vac (H) around (any) given epoch H, where G(H) is the renormalized gravitational coupling, which is also running, although very mildly (logarithmically) with H, i.e. G = G(ln H). At present, G(H 0 ) = G N and ρ vac (H 0 ) = ρ 0 vac , and this defines Λ phys = 8πG N ρ 0 vac in a precise way in QFT in FLRW spacetime (within our renormalization framework).
The longstanding and widespread confusion in the literature about the notion of cosmological constant, Λ, and that of vacuum energy (density), ρ vac , has prevented from achieving a proper treatment of the renormalization of these quantities in cosmological spacetime. In particular, the attempts to relate these concepts in the context of flat spacetime calculations are meaningless and their repeated iteration has been highly counterproducing [41].
In the simplified scenario considered here, where only interactions with the gravitational background are allowed, the VED is the sum of two contributions, a parameter in the effective action, ρ Λ , and the ZPE of the quantized fields. After renormalization, the VED depends on a scale M , and the setting M = H at the end of the calculation allows us to compare the VED at different epochs of the cosmic history, in a manner similar to the standard association made of the renormalization point with a characteristic energy scale of a given process in ordinary gauge theories. Thus the difference between the VED values at any two points of the cosmological expansion, say H(t 1 ) and H(t 2 ), provides a smooth running of the VED. Remarkably, such an evolution turns out to be free from the undesirable ∼ m 4 contributions that emerge from the quantized matter fields in other frameworks. As a result there is no fine tuning involved in the evolution of the VED in the present calculation. The VED, in fact, adopts the standard form of the RVM, which in the late universe reads ρ vac (H 2 ) = ρ vac (H 1 ) + 3ν eff /(8π)m 2 Pl (H 2 2 − H 2 1 ), where H 1 and H 2 can be, for example, the current value, H 0 , and another value H near our past. Finally, ν eff is a small parameter related to the β-function of the renormalization group running of the VED, whose value has been explicitly computed in this work from the fluctuations of the quantized matter fields. Depending on the sign of ν eff , the VED can mimic quintessence or phantom-like behavior.
Much earlier in the cosmic history, the higher powers of H (larger than H 2 and of even adiabatic order to preserve covariance) took their turn and could be relevant for the inflationary regime, in the sense that they had the capacity to trigger inflation through a mechanism that has been called 'RVM-inflation' [15][16][17]. While the scalar field contribution to this inflationary mechanism had been computed in the previous references, in this work we have accounted for the spin-1/2 fermionic contribution and combined the two types of effects for an arbitrary matter content. In both cases (scalar and fermion fields) the sixth order adiabatic terms O(H 6 ) had to be computed. Finally, the renormalized vacuum fluid's pressure, P vac , has been determined using the same QFT techniques as for ρ vac . Equipped with these nontrivial results the equation of state (EoS) of the quantum vacuum can be computed from first principles. The entire contribution from quantized matter fields (bosons and fermions) can be encoded in the effective ν eff parameter. We find that the EoS function w vac = P vac /ρ vac deviates from the traditional result -1, a fact which is worth emphasizing. This is true in most of the cosmological history, especially after inflation (which is the only period in our past where the vacuum EoS stayed very close to −1). It is no less noteworthy, as previously mentioned, that in the late universe, and most particularly near our time, the vacuum EoS behaves as quintessence for ν eff > 0, the latter being the sign preferred by the existing phenomenological fits to the overall cosmological data -see [62], for example. For higher and higher redshifts during the FLRW regime, the vacuum EoS mimics the equation of state of the dominant matter component (relativistic or non-relativistic) at the corresponding epoch. Such a peculiar behavior of the running vacuum energy density was referred to as "chameleonic" in [17]. The tracking of the EoS of matter by the vacuum ceases to exist in the late universe, where the DE epoch breaks through and w vac behaves as effective quintessence, the reason being that the EoS is then in the process to asymptote towards −1 in the remote future, exactly as it was in the primeval inflationary time. In fact, the inflationary process in the late universe is eventually resumed, but very slowly.
Overall, by combining the results from an arbitrary number of quantized matter fields we find that the main cosmic running of ρ vac depends on the quadratic terms in the boson and fermion masses times the square of the Hubble function, i.e. ∼ m 2 ψ H 2 and ∼ m 2 φ H 2 . These effects are obviously much softer than the naively expected (hard) contributions ∼ m 4 ψ and ∼ m 4 φ . As remarked, the soft terms have been profusely tested in phenomenological works on the RVM existing in the literature, see e.g. [62][63][64][65][66][67][68] -and the upcoming [129]. The QFT effects that we have computed here and in the preceding studies [15][16][17] provide a solid theoretical underpinning of those phenomenological analyses. They even bring to light new relevant features, such as the dynamical character of the EoS of the quantum vacuum, which is unprecedented in the literature to the best of our knowledge. In particular, they suggest that if in future cosmological observations we can collect clear signs that the EoS of the dark energy deviates from −1, such a feature could be explained by the running vacuum. It therefore opens the possibility that such observations (if confirmed) may be accounted from fundamental properties of QFT attributable to the fluctuations of the quantized matter fields in curved spacetime rather than to the existence of ad hoc quintessence fields and the like. This could be an extremely interesting smoking gun of this approach. The EoS dynamics is prompted here by the virtual quantum effects produced by quantized fermion and boson fields, the same kind of effects which trigger a smooth evolution of the vacuum energy density in cosmological spacetime. As it turns out from the above considerations, in the RVM framework the need for fundamental quintessence fields and also for inflaton fields subsides dramatically, for they can both be replaced by the running effect of the quantum vacuum.
On pure cosmological/observational grounds, the physical outcome of the theoretical framework presented here can be summarized in the following way. The renormalization of the vacuum energy density in QFT in the FLRW background leads to the non-constancy of the 'cosmological constant', Λ, in Einstein's equations and predicts a slow time variation of the vacuum energy density and of its equation of state at the present time, which departs slightly from the traditional EoS value −1 for the vacuum. This conclusion emerges from explicit QFT calculations in our approach and may point to a possible explanation for a wide range of cosmological problems that have been dealt with phenomenologically in the past in terms of ad hoc quintessence or phantom fields. In our context, the currently measured cosmological 'constant' is neither mimicked nor supplanted by any ersatz entity from the already too crammed black box of the dark energy. We could simply phrase it in a nutshell as follows: the quantum vacuum shows up here as if it were a form of dynamical dark energy, but it is (quantum) vacuum after all. In fact, it is the same vacuum producing inflation in the early universe by means of higher (even) powers of the Hubble rate H n beyond n = 2, and still leaving a smoothly evolving cosmological term in the late universe through the lowest possible power compatible with general covariance, which is H 2 . Formally, all these powers of the Hubble rate emerge as quantum effects that are part of the effective action of vacuum. Today's physical cosmological term appears as a quantity directly connected with the (properly renormalized) vacuum energy density in QFT in curved spacetime: Λ phys = 8πGρ vac . As a running parameter, it is sensitive to the fluctuations of the quantized matter fields. Taking into account that in our framework the scale of renormalization is linked with the cosmological expansion, represented by the Hubble rate, it turns out that Λ phys , despite it appearing as an approximately rigid term during a typical cosmic span around any given epoch, it is actually a physical observable in evolution during the entire cosmic history. Its time variation, which is ultimately of quantum origin, is very small at present but it helps in relieving the current tensions of the ΛCDM and it might eventually provide an explanation for the cosmic acceleration observed in our Universe from first principles.
In this case, the non-vanishing Christoffel symbols in the conformal frame are On the other hand, the Ricci scalar is and the 00th components of the Ricci and Einstein tensors are The renormalization program requires taking into account higher derivative (HD) terms in Einstein's equations [1]. In the particular case of FLRW spacetime, it is enough to consider just one of the characteristic higher order tensors, (1) H µν , given by the metric functional derivative of R 2 in the effective action (we refer once more to Appendix A.1 of [16] for more details): (1) Its 00th component is (1) and its 11th component is (1) We will also need the invariants which hold good for flat three-dimensional FLRW spacetime. For gamma matrices (in flat spacetime), the standard Dirac basis is chosen for our calculations with spin-1/2 fermions: where σ k (k = 1, 2, 3) are the usual Pauli matrices. In terms of the above γ α , the curved spacetime γ-matrices read γ µ (x) = e µ α (x)γ α , where e µ α (x) is the vierbein (cf. Sec.3).

B Appendix: Adiabatic expansion of the spin-1/2 field modes
In the main text (cf. Sec. 3) we have presented an iterative procedure which allows us to determine the two types of field modes h I k and h II k which are necessary to construct the 2-component spinor fields. They are functions of both momentum (k) and conformal time (τ ) and have the following structure: where Here W k is a real function playing an analogous role to the effective frequency introduced (with the same notation) in the scalar field case, Eq. (12). The superscript (n = 1, 2, 3, ...) indicates that the corresponding quantity is of adiabatic order n. The modes (152) are constrained to satisfy the normalization condition Some of the notation is similar to that of [91][92][93], although we use conformal metric and different conventions, and moreover we deal with FLRW spacetime rather than with de Sitter (where the EMT takes a simpler form). In addition, as explained in the main text, we perform the renormalization at the arbitrary scale M (not at the on-shell point). This is important in order to test the scaling dependence of the renormalized VED, which is the main aim of our calculation.
In what follows we use the notation ω k ≡ √ k 2 + M 2 a 2 . Recall that unless it is explicitly noted the mass scale involved is the off-shell point M . When the subtraction (79) is implemented within our renormalization procedure, one just sets M = m ψ in the subtracted part. The mass m ψ can be conveniently expanded in even adiabatic orders as After completing ℓ ≥ 1 steps in the process described in Sec. 3, we end up with Eq. (51) in the main text, which depends on the following expression: where ω (j) k represents a function of adiabatic order j. Functions W k (τ, M ), F (τ, M ) and G(τ, M ) in the ansatz (152) can be estimated with the help of this product. However the following clarifications may be necessary to better understand this process, together with the explanations already given in the main text, see Sec. 3: • For ℓ = 1 we only have performed one iterative step, and at this point we need to deal with the square root of as defined in (40). The dots "..." in (159) account for terms of adiabatic order 4 and beyond.
The square root of the previous result yields with From the r.h.s of (160), the first two terms, ω k and ω (1) k , are used in the first order approximation of the modes (see (59) and previous equations in the main text). Now suppose that we proceed with a further step in the iterative process, ℓ = 2. We have to deal with the square root of The introduction of whose expression can be seen in (63), does not alter the 0th nor the 1st orders of (159) and (160) since ǫ 2 is a sum of terms of adiabatic order 2 and higher: Nonetheless, the 2nd, 3rd, . . . adiabatic orders of (164) do not coincide with the ones of (160). Similarly, by going to the next iterative step, ℓ = 3, implies working with the square root of the product with The introduction of ǫ 4 does not alter neither the 0th, 1st, 2nd nor the 3rd adiabatic orders of (162) or (164), since ǫ 4 is a sum of terms of adiabatic order 4 and higher. By the same token, the 4th and 5th adiabatic orders and beyond in (162) do not coincide with the ones in (165). Similar considerations apply to the square root of these quantities.
We can sum up this by saying that for each iterative step we can compute two consecutive adiabatic orders of (158) that will not be altered by the subsequent steps. Then, after ℓ steps, the 0, 1, . . . , 2ℓ − 1 adiabatic orders of the product (158) are trustworthy for the estimation of the modes.
• The r.h.s. of (158) has a pure imaginary part conformed by the odd orders, precisely those which do no take part in the computation of W k : However, because of the factor −i in front of the integral, the imaginary terms in the integrand become real and are then necessary for the computation of F and G in (152).
• We did not specify the limits of integration in (51) for the following reasons. On the one hand, even terms take part in the pure imaginary exponential of (152). Now because the imaginary exponential does not appear in the final result of the relevant quantities that we compute in the main text (since they cancel in the products of a function times its complex conjugate) one might wrongly be led to conclude that W k does not influence the calculation of the EMT. This would, however, be incorrect since the derivatives of the modes h I,II k are present in these calculations. On the other hand, after integrating the odd terms without specifying the limits in the integral, there exists some residual freedom in the form of a set of functions of the momentum only (i.e. not depending on time). These are called f k , . . . where the superscript indicates the adiabatic order. If our goal is to compute an adiabatic expansion of the modes up to 6th order we need 7 arbitrary constants for h I , namely f (0) k , . . . , f (6) k . Similarly for h II , which we denote g (0) k , . . . , g (6) k .
• From Eqs. (37) and (38), it is clear that With this considerations in mind let's us put forward the adiabatic expansion of W k explicitly. As said, W k can be specified though the even terms of (158). We are interested to compute at least up to 6th adiabatic order (that means, at least, ℓ = 4 steps). Therefore we find: with ω The odd terms in the expansion (158) yield a real exponential contribution in the integrals involved in (152) and hence do not contribute to the expansion of W k in (155), but are nevertheless necessary to compute the amplitude of the modes. Notice that after computing the integral, the adiabatic order decreases by one unit, so in order to estimate the amplitude up to 6th order is mandatory to compute up to ω k . The corresponding integrals for these terms are listed below: −iˆτ ω Use of Mathematica [99] has been helpful to work out the above integrals. The computational strategy consists in using a sufficiently general ansatz for the structure of the result that is compatible with the adiabaticity order of the calculation, and then solve for the coefficients (form factors) of the ansatz from pure algebraic manipulations. The procedure has been illustrated with a specific example in Sec. 3.
Finally, by applying the normalization condition (156) for the modes at each order, there exists a constraint to fix the residual freedom mentioned earlier, which is parametrized by the time independent factors f k and g k (only depending on the momentum k): k Imf k Imf (4) k + 2Imf (2) k Imf k Imf (5) k + 2Imf (2) k Imf Notice that, imposing the former conditions is equivalent to claim that where the departure from 1 are just terms of adiabatic order 7 or bigger. Similarly for the functions g k . It can be shown that the observables made out of quadratic terms in the modes h I k , h II k (such as e.g. the EMT), depend on the particular values of f k in the form (178). It follows that they are not actual degrees of freedom if they satisfy the conditions (177). It is then safe to set particular values for the functions as long as quantities are computed up to 6th adiabatic order. The simplest solution satisfying the normalization conditions (177) is f Equipped with these results we are able to calculate the different orders of F (τ, M ) up to 6th order. A comparison between the general equation (51) and the ansatz (152), the different orders of F are conformed by the rightful combinations of terms of the denominator Ω k Ω k,1 Ω k,2 Ω k,3 and the real factors of the exponential exp −i´τ Ω k Ω k,1 Ω k,2 Ω k,3 dτ . Now the different orders of F are For h II k , one can make use of the relation G (n) (M ) = F (n) (−M ).

C Appendix: Adiabatic expansion of T µν for spin-1/2 fields
The unrenormalized components of the VEV of the EMT, T δψ µν , for spin-1/2 fermions can be obtained through the adiabatic expansion, which we compute up to 6th order 14 The various terms in the expansion (185) read as follows: T δψ 00 (0) T δψ

00
(2) To obtain the component T δψ 11 a similar expansion as in (185) holds. However, it is possible to use the following relation with the previously computed T δψ 00 component: As a result we find: T δψ

11
(2) We refer the reader to the master formula in Appendix A.2 of [16] for the explicit computation of the integrals in the above results. We refrain from providing additional details here, as the remaining calculations can be handled straightforwardly with the help of the mentioned formula.