In Ref. [278], the results from Refs. [89,90,279] were used to analyze the early universe EoS *p*(). The extracted EoS depicted on the left panel of Figure 8 covers the broad interval in energy density 10−2 ≤ ≤ 1016 GeV fm−3 and corresponds to the evolution periods down from the GUT era, through the EW era and the QGP era, into the hadron era. The apparent smoothness of the function *p*() hides in the complicated temperature dependence (*T*) = *g*eff(*T*)0(*T*) depicted on right panel of Figure 3.

While the GUT EoS *p*GUT = (0.330 ± 0.024) valid for 108 ≤ 1016 GeV·fm−3 (the red triangles in Figure 8) appears only slightly below the ideal gas limit, the hadronic-era EoS *ph* = (0.003 ± 0.002)+(0.199 ± 0.002) characterizes the region 1 GeV·fm−3. Most interesting is the intermediate region containing both the QCD and the EW epochs, which can be described by a single function *p*SM = *a* + *b* + *cd* with *a* = 0.048 ± 0.016, *b* = 0.316 ± 0.031, *c* = −0.21 ± 0.014, *d* = −0.576 ± 0.034 [278]. It is worth noting that the critical energy density *c* defined implicitly by the equality of the pressures in hadronic and QGP phases *ph*(*c*) = *p*SM(*c*) reads: *c* (1.2 ± 0.2) GeV·fm−3. As can be seen from the parametrization of the fit,

$$p\_{\rm SM} = p\_1(\varepsilon) + p\_2(\varepsilon), \quad p\_1(\varepsilon) = b\varepsilon, \quad b > 0, \quad p\_2(\varepsilon) = a + c\varepsilon^d, \quad a > 0, \quad c < 0, \quad d < 0,\tag{62}$$

throughout the whole QCD and EW eras, there are two independent contributions *p*1() and *p*2() to the overall pressure of the universe. While *p*1 is always positive, the second pressure *p*2 is negative up to (7 − 13) GeV·fm−3. Although the corresponding value of the trace anomaly can not be directly deduced from the fitted EoS (62)

$$
\Theta = \frac{\epsilon - 3p}{T^4} = \frac{\epsilon(1 - 3b) - 3a - 3c\epsilon^d}{T^4},
\tag{63}
$$

it is positive for - (3 − 4) GeV·fm−3, and for very large energy densities, it falls as Θ ∼ *g*eff−3/2, cf. Equation (2).

**Figure 8.** (**Left**) The combined EoS *p*(), ≡ *ρ* of QCD and EW matter, using non-perturbative results [89] extended to include other DoFs such as *γ*, neutrinos, leptons, EW, and Higgs bosons as well as perturbative results [90,279]. Adapted from Ref. [278]. (**Right**) Bulk viscosity *ζ* at *μB* = 0 as a function of the energy density . The top symbols stand for the SM contributions, while the bottom ones stand for the QCD contributions, only. Adapted from Ref. [280].

Note that the sound velocities

$$c\_{s,1}^2(\varepsilon) = \frac{dp\_1}{d\varepsilon} = b > 0, \qquad c\_{s,2}^2(\varepsilon) = \frac{dp\_2}{d\varepsilon} = c d\varepsilon^{d-1} > 0,\tag{64}$$

are well defined, making it possible for each of two components to represent the EoS of some substance. Let us add that analytical expressions for the scale factor of the universe *a*(*t*) and Hubble parameter *H*(*t*) deduced from the EoS (62) were recently discussed in Ref. [281].

While the first component of pressure (8) corresponds roughly to the EoS of the massless gas of non-interacting particles, the second one, in agreement with the asymptotic freedom of QCD [33], dies out with increasing . Interestingly, *p*2(), up to the constant additive term *a*, coincides with the generalized Chaplygin EoS used in Ref. [282] to describe the evolution from a phase dominated by non-relativistic matter to a phase dominated by the Cosmological Constant (or DE)5. For *d* = −1 in Equation (62), we obtain the ordinary Chaplygin gas [283,284] with EoS *p* = *c*/.

Another possibility of how to read the term *p*2() in Equation (62) is, in analogy with Refs. [264,265] and Equation (54), to interpret it as the density-dependent bag function B() = −(*a* + *cd*), cf. Equation (62). The latter may account for the density-dependent character of the physical vacuum due to, e.g., the instanton liquid [210,285]. Instantons [272], classical solutions to the Euclidean equations of motion, are localized in all the four dimensions and correspond to tunneling events between degenerate classical vacua in Minkowski space. Since tunneling lowers the ground-state energy, the instantons provide a simple understanding of the negative non-perturbative vacuum energy density. Yet, the Euclideanbased instanton model is not the only solution representing the QCD vacuum. It remains, in fact, questionable to what extent it represents the reality due to non-analyticity of the gluonic field operators and the associated color confinement property. Below, in Section 4, we elaborate on a recently proposed alternative picture and new gluonic vacuum solutions, which are readily formulated in Minkowski and FLRW spacetimes.

#### *3.3. Hydrodynamical Description of Dissipative Effects and the Early Universe*

According to the currently accepted scenario, see, e.g., Refs. [84,286], the evolution of the early universe must include a number of dissipative processes in order to explain the current large value of the entropy per baryon. Some of them, such as the decoupling of neutrinos during the radiation era [120] or different cooling rates of the fluid components in the expanding universe [287], can result from the conventional physics; others, involving more exotic mechanisms, assume entropy production via string creation [288] or the GUT phase transitions [289]. The hydrodynamical description of dissipative effects is summarized in Appendix B.

Let us now follow the evolution of the early universe in terms of the EoS including bulk viscosity *ζ* defined in Equation (A11). In Ref. [280], data from non-perturbative [89] and perturbative [90,279] SM simulations were used to study behavior *ζ* over a wide range of temperatures *T*, entropy densities *s* and energy densities *ε*. It was found that *ζ*/(*Ts*) decreases exponentially with *T* increasing. The bulk viscosity dependence on the energy density at zero baryon chemical potential *μB* = 0 is displayed on the right panel of Figure 8. Looking first on the QCD contributions only, it is apparent that the non-monotonic dependence of *ζ*() can be divided into four regions. The first, spanning 100 GeV/fm3, corresponds to the hadron-QGP phase. The second region, up to ∼5 × 107 GeV/fm3, contains both the QCD and the EW phases of matter. The third one seems to form an asymmetric parabola with focus at the critical energy density of the universe, *ρc* 1012 GeV/fm3 [280]. The fourth region shows a rapid increase in emerging from a non-continuous point.

In the SM contributions, as shown in the upper curve of Figure 8, besides gluons and (2 + 1 + 1 + 1) quarks, the contributions of the gauge bosons: photons, *W*±, and *Z*0, the charged leptons: neutrino, electron, muon, and tau, and the Higgs bosons: scalar Higgs particle, were also taken into account. An overall conclusion drawn in Ref. [280] is that over the entire range of energy densities, the SM contributions are very significant. It is worth mentioning that the characteristic structures observed with the QCD contributions only are almost removed when adding also the EW contributions.

It is worth noting that the bulk viscosity could be important even outside the realm of the Hot Big Bang. In Ref. [290], the theory of the inflationary epoch, covering cold and warm inflation, as well as the models of late universe expansion in the presence of bulk viscosity, were analyzed. Assuming that the viscous effects during the inflationary epoch can be represented by a generalized and inhomogeneous EoS of the form:

$$p = -\epsilon + A\epsilon^{\beta} + \zeta(H) = -\epsilon + A\epsilon^{\beta} + \overline{\zeta} \left(\frac{8\pi G\epsilon}{3}\right)^{\gamma/2},\tag{65}$$

where *A*, *β*, ¯ *ζ*, and *γ* are positive constants and *ζ*(*H*) = ¯ *ζHγ* is the Hubble parameterdependent bulk viscosity, the authors have studied the behavior of various inflationary observables. Let us turn to important implications of non-perturbative QCD vacuum dynamics in cosmological evolution.

#### *3.4. Theory of Hot Meson Plasma Interacting with the QCD Vacuum*

Shortly after the confinement phase transition, the universe enters the state of hot meson plasma whose thermal evolution has been thoroughly explored in the framework of the Linear Sigma Model (L*σ*M) in Ref. [277]. This analysis exploits the non-perturbative method of a generating functional derived from the effective Lagrangian at finite temperatures (for further references on this method, see Refs. [291–294]) and accounting for the quartic self-interactions of the *σ*-meson only. The latter approximation corresponds to a realistic well-motivated configuration of the "hadron gas" interacting with the non-linear *σ*-field and reproduces some of the basic thermodynamical characteristics of the hot meson plasma observed also in other approaches, in particular, in L*σ*M-based scenarios [293,295] and Polyakov-loop extended Nambu–Jona-Lasinio (PNJL) models [296–298], but also features additional properties such as a possibility for *σ* → *ππ* decays in the plasma above the critical temperature.

The effective L*σ*M chiral Lagrangian accounting for the lightest scalar and pseudoscalar degrees of freedom *π*±, *π*0, *K*±, *K*0, *K*¯ 0, *η*, *η* , and the *σ*-meson, in the hadron gas approximation at *T* = 0, reads [277],

$$\begin{split} \mathcal{L}\_{\mathrm{eff}} &= \frac{1}{2} \partial\_{\mu} \sigma \partial^{\mu} \sigma + 2 \, \mathrm{g}^{2} \sigma\_{0}^{2} \sigma^{2} - \mathrm{g}^{4} \sigma^{4} + \\ &\frac{1}{2} (\partial\_{\mu} \pi\_{\mathrm{a}} \partial^{\mu} \pi\_{\mathrm{a}} + \partial\_{\mu} \eta \partial^{\mu} \eta + \partial\_{\mu} \eta' \partial^{\mu} \eta') + \partial\_{\mu} \mathcal{K} \partial^{\mu} K - \\ &\frac{1}{2} \Big[ 2 \kappa g^{2} (m\_{\mathrm{u}} + m\_{\mathrm{d}}) \sigma^{2} \pi\_{\mathrm{a}} \pi\_{\mathrm{a}} + \frac{2}{3} \kappa g^{2} (m\_{\mathrm{u}} + m\_{\mathrm{d}} + 4m\_{\mathrm{s}}) \sigma^{2} \eta^{2} + \\ &\frac{4}{3} \kappa g^{2} (m\_{\mathrm{u}} + m\_{\mathrm{d}} + m\_{\mathrm{s}} + \Delta\_{\mathrm{an}}) \sigma^{2} \eta'^{2} \Big] - \kappa g^{2} (m\_{\mathrm{u}} + m\_{\mathrm{d}} + 2m\_{\mathrm{s}}) \sigma^{2} \mathcal{K} \, , \end{split} \tag{66}$$

where *mu*,*d*,*s* are the constituent up, down and strange quark masses, respectively, *g* is the *σ* quartic coupling, and Λan 0.5 GeV is the gluon anomaly term that provides an explicit breaking of U(1)L × U(1)R symmetry. The QCD order parameter of the hadron mater *v*0 = 265 ± 15 MeV represents the amplitude of the quark-gluon (quantum-topological) condensate, such that

$$\epsilon\_{\rm top}(T=0) = -\left(\frac{b}{32} + \frac{(m\_{\rm ul} + m\_{\rm d} + m\_{\rm s})l\_{\rm g}}{4}\right) \langle 0| \frac{a\_{\rm s}}{\pi} G\_{\mu\nu}^{\rm a} G\_{\rm a}^{\mu\nu} |0\rangle \equiv -v\_0^4 \tag{67}$$

given in terms of the gluon correlation length *lg* (1.2 GeV)−1 [210,299], and the coefficient of the one-loop *β*-function of three-flavor QCD, *b* = 9. The quark-gluon condensate contributes together with the perturbative hadronic vacuum had vac emerging due to regularized contributions from meson fluctuations to the net QCD ground-state energy density,

$$\begin{split} \epsilon\_{\text{vac}}^{\text{QCD}} & \equiv \quad \epsilon (T=0) = \epsilon\_{\text{top}} + \epsilon\_{\text{vac}}^{\text{had}} = -v\_0^4 - \frac{1}{128\pi^2} (m\_{\sigma(\text{vac})}^4 + 3m\_{\pi(\text{vac})}^4) \\ &+ \quad m\_{\eta(\text{vac})}^4 + 4m\_{K(\text{vac})}^4 + m\_{\eta'(\text{vac})}^4) \simeq -7 \times 10^9 \,\text{MeV}^4. \end{split} \tag{68}$$

that satisfies the vacuum equation of state in the zero-temperature limit, (*T* = 0) = −*p*(*T* = 0). The hadronic vacuum term is negative had vac < 0 and appears to have a relatively small magnitude, i.e., had vac / QCD vac ≈ 0.15 [277]. The current *u*, *d*,*s* quark masses break the global chiral SU(3)L × SU(3)R symmetry explicitly, while it is also broken spontaneously by means of a *σ*-field expectation value,

$$
\sigma = \langle \sigma \rangle + \overline{\sigma} \, , \qquad \langle \sigma \rangle \equiv \frac{\overline{\sigma}}{\mathcal{g}} \, , \tag{69}
$$

In order to generalize the effective Lagrangian approach to finite temperatures, following theoretical foundations laid out in, e.g., Refs. [292,300–302], one should resume the daisy and superdaisy contributions. This is effectively achieved by utilizing an approximation that the expectation values of different fields are independent of each other, while omitting odd-point correlation functions and factorizing the four-point correlation functions into a product of two-point ones, for instance, *η*2*σ*˜ 2 = *η*2*σ*˜ 2, etc. Then, as a result of the minimization procedure of the non-equilibrium vacuum potential, one obtains the equation of state for the *σ*-condensate,

$$\begin{aligned} v^2 &=& v\_0^2 - 3g^2 \langle \hat{\sigma}^2 \rangle - \frac{1}{2} \kappa (m\_\iota + m\_d) \langle \pi\_a \pi\_a \rangle \\ &- \frac{1}{6} \kappa (m\_\iota + m\_d + 4m\_s) \langle \eta^2 \rangle - \frac{1}{3} \kappa (m\_\iota + m\_d + m\_s + \Lambda\_{\text{an}}) \langle \eta'^2 \rangle \\ &- \frac{1}{2} \kappa (m\_\iota + m\_d + 2m\_s) \langle \mathcal{R}K \rangle \, , \end{aligned} \tag{70}$$

as well as the equations of motion for the (pseudo)scalar field fluctuations about the evolving non-trivial ground state, for example, *∂μ∂μσ*˜ + *m*2 *σσ*˜ = 0, etc. Those fluctuations after quantization correspond to physical mesons with masses *m*2 *σ* = 8*g*2*v*2, etc. that are, in general, dependent on temperature. The vacuum values of the pseudo-scalar pseudo-Goldstone meson masses are found in terms of the light quark condensates via the Gell-Mann–Oakes–Renner relation [303–305], while the *σ*-meson has been identified phenomenologically with the scalar *f*0(500) state with mass, *mσ*(vac) 400 − 500 MeV.

The generating functional in the case of zeroth chemical potential is the free energy density of the considered meson plasma found in terms of the spatial part of the energymomentum tensor as follows,

$$\mathcal{F}(T, v, m\_{\sigma\prime}^2 \mathcal{M}^2) \equiv \frac{1}{3} \langle \mathcal{T}\_i^{\bar{i}} \rangle\_{\sigma}$$

whose minimization over the independent variables *m*2 *σ* and M2 ≡ *v*2 + *g*2*σ*˜ 2 provides the physical meson masses and the equation of state for the condensate *v*2. Substituting the meson masses expressed in terms of the temperature *T* and the order parameter *v* into F, one obtains the non-equilibrium Landau functional,

$$\mathcal{F}\_{\rm NE}(T, v) \equiv \mathcal{F}(T, v, m\_{\mathcal{F}}(T, v), \mathcal{M}(T, v))\,,\tag{71}$$

that provides the critical temperature of the chiral phase transition, where the finitecondensate phase becomes unstable, by means of

$$\frac{d^2 \mathcal{F}\_{\rm NE}}{dv^2}\Big|\_{T=T\_\varepsilon} = 0 \,, \tag{72}$$

while the stability condition of the low-symmetry phase reads *d*2FNE/*dv*2 ≥ 0. Finally, resolving the order parameter as a function of temperature, i.e., *v* = *v*(*T*), one finds the so-called equilibrium Landau functional as,

$$\mathcal{F}\_{\rm E}(T) \equiv \mathcal{F}\_{\rm NE}(T, v(T)) \; , \tag{73}$$

that matches the usual free energy definition. A variety of thermodynamic observables of the meson plasma are then derived in terms of FE(*T*) such as pressure, entropy density, energy density, heat capacity and the speed of sound squared,

$$p(T) = -\mathcal{F}\_{\rm E}(T), \quad \sigma(T) = -\frac{d}{dT}\mathcal{F}\_{\rm E}(T), \quad \varepsilon = \mathcal{F}\_{\rm E} + T\sigma, \quad \varepsilon\_V = T\frac{d\sigma}{dT} = \frac{d\varepsilon}{dT}, \quad u^2 = \frac{\sigma}{\varepsilon\_V}, \tag{74}$$

respectively. While at *T* = 0, the QCD vacuum equation of state, = −*p*, is satisfied, both pressure and energy density grow with *T* due to positive particle contributions. The total energy density of the "plasma + condensate" system vanishes at *T*=0 237 MeV for *mσ*(vac) 500 MeV.

Using the above formalism, in Ref. [277], the properties of different phases and transitions between them are explored in detail. For instance, at temperatures *T* > *Tc*, hadrons deconfine into quarks and gluons, while the condensates melts away, yielding a deconfined (or zero-condensate) phase of the QCD matter. Such a phase becomes metastable at temperature *T*0 < *Tc*, when the *σ*-meson fluctuation becomes massless *mσ*(*T*0) = 0, while *v*(*T*0) = 0 is still valid, and *T*0 is then found by resolving the extremum conditions on the generating functional. At another *T*1, two minima of *F*NE(*T*, *v*) become equal such that the zero-condensate phase stabilizes and the corresponding temperature is found from the following equation, *F*NE(*T*1, 0) = *F*E(*T*1). The first-order chiral phase transition to the zero-condensate phase then occurs at some temperature between *T*1 and *Tc*, and its strength increases with the *σ*-meson mass in the vacuum. Note, as expected from the first-order nature of this transition, the entropy density of the meson plasma grows with *T* and remains finite at all values of *T*, while the heat capacity appears to have a singularity, and the speed of sound squared vanishes at the critical temperature *T* = *Tc*.

The thermal evolution of the meson mass spectrum and the condensate of the L*σ*M is shown in Figure 9. Interestingly enough, both the condensate and the masses of all the mesons decrease with temperature for *T* < *Tc* = 438 MeV. The critical temperature becomes reduced if the fermions (quarks and baryons) are introduced [293]. As a result of the "hadron gas" approximation, the *σ*-mass rapidly falls to zero at *T*0 = 402 MeV and then grows much faster than the masses of other mesons (such that *mσ* > 2*mπ* almost for any values of *T*) in contrast with the corresponding predictions of other existing approaches such as PNJL. As a result, at low temperatures, the hadronic plasma is dominated by pions. There is also a significant phase co-existence domain of size (*Tc* − *T*0)/*Tc* 0.1. Pressure *p*(*T*), energy density (*T*) and the EoS ((*T*) − 3*p*(*T*) − *A*)/*T*4, where *A* = (*T* = 0) − 3*p*(*T* = 0) is the net vacuum contribution, are shown in Figure 10 from left to right, respectively. Due to the positive "hadron gas" contribution, both *p*(*T*) and (*T*) rise with temperature. Their profiles can be approximately reconstructed as a sum of the negativelydefinite constant QCD vacuum term (*T* = 0) and the contribution of the relativistic hadron plasma ∝ *T*4 (dashed lines), with the numerical coefficient *α* 3.5 and the effective number of DoFs, *gi* 9.

**Figure 9.** The condensate and meson masses as function of *T*. Here, *mσ*(vac) = 500 MeV, *g*2 = 0.4, *Tc* = 438 MeV, *T*0 = 402 MeV, *T*1 = 430 MeV.

**Figure 10.** (**Left**) Pressure *p*(*T*) as a function of temperature for the finite-condensate *v*(*T*) = 0 phase compared to that in the zero-condensate *v*(*T*) = 0 phase, *p*0(*T*) (solid lines), and to the approximated result (dashed line). (**Middle**) The same but for the energy density. (**Right**) The normalized EoS ((*T*) − 3*p*(*T*) − *A*)/*T*4 as a function of temperature, where *A* ≡ (*T* = 0) − 3*p*(*T* = 0) is the net vacuum contribution.

#### *3.5. Cosmological Constant and Vacuum Catastrophe*

The tight observational constraint on the DE EoS

$$w\_{\rm DE} = -1.03 \pm 0.03\,,\tag{75}$$

comes through a combination of the cosmological data from various sources such as the Type Ia supernovae, the baryon acoustic oscillations, the CMB anisotropies, and the weak gravitational lensing, etc.; for details on the confidence level and datasets, see Ref. [306]. These constraints are consistent with the standard cosmological model known as the ΛCDM (Λ term plus DM in the form of CDM, as two dominant components of the universe). Specifically, the DE is considered to be in the form of Cosmological Constant or Λ-term density,

$$
\epsilon\_{\rm DE} = \epsilon\_{\Lambda'} \qquad \epsilon\_{\Lambda} \equiv \frac{\Lambda}{\kappa} \; , \qquad \kappa = 8\pi G \; , \tag{76}
$$

expressed in terms of Λ-term in conventional normalization, Λ; see Equation (37). The latter satisfies the EoS *w*DE = −1 exactly. For a detailed recent review on achievements and challenges of the ΛCDM, see, e.g., Ref. [307].

Classically, an arbitrary Λ-term density 0 can be readily added to the right-hand side of the Einstein equations of GR that determine the macroscopic evolution of the universe,

$$R\_{\mu\nu} - \frac{1}{2}g\_{\mu\nu}R = \kappa(\epsilon\varrho g\_{\mu\nu} + T\_{\mu\nu})\,, \qquad T\_{\mu\nu} = -\frac{2}{\sqrt{-g}}\frac{\delta S\_m}{\delta g^{\mu\nu}}\,, \qquad S\_m = S\_m[\phi\_\nu \psi\_\nu A\_{\mu\nu} g\_{\mu\nu}]\,,\tag{77}$$

in terms of action of matter fields *Sm* and their energy-momentum tensor *Tμν*. In quantum theory, a non-trivial contribution to the ground-state energy density emerges as an average of the energy-momentum tensor over the Heisenberg vacuum state [308,309]

$$<\langle 0 \vert T\_{\mu\nu} \vert 0 \rangle = \mathfrak{e}\_{\text{vac}} \mathfrak{g}\_{\mu\nu} \,, \qquad \mathfrak{e}\_{\text{vac}} \not\equiv 0 \,. \tag{78}$$

The latter is proportional to the trace of the energy-momentum tensor; hence, it represents an effect of conformal symmetry breaking in a given fundamental QFT through either the formation of a Bose-Einstein condensate in a massless theory or through nonzero mass-dimensional terms in the original Lagrangian. It is generically ill-defined and should be renormalized, with the classical ("bare") 0 being treated as a counter-term in the initial Lagrangian, such that the divergences are cancelled between the two yielding a finite, but renormalization scale *μ* dependent, vacuum energy density, Λ(*μ*). This is the physical vacuum energy density that emerges in cosmological measurements performed at some fixed scale *μ* = *μIR* in the present universe, such that

$$
\epsilon\_{\Lambda}(\mu\_{IR}) \equiv \epsilon\_0 + \epsilon\_{\text{VAC}} \,. \tag{79}
$$

A macroscopic Cosmological Constant effect causing the universe to expand with acceleration (de-Sitter phase) is usually identified with energy density of the quantum vacuum that acquires contributions from all the incident vacuum subsystems existing in the SM and beyond. These would correspond to all quantum fields existing at energy scales ranging from the quantum gravity (Planck) scale, *M*PL ∼ 1019 GeV, down to the QCD confinement scale, *M*QCD ∼ 1 GeV—the maximal and minimal energy scales of particle physics, respectively. The current vacuum state of the universe is considered to be produced in the aftermath of the latest QCD phase transition associated with hadronization of the cosmological plasma.

In the framework of SM, besides the zero-point energy contributions to the ground state coming from each elementary particle, there are two major vacuum condensates whose characteristics are well established in particle physics—the weakly coupled classical Higgs condensate responsible for spontaneous EW symmetry breaking giving masses to the SM vector bosons and fermions and the strongly-coupled quantum-topological quark-gluon condensate in QCD.

One of the important aspects of the cosmological QCD transition epoch concerns the formation of the negatively-definite (CM) contribution to the ground-state of the universe that has received little attention in the literature so far. For illustration of this effect, let us consider the conformal anomaly term in the trace of the effective QCD energy-momentum tensor [310–312],

$$T^{\mu, \text{QCD}}\_{\mu} = \frac{\beta(g\_s^2)}{2} F^{a}\_{\mu\nu} F^{\mu\nu}\_{a} + \sum\_{q=u,d,s} m\_q \bar{q} q \,\,\,\,\,\tag{80}$$

where *mq* are the light (sea) quark masses *q* = *u*, *d*,*s*, *gs* and *β* are the QCD coupling constant and the *β*-function, respectively, and *Fa μν* is the gluon field stress tensor. Taking the vacuum average, we obtain

$$
\begin{split} \langle 0|T^{\mu,\text{QCD}}\_{\mu}|0\rangle &= \left. -\frac{9}{32}\langle 0| : \frac{\mathfrak{a}\_{\text{S}}}{\pi} F^{\mu}\_{\mu\nu}(\mathbf{x}) F^{\mu\nu}\_{\mu}(\mathbf{x}) : |0\rangle + \frac{1}{4} \left[ \langle 0| : m\_{\text{u}} \mathrm{d}u : |0\rangle + \langle 0| : m\_{\text{d}} \mathrm{d}d : |0\rangle \right] \\ &+ \quad \langle 0| : m\_{\text{s}} \mathrm{s}s : |0\rangle \right] \simeq -(5 \pm 1) \times 10^{-3} \,\text{GeV}^{4}, \qquad \mathfrak{a}\_{\text{S}} = \frac{g\_{\text{s}}^{2}}{4\pi}, \tag{81}
$$

representing the maximal value of the averaged quantum-topological QCD contribution to the physical vacuum energy density, whose spacetime dynamics is not fully understood and yet to be established. This contribution, also known as the quark-gluon condensate, is predicted by the theory of QCD instantons [210,299] and plays an important role in the chiral symmetry breaking and in dynamics of color confinement as well as in generation

of mass of light mesons in hadron physics as suggested by the Gell-Mann–Oakes–Renner relation [303–305].

The so-called "Vacuum Catastrophe" reflects the fundamental problem of consistent matching between the macroscopic observable Λ value, close to the critical energy density of the universe *ρc*,

$$
\varepsilon\_{\Lambda} \simeq 0.7 \rho\_{\varepsilon} \simeq 2.5 \times 10^{-47} \,\text{GeV}^4 > 0, \qquad \rho\_{\varepsilon} \equiv \frac{3H\_0^2}{\kappa} \,\text{}\,\text{}\tag{82}
$$

and the characteristic sizes of microscopic QFT predictions for the energy scale of each separate vacuum condensate [313,314]. Indeed, considering the topological QCD vacuum energy density QCD vac alone at macroscopic time and space separations typical for cosmological measurements, see Equation (81), its value is off by over forty orders of magnitude and has a wrong sign compared the observable Λ.

Indeed, such a large and negative contribution to the vacuum density creates a big problem for existence of the spatially flat universe, as the right-hand side of the corresponding Friedmann equation must be positive at all times (see, e.g., Refs. [67,307,315,316]). The presence of a negative cosmological constant of the QCD scale would necessarily trigger a fast collapse of the universe at the time scale of a microsecond, as the other components of the cosmological plasma energy-density decay as ∝ 1/*an* (with *n* = 4 for relativistic and *n* = 3 for non-relativistic media). This would prevent the universe from traversing the QCD horizon scale such that no macroscopic evolution would be possible. For a recent review on the status of this problem and the existing approaches, see, e.g., Ref. [67] and references therein.

A consistent resolution of the so-called "old" Cosmological Constant problem (why is Λ small and positive?) and the "new" Cosmological Constant problem (why is Λ non-zeroth and exists at all?) [307] may require a dynamical mechanism for compensation of different short-distance vacuum configurations in the infrared limit of the corresponding QFT. Such a vacuum self-alignment in the non-perturbative regime would be desired in order to avoid a major fine tuning of different parameters of the fundamental theory [317], and it may be considered as a new physical phenomenon [65,67,315,318].

As has been pointed out in Ref. [27], in the thermal SU(2)/SU(3) YM theories, the confining phase at low temperatures is expected to be void of energy density and pressure. Along these lines, one introduces a natural hypothesis about a heterogeneous structure of the non-perturbative QCD vacuum in the infrared limit of QCD [67]. Such a structure is characterized by the presence of at least two distinct vacuum subsystems which contribute with opposite signs to the net QCD vacuum energy density and mutually eliminate each other on average at large space and time separations, Δ*x* ∼ Δ*t*  1/ΛQCD. Then, it is reasonable to assume that a phase transition in the QCD vacuum has occurred in the course of cosmological evolution. Such a transition has led to a dramatic drop in the net vacuum energy density due to an (almost) exact cancellation between the different vacuum subsystems in the IR limit of QCD.

Generically, an ultimate goal would be to develop a universal framework that consistently describes quantum vacua (condensates) dynamics in real cosmological time at both macroscopic (IR limit) and microscopic (UV limit) separations and then identify the phase transitions between those. This remains one of the major unsolved problems of fundamental physics [319,320] (see also Refs. [67,314]).

Since the effect of the negative QCD condensate term must be eliminated somehow beyond the Fermi scale of QCD, Ref. [277] explores the simplest scenario for cosmological evolution by invoking an additional positively-definite cosmological constant in the framework of the effective meson plasma model, and also using the Bag-like model of the QCD crossover transition [321] for comparison. The basic working assumption adopted in this work was that a positive cosmological constant has been formed (stochastically) at the QCD scale together with the negative (topological) term, which means that the QCD vacua effects are dynamically eliminated at distances beyond the typical hadron scale. In Figure 11

(left), the net QCD vacuum energy density is shown as a function of the normalized scale factor in both scenarios with and without an additional positive Λ-term "compensator". In this scenario, as the QCD vacuum evolves with temperature and the universe eventually collapses, a "backward" QCD transition from the meson plasma to QGP may occur above the critical QCD temperature. Provided that there exists a mechanism for a bounce from the singularity (for possible scenarios for such a bounce, see, e.g., Refs. [322–325] and references therein), a possible series of such sequential "direct" and "backward" QCD transitions implies that the universe may, in principle, oscillate around the QCD epoch for some time. Eventually, the negative QCD vacuum effect is eliminated, the universe enters the phase of unbound expansion, and the standard cosmological evolution takes off [277] (see Figure 11, right).

**Figure 11.** (**Left**) The QCD vacuum energy density as a function of the (normalized) scale factor *a* in the meson plasma model (solid line) and the same quantity but with an extra positive Λ-term (dashed line) that exactly compensates the negative (topological) QCD term at large time-scales. (**Right**) A scenario of the universe oscillating during the QCD phase transition epoch with stochastic generation of a positively-definite QCD-scale Λ-term "compensator".

#### **4. Dynamics of Ground State in YM Theories**

Let us discuss the properties of quantized YM theories and their major implications in cosmology while focusing primarily on QCD-like strongly-coupled dynamics and its connections to confinement and to the nearly vanishing value of the cosmological constant.

#### *4.1. YM Ground State as a Time Crystal*

While Gross, Wilczek and Politzer proved the asymptotic freedom of non-Abelian gauge theories at large momentum transfers [31,32,138], Savvidy showed that the perturbative QCD vacuum at zero field strength is unstable [326] and thereby demonstrated the existence of the vacuum condensate (see also Batalin, Matinyan and Savvidy [327]). Nielsen and Olesen worked out an argument for why the explanation of vacuum condensates in terms of a homogeneous field filling the vacuum is problematic: such a field mode is unstable, at least on a static Minkowski background [328]. The ground state of QCD is a non-perturbative quantum-topological state of the YM theory that is of primary importance for the understanding of color confinement dynamics [329] as well as of hadronic and effective quark masses. For a thorough discussion on the QCD vacuum and its implications, see, e.g., Ref. [330] and references therein.

Despite the possible theoretical frameworks that are available for the description of the quantum ground state in YM theories, an interesting case of the formation of a metastable condensate in the Savvidy approach, with gravitational back-reaction providing a stationary stabilization, was studied in Ref. [74]. Intriguingly, the authors showed that the relaxation process induced by the QCD phase transition provides a novel mechanism for the production of GWs in the early universe. Such production is enabled through the SSB of time translation invariance that is reminiscent of what happens in the time-crystals that were theoretically predicted by Wilczek in Refs. [331,332] and observed by the Monroe group [333] (for a detailed review, see Ref. [334]). Within the setting of early cosmology, as discussed in Ref. [74] that included the perturbative all-order effective action, it was derived that the energy density of the quark-gluon mean-field decays monotonically in time, while the pressure density undergoes violent oscillations at the characteristic QCD scale. This mechanism entails the generation of a primordial multi-peaked GW signal that eventually is shifted into the radio frequencies' domain. If detected, such a signal would represent an unprecedented echo of the QCD phase transition and it is, in principle, observable through forthcoming measurements at the FAST and SKA telescopes.

The scenario depicted in Ref. [74] requires the emergence of a quark-gluon condensate characterized, as mentioned above, by violent oscillations of the pressure density. Such oscillations would be periodic in multiples of the inverse characteristic scale ΛQCD 0.1 GeV. A very similar result that is consistent with this analysis was derived in Ref. [103] in which the authors deployed holography with the aim of analyzing relativistic collisions in a one-parameter family of strongly coupled gauge theories that were undergoing thermal phase transitions. An oscillating behavior of the pressure density was discovered also in this latter work, and again, the period of the oscillations was found to be a multiple of Λ−1 QCD. It was concluded that out-of-equilibrium physics smoothes out the details of the transition.

Most parts of analyses developed that involve lattice QCD tacitly assume an analytical continuation of the results obtained on the Euclidean space to the Minkowski spacetime. The underlying argument, as has been commonly advocated, relies on the idea that locally, on any FLRW spacetime, the dynamics of QCD can be studied on a "frozen" Minkowskilike background and, thus, results may be analytically continued to the Euclidean space. Relying on this assumption, a notable theorem due to Maiani and Testa [335] showed that the scattering of asymptotic states can be counter-Wick rotated only in the infinite volume limit (with the scale being the threshold amplitude) and for time scales much smaller than the level spacing due to momentum discretization.

Specifically, the authors of [335] started out with the Osterwalder-Schrader theorem that ensures that the Euclidean correlation functions can be analytically continued back to the Minkowski spacetime. However, this theorem heavily relies on the so-called reflexion positivity condition [336], and this condition is not fulfilled in the aforementioned cases when the FLRW dynamics are studied at time scales that exceed the Hubble time by one order of magnitude. Indeed, for the cases discussed, for instance in Refs. [74,103], nonperturbative effects that are originating from violent oscillations of the pressure density definitely spoil the time reflexion positivity condition. This insight suggests that the Maiani-Testa theorem is inapplicable in the context of the current discussion.

In fact, recent analyses developed according to different frameworks surprisingly confirm a quite different picture than the one on which the lattice QCD analyses are based: for instance, the determined period of oscillations in [74,103] that exceeds the Hubble time by one order of magnitude (in Planck units), as was mentioned above. Thus, effects of the spacetime curvature cannot be neglected due to the influence of non-trivial nonperturbative effects that are recovered beyond the characteristic-time of QCD, Λ−1 QCD. This provides an argument that should influence the confidence of the application of lattice QCD methods in cosmology; specifically when such methods aim to determine the order of the phase transition.

Finally, we mention, as a further approach to the problem of determining the order of the QCD phase transition in the early universe, a method provided by conformal field theories prescribed by the foundation of the modern understanding of quantum field theory and particle physics (for a comprehensive overview of the basic concepts, see Ref. [337]). Having been hitherto deepened so as to unveil the universal properties of scale invariant critical points, these frameworks were applied to describe continuous phase transitions in fluids and magnets as well as in many other materials. Substantial efforts were devoted to the study of non-perturbative strongly coupled conformal field theories, especially concerning their symmetries and theoretical constraints. Such work has opened the pathway to the development of the so-called conformal bootstrap [338,339], which has proven to be a successful framework in two dimensions and which finally was extended within the last decade so as to account for higher dimensionality including the physically relevant cases of three and four dimensions (for a detailed pedagogical review on the conformal bootstrap approach in *d* dimensions, see, e.g., Refs. [340,341]). Notably, significant progress in analytical methods has been achieved, shedding light on the possibilities of how to cast the bootstrap equations and, in parallel, more powerful numerical techniques were developed in attempts to find their solution [342]. This has brought about ground-breaking results including the determination of critical exponents and correlation function coefficients in the Ising O(*N*) models in three dimensions [343].

#### *4.2. Effective Action Approach*

In this section, the effective action approach is outlined in order to provide a background to the discussion of the contribution to the vacuum energy density through the trace of the energy-momentum tensor (EMT). In 1977, Matinyan and Savvidy [344] and Savvidy [326] investigated the asymptotic behavior of the effective Lagrangian density in gauge theories building on earlier work by Heisenberg and Euler and by Schwinger (see references therein). The behavior was studied using RG methods in order to relate the effective picture of strong fields to the short-range properties of gauge theories. The quantum corrections to the classical action as found by Schwinger were discussed both in these two publications and further summarized in a recent review [345].

The investigation begins with an examination of corrections to the classical action as suggested by Schwinger, i.e., corrections that allow for an expansion of the effective YM action in the gauge fields *A*¯ *a μ* (also known as connections):

$$\begin{split} \Gamma[A] &= \int \mathrm{d}x \, \mathcal{L}\_{\mathrm{eff}} \\ &= \sum\_{\boldsymbol{n}} \int \mathrm{d}\mathbf{x}\_{1} \cdots \cdot \mathrm{d}\mathbf{x}\_{\boldsymbol{n}} \, \Gamma^{(n)}\_{\mu\_{1} \cdots \mu\_{n}} A^{a\_{1}}\_{\mu\_{1}}(\mathbf{x}\_{1}) \cdots A^{a\_{n}}\_{\mu\_{n}}(\mathbf{x}\_{n}) \\ &= \mathcal{S}\_{\mathrm{cl}} + \mathcal{W}^{(1)} + \mathcal{W}^{(2)} + \dots \, \end{split} \tag{83}$$

Here, the effective Lagrangian density Leff undergoes a perturbative expansion so that the *n*-loop corrections provide a deviation from the classical action *S*cl. The charged vector connection *A*¯ *a μ*(*x*) ≡ 0|*Aa μ*(*x*)|0 is the vacuum expectation value of the field operator and Γ(*n*) is the one-particle irreducible (1PI) vertex function. To each order, *W*(*n*) provides the *n*-loop correction to the classical action.

The effective Lagrangian at all-loop order in an SU(*N*) YM theory can be defined in terms of an order parameter J and a running coupling *g*¯(J ). It should be noted that the latter is different from the bare coupling *g*YM of the classical theory [66,316]. In a non-stationary cosmological background characterized by the FLRW metric, the conventional effective (quantum) YM Lagrangian can be written, through a rescaling of the fields according to Equation (17), as

$$\mathcal{L}\_{\rm eff} = \frac{\mathcal{J}}{4\overline{\mathcal{J}}^2}, \quad \overline{\mathcal{g}}^2 = \overline{\mathcal{g}}^2(\mathcal{J}), \quad \mathcal{J} = -\frac{\mathcal{F}^{\rm a}\_{\mu\nu}\mathcal{F}^{a\,\mu\nu}}{\sqrt{-\mathcal{g}}},\tag{84}$$

where A*a μ* are the rescaled SU(*N*) connections, F*a μν* are the rescaled field-strength components, and the equality entering the covariant field-strength F*a μν* in a curved background ∇*μ*A*a ν* − ∇*ν*A*a μ* = *∂μ*A*a ν* − *∂ν*A*a μ* has been accounted for. Furthermore, *g* ≡ det(*gμν*), where *gμν* = *a*(*η*)2 diag(1, −1, −1, −1) is the FLRW metric as a function of the conformal time *η*. Now, J simplifies to

$$\mathcal{J} = \frac{2}{\sqrt{-\mathcal{g}}} \sum\_{a} \left( \vec{E}\_{\vec{a}} \cdot \vec{E}\_{\vec{a}} - \vec{B}\_{\vec{a}} \cdot \vec{B}\_{\vec{a}} \right) \equiv \frac{2}{\sqrt{-\mathcal{g}}} \left( \vec{E}^2 - \vec{B}^2 \right), \tag{85}$$

which is an expression that emphasizes the dependence on the components of the CE and CM fields: *Ea*, *Ba*. The running of the coupling *g*¯ as a function of the order parameter J in Equation (84) fully determines the dynamics of the effective YM theory and is given by the solution of the RG evolution equation [66,345]

$$2\mathcal{J}\frac{\mathrm{d}\mathcal{g}^2}{\mathrm{d}\mathcal{J}} = \mathcal{g}^2\beta(\mathcal{g}^2) \,. \tag{86}$$

where *β* is the standard beta-function of the YM theory [346,347].

As an aside, it should be noted that apart from J , a second, independent, invariant may be constructed using the dual field strength [344]:

$$\mathcal{G} = -\frac{\mathcal{F}^{a}\_{\mu\nu}\mathcal{F}^{\*a}\mu^{\mu\nu}}{\sqrt{-g}} = \frac{4}{\sqrt{-g}}\mathcal{E}\cdot\mathcal{B}, \qquad \mathcal{F}^{\*a\,\mu\nu} = \frac{1}{2}\epsilon^{\mu\nu\rho\sigma}\mathcal{F}^{a}\_{\rho\sigma}.\tag{87}$$

This quantity is usually disregarded in the effective YM approach, since it vanishes for all fields that have orthogonal electric and magnetic components but it may as well, in principle, be incorporated into the effective Lagrangian.

In order to study the vacuum dynamics on cosmological scales, the spatially averaged quantity J should be considered, and two cases are distinguished in which: (i) J is positive, meaning that the averaged CE components *E*2 dominate over the averaged CM terms *B*2; (ii) vice versa, that is the case of a CM-dominated state with J < 0 that corresponds to a CM condensate. For the purpose of studying the basic features of the cosmological evolution of the CM and CE condensates in pure gluodynamics, it is sufficient to consider the effective SU(2) YM theory, since SU(2) subgroups can always be picked out of the SU(*N*) YM theory, and such a subgroup is the part that accounts for the cosmological application [66]. The explicit brackets ... will be dropped for the remainder of the discussion.

Applying the variational principle to the effective YM action as in the classical field theory, one straightforwardly obtains the effective YM equations of motion as described in Appendix C. Similarly, the EMT of the effective YM theory can be found as

$$T^{\nu}\_{\ \mu} = \frac{1}{\mathcal{S}^2} \left[ \frac{\beta(\mathcal{J}^2)}{2} - 1 \right] \left( \frac{\mathcal{F}^a\_{\mu\lambda}\mathcal{F}^{\mu\nu\lambda}}{\sqrt{-\mathcal{g}}} + \frac{1}{4}\delta^{\nu}\_{\ \mu}\mathcal{J} \right) - \delta^{\nu}\_{\ \mu}\frac{\beta(\mathcal{J}^2)}{8\mathcal{g}^2}\mathcal{J} \,. \tag{88}$$

which is particularly useful for our discussion of cosmological evolution of the quantum YM system in what follows.

#### *4.3. Mirror Symmetry of the Ground-State Solutions*

The one-loop ground-state solution behaves differently depending on the sign of J and of the running coupling *g*¯1. In the case of the CM solution, it is again noted that J < 0 and, hence, F > 0. In addition, one refers to a positive *g*¯2 > 0 in this case so that an absolute value in the one-loop effective Lagrangian, see Equation (A32), may be removed. Considering the CM branch to one-loop order, which corresponds to a choice of the initial condition in the RG equation such that *g*¯2 1(*μ*4 0) > 0, the RG solution as derived in Equation (A29) can be written on the compact form

$$\mathfrak{g}\_1^2(\mathcal{I}) = \frac{96\pi^2}{bN\ln(-\mathcal{I}/\lambda^4)}\,. \tag{89}$$

Here, note that *g*¯2 1(J ) > 0 when −J > *λ*4, with

$$
\lambda^4 \equiv \mu\_0^4 \exp\left[ -\frac{96\pi^2}{bN\mathfrak{g}\_1^2(\mu\_0^4)} \right],
\tag{90}
$$

and this gives yet another frequently used representation for the CM SU(*N*) Lagrangian [66]

$$\mathcal{L}\_{\rm eff, CM}^{(1)} = \frac{bN}{384\pi^2} \mathcal{J} \ln\left(\frac{-\mathcal{J}}{\lambda^2}\right). \tag{91}$$

The minimum of the effective Lagrangian at J∗ is taken to be the physical scale of the quantum YM theory, which in the case of a CM vacuum reads

$$-\mathcal{J}\_\* = \mu\_0^4,\tag{92}$$

which is the well-known phenomenon of dimensional transmutation. Readily from Equation (91), the minimal value of the effective CM Lagrangian reads

$$\mathcal{L}\_{\text{eff, CM}}^{(1)}(\mathcal{J}\_{\ast}) = \frac{\mathcal{J}\_{\ast}}{4\bar{\mathfrak{g}}\_{1}^{2}(\mathcal{J}\_{\ast})} < 0,\qquad \text{with } \bar{\mathfrak{g}}\_{1}^{2}(\mathcal{J}\_{\ast}) > 0,\tag{93}$$

and this is a negative value. In the standard notation, the CM minimum corresponds to

$$2g\_{\rm YM}^2 \mathcal{F}^\* = \varepsilon \mu^4 \equiv \Lambda\_{\rm QCD}^4 \quad \rightarrow \quad -\mathcal{J}\_\* \equiv 2\Lambda\_{\rm QCD}^4, \qquad \mathcal{L}\_{\rm eff, CM}^{(1)}(\mathcal{J}\_\*) = \frac{-\Lambda\_{\rm QCD}^4}{2g\_{\rm YM}^2}, \tag{94}$$

expressed in terms of the conventional scale, ΛQCD.

The exact ground-state solution for positive J∗ > 0 may be found immediately by inspection of the all-loop YM equation of motion; see Equation (A25). This is the CE condensate characterized by

$$\beta(\mathfrak{g}\_\*^2) = 2, \quad \mathfrak{g}\_\*^2 = \mathfrak{g}^2(\mathcal{J}\_\*), \quad \mathcal{J}\_\* > 0,\tag{95}$$

from which it follows that the equation of motion is trivially satisfied. It should be pointed out that such a special solution is universal for any SU(*N*) symmetry, i.e., it is independent of *N*.

The contribution to the vacuum by the CE condensate comes from the trace of the EMT (see e.g., Ref. [348]). The trace at the minimum as given by Equation (88) is

$$T^{\mu}\_{\ \mu} = -\frac{\beta(\tilde{\mathcal{G}}\_\*^2)}{2\tilde{\mathcal{G}}\_\*^2} \mathcal{J}\_\* = -\frac{1}{\tilde{\mathcal{G}}\_\*^2} \mathcal{J}\_\* \,. \tag{96}$$

The CE condensate hence contributes positively to the vacuum energy density, and this observation is intriguing when remembering that the CM condensate of Savvidy theory comes as a negative-definite contribution [345].

A striking and very interesting property of the YM effective Lagrangian will be discussed here, namely a mirror symmetry. It is apparent that the Lagrangian of Equation (84) is Z2-symmetric under simultaneous sign changes of J and *g*¯2. The invariance of the Lagrangian under this Z2 symmetry results in that the two condensates (CE and CM) are associated with two, apart from the overall sign, equal minima of the Lagrangian. Since the running of the coupling is a non-linear function of J in general, this symmetry may only be realized close to the ground state given in Equation (95) [66]. Therefore, in the vicinity of the ground state, the action is symmetric under the simultaneous transformation

$$\mathbb{Z}\_2 \colon \qquad \mathcal{J}\_\* \longleftrightarrow -\mathcal{J}\_\* \,, \qquad \mathbb{g}\_\*^2 \longleftrightarrow -\mathbb{g}\_\*^2 \,. \tag{97}$$

As implicated by the form of the *β*-function, see Equation (A28), the imposed symmetry forces an additional change of sign in precisely *β*: *β*(*g*¯2 ∗) ←→ −*β*(−*g*¯2 ∗). There are important consequences of this symmetry of the action. Mainly, the conventional CE condensate effectively becomes mapped onto the CM gluon condensate with J∗ < 0 and *g*¯2 ∗ > 0 [345].

Taking the order parameter at the ground state to be the physical scale of the YM theory, one has, for the two condensates *μ*4 0 = |J∗|. Then, for the CM branch, with *g*¯2 ∗ > 0, the running coupling for the one-loop solution may be rewritten as

$$g\_1^2(\mathcal{I}) = \frac{96\pi^2}{bN\ln(|\mathcal{I}|/\lambda^4)}, \qquad \lambda^4 = |\mathcal{I}\_\*|\exp\left[-\frac{96\pi^2}{bN\overline{g}\_1^2(\mathcal{I}\_\*)}\right],\tag{98}$$

c.f. Equation (89) and Equation (A29). Insertion of this expression back into the effective Lagrangian results in

$$\mathcal{L}\_{\text{eff, CM}}^{(1)} = \frac{bN}{384\pi^2} \mathcal{J} \ln\left(\frac{|\mathcal{J}|}{\lambda^2}\right),\tag{99}$$

c.f. Equation (91).

Due to the Z2 mirror symmetry, the minima of the effective Lagrangian on the two symmetry branches where *g*¯2 ∗ is either positive or negative come with the same value but with different scales: for the two condensates; the scale is modified by the exponential in Equation (98), so that

$$
\lambda\_{\pm} = |\mathcal{J}\_{\ast}| \exp\left[\mp \frac{96\pi^2}{bN|\mathcal{J}\_1^2(\mathcal{J}\_{\ast})|}\right].\tag{100}
$$

The upper sign stands for the CM condensate while the lower sign is the CE branch with *g*¯2 1(J∗) < 0.

As is clear from Equation (100) and from the discussion in Ref. [66], the minimum for which J∗ > 0 appears in the non-perturbative region defined by 0 < J∗ < *λ*4. Therefore, this minimum corresponds to the CE condensate found in Equation (95). The mirror minimum is then swiftly found by applying the Z2-symmetry transformation. Note that the physical scale for the mirror minimum associated with the CM condensate is exponentially suppressed relative to the minimum point |J∗| with the consequence that the CM condensate, with J∗ < 0, appears in the perturbative region where |J∗| > *λ*4.

Finally, we may return to the contribution of the CM condensate to the vacuum energy density. Applying the mirror symmetry to the CE minimum gives *β* → −2, and the equations of motion, as presented in Equation (A25), are no longer trivially satisfied. However, as pointed out in Ref. [66], the dynamical equation in the vicinity of the CM condensate becomes

$$\mathcal{D}\_{\nu}^{ab} \left[ \frac{\mathcal{F}^b \mu \nu}{\mathcal{g}^2 \sqrt{-g}} \right] = 0 \,. \tag{101}$$

This expression bears a close resemblance to the classical YM equations of motion that are valid in the vicinity of the ground state. The EMT for the CM minimum involves extra terms in comparison to Equation (96):

$$T^{\nu}\_{\ \mu} \left( \beta = -2 \right) = \frac{-2}{\bar{\mathcal{g}}\_{\ \ast}^{2}} \left( \frac{\mathcal{F}^{a}\_{\mu \lambda} \mathcal{F}^{a \nu \lambda}}{\sqrt{-\mathcal{g}}} + \frac{1}{4} \delta^{\nu}\_{\ \mu} \mathcal{J}\_{\ \ast} \right) + \delta^{\nu}\_{\ \mu} \frac{1}{4 \bar{\mathcal{g}}\_{\ \ast}^{2}} \mathcal{J}\_{\ \ast}. \tag{102}$$

However, in the trace of this expression, the terms within the parenthesis cancel out, yielding

$$T^{\mu}\_{\ \mu} = +\frac{1}{\mathcal{J}\_{\ast}^{2}} \mathcal{I}\_{\ast} \,. \tag{103}$$

The conclusion of Ref. [66] is therefore that the contribution of the two condensates to the vacuum energy density cancel each other with equal magnitude and opposite signs as long as the averaging over macroscopic volumes that contains many CE and CM vacuum "pockets" of typical microscopic length-scales of ∼ *λ*−1 ± ∼ Λ−1 QCD is considered:

$$\mathfrak{e}\_{\rm vac} = \frac{1}{4} \langle T^{\mu}\_{\ \mu} \rangle = \mp \mathcal{L}\_{\text{eff}}(\mathcal{J}\_{\ast}) \,. \tag{104}$$

Provided that the energy scales of "electric gluon" and "magnetic gluon" condensation are not the same, as was elaborated above, the condensates are formed at different spacetime separations. However, they evolve toward the same absolute value of the energy density, but with opposite signs, due to their cosmological attractor nature (see Section 4.4 below). This effectively causes the cancellation of "electric gluon" and "magnetic gluon" contributions to the QCD ground state at sufficiently large separations, i.e., in the deep infrared limit of the theory, although at the expense of a loss of homogeneity at typical length (Fermi) scales of QCD ∼ Λ−1 QCD. As the QCD vacuum appears to be locally inhomogeneous at these length-scales, gravity is expected to react to its electric and magnetic pockets in opposite ways such that the local metric fluctuations become averaged out beyond the length-scale of QCD to a small net effect compatible to that of the global observed cosmological Λ-term [65].

In the framework of the YM effective action approach, one may reach an intriguing conclusion that the exact compensation of the CE and CM gluon condensate components in the QCD vacuum (averaged over spacetime volumes above the QCD Fermi scale) is the necessary and sufficient condition for confinement in QCD. Indeed, as will be discussed below, the CE and CM vacuum pockets are always separated by non-analytic domain walls effectively blocking the gluon field from propagating over length-scales beyond the Fermi scale of each such pocket. The domain walls that separate different CM and CE pockets of the QCD vacuum prevent the color fields from propagating over macroscopic distances and, thus, effectively confine them within such pockets. An exact cancellation of their contributions to the net vacuum energy-density emerges in experimental observations as a complete disappearance of the gluon DoFs in the IR limit of the theory, i.e., beyond the QCD Fermi scale. In this picture of confinement, it would be natural to consider such pockets (with quarks and gluons being locked inside of them) as hadronic vacuum excitations. This is fully compatible with the classical limit of the YM theory where the conformal anomaly is absent and where only the radiation-like medium (hadron gas) remains at large separations. It remains to be seen exactly how the domain-wall picture of confinement readily formulated in Minkowski spacetime relates to a more standard centervortex mechanism of confinement [349–355] strongly supported by lattice simulations in Euclidean spacetime (for a recent review of the current status of this research field, see, e.g., Ref. [22,356]).

Note, with respect to the exact CM/CE cancellation and, hence, color confinement, the restoration of a discrete (mirror) symmetry between "electric gluon" and "magnetic gluon" contributions at the level of the ground state is an intrinsic property of the pure YM theory and the RG flow equations. This generic property is intricately connected to the fact that the QCD-induced component of the cosmological constant term vanishes for averages over macroscopic volumes of physical spacetime. A residual effect of the CE/CM cancellation, emerging due to an effective dynamical breaking of the mirror symmetry by gravitational interactions in the QCD vacuum, rigorously matches (in both the sign and an order of magnitude) the observed value of the cosmological constant [65].

Thus, the color confinement phenomenon and the tiny value of the cosmological constant are the direct and closely connected consequences of the mirror symmetry of the QCD vacuum in the infrared regime. An important implication of the domain walls in the QCD vacuum is that no analyticity of the scattering amplitudes can be assumed in such a case, causing potential problems with the standard imaginary time and Euclidean formulations such as lattice QCD that is relying on the analyticity properties and vacuum triviality of the theory.

#### *4.4. YM Cosmological Attractors*

The temporal evolution of the condensate(s) discussed in the subsections above may be described on cosmological scales where short-distance fluctuations are averaged (integrated) out. A simple background can be obtained by splitting the full gauge field into a background field component *A*¯*μ* and a fluctuating field *aμ*: A*μ* = *A*¯*μ* + *aμ*. This scheme

was developed in general by Wetterich in, e.g., Ref. [357,358], first for scalar fields in *SO*(*N*) and later for gauge fields, and it was further studied by, e.g., Gies [359] and Eichhorn et al. [360].

The (up to a rescaling) unique SU(2) pure YM theory will be parameterized in terms of a scalar time-dependent but spatially homogeneous field component (the background) on large scales due to the local isomorphism of the isotropic SU(2) gauge group and the SO(3) group of spatial 3-rotations [361–365]. One may therefore obtain a unique decomposition of the gauge field into this spatially homogeneous isotropic part (describing the YM condensate) and a non-isotropic/inhomogeneous component (accounting for the YM waves), the latter being the fluctuations, according to

$$\mathcal{A}^a{}\_k = \mathcal{U}(t)\delta^a{}\_k + \vec{A}^a{}\_k(t, \vec{x})\,. \tag{105}$$

Here, the decomposition has been performed using the gauge condition A*a* 0 = 0 [66]. It should be stressed that the fluctuations that parameterize the inhomogeneous YM wave modes interpreted as gluons in the field-theoretical framework average out over large distances by definition in the sense that *A*˜ *a k*(*t*,*x*) ≡ d3*x A*˜ *a k*(*t*,*x*) = 0 . Further, the homogeneous YM condensate itself can be considered positively definite, *U*(*t*) > 0, and it contributes to the ground state of the theory. The parameterization of the gauge field in SU(2) as a spatially homogeneous isotropic condensate and wave modes may be generalized to SU(3) for an application to QCD.

A quasi-classical theory of SU(2) YM quantum-wave excitations of the classical homogeneous condensate (i.e., without accounting for the vacuum polarization effects) has been thoroughly discussed in Ref. [276]. The formalism enables a proper extension to an arbitrary gauge and symmetry group with at least one SU(2) subgroup. Among the key results: an excitation of longitudinal wave (plasma) modes as well as an energy swap between the evolving homogeneous condensate and waves have been established in the linear and next-to-linear approximations. As is shown in Figure 12, the condensate tends to loose its energy, leading to the growth of YM wave amplitudes denoted as "particles". This represents a possible mechanism of particle production due to the dynamical vacuum decay which can be particularly relevant for cosmology and also in QGP production in heavy ion collisions. The effect has further been observed in the maximally supersymmetric N = 4 YM theory and in the more complicated two-condensate SU(4) gauge theory. As the next step, it would be important to perform an analogical study of quantum-wave dynamics in the effective action approach, i.e., in the case of a quantum YM vacuum, in order to study the impact of vacuum polarization phenomena on the energy balance in the "condensate + waves" system and hence on the growth of the wave modes.

Now, the dynamical behavior of the homogeneous YM condensate *U*(*t*) introduced in Equation (105) will be discussed. The Einstein equations for the pure YM theory in a non-trivial spacetime are obtained through the principle of variation starting from the effective action [315,316] and read as follows:

$$\begin{split} \frac{1}{\kappa} \left( \boldsymbol{\mathcal{R}}^{\boldsymbol{\nu}} - \frac{1}{2} \boldsymbol{\delta}^{\boldsymbol{\nu}} \boldsymbol{\mathcal{R}} \right) &= \\ \boldsymbol{\varepsilon} \boldsymbol{\delta}^{\boldsymbol{\nu}} \boldsymbol{\mu} + \frac{b}{32\pi^{2}} \frac{1}{\sqrt{-\mathsf{g}}} \Biggl[ \left( -\boldsymbol{\mathcal{F}}\_{\mu\lambda}^{\boldsymbol{\mu}} \boldsymbol{\mathcal{F}}^{\boldsymbol{\nu}\boldsymbol{\upsilon}\lambda} + \frac{1}{4} \boldsymbol{\delta}^{\boldsymbol{\nu}} \boldsymbol{\mathcal{F}}\_{\sigma\lambda}^{\boldsymbol{\mu}} \boldsymbol{\mathcal{F}}^{\boldsymbol{\mu}\boldsymbol{\varepsilon}\lambda} \right) \ln \frac{e|\boldsymbol{\mathcal{F}}\_{\mu\boldsymbol{\beta}}^{\boldsymbol{\mu}} \boldsymbol{\mathcal{F}}^{\boldsymbol{\nu}\boldsymbol{\alpha}\boldsymbol{\beta}}|}{\sqrt{-\mathsf{g}}\lambda^{4}} - \frac{1}{4} \boldsymbol{\delta}^{\boldsymbol{\nu}} \boldsymbol{\mathcal{F}}\_{\sigma\lambda}^{\boldsymbol{\mu}} \boldsymbol{\mathcal{F}}^{\boldsymbol{\nu}\boldsymbol{\alpha}\boldsymbol{\beta}} \Bigg], \\ \left( \frac{\boldsymbol{\delta}^{ab}}{\sqrt{-\mathsf{g}}} \boldsymbol{\partial}\_{\nu} \sqrt{-\mathsf{g}} - \boldsymbol{f}^{abc} \boldsymbol{\mathcal{A}}\_{\boldsymbol{\nu}}^{\boldsymbol{\epsilon}} \right) \left( \frac{\boldsymbol{\mathcal{F}}^{\boldsymbol{\mu}\boldsymbol{\alpha}}}{\sqrt{-\mathsf{g}}} \ln \frac{e|\boldsymbol{\mathcal{F}}\_{\mu\boldsymbol{\beta}}^{\boldsymbol{\mu}} \boldsymbol{\mathcal{F}}^{\boldsymbol{\alpha}\boldsymbol{\beta}}|}{\sqrt{-\mathsf{g}}\lambda^{4}} \right) = 0. \end{split} \tag{107}$$

**Figure 12.** The time dependence of the Hamiltonian corresponding to the YM condensate, *H*U, YM wave modes (or particles), *H*particles, as well as the interaction term between them, *H*int of the YM "condensate + waves" system (with total energy *H*) in the quasi-classical approximation of small wave amplitudes. Adapted from Ref. [276].

Here, *λ* = *ξ*ΛQCD, with ΛQCD ∼ 0.1 GeV being the QCD scale and where *ξ* has been introduced for scaling purposes, *e* is the base of the natural logarithm and *κ* is the gravitational constant. Finally, ¯ describes the ground-state energy density s.t. ¯ = QCD top + CC in terms of the quantum-topological contribution from QCD and the contribution from the Cosmological Constant. For confined QCD, QCD top ∼ −5 × 109 MeV4 for an SU(3) color symmetric theory, and this value may be extracted from the evaluation of non-perturbative quantumtopological fluctuations of the quark and gluon fields. The contribution from the Cosmological Constant as obtained from astrophysical measurements is CC ∼ 3 × 10−35 MeV4, which is a comparatively minuscule and positive value. The fact that QCD top contributes to the ground-state energy of the universe with a large negative value is a severe problem for all existing cosmological paradigms, or more accurately, for the particle physics theories from which it arises. This is because the large negative contribution must be compensated for, in any first-order approximation, to a remarkable precision, resulting in the observed cosmological constant value to an accuracy of a few tens of decimal digits.

The conformal dynamics of the ground state (the condensate) *U* = *U*(*η*) and of the scale factor *a* = *a*(*η*) are described by the following equations of motion as derived from Equations (106) and (107):

$$\frac{\hbar}{\kappa} \frac{a^{\prime\prime}}{a^3} = 4\vec{\varepsilon} + T^{\mu,II}\_{\mu\nu} \, , \qquad T^{\mu,II}\_{\mu} = \frac{3b}{16\pi^2 a^4} \left[ (\mathcal{U}^{\prime})^2 - \frac{1}{4} \mathcal{U}^4 \right] \, , \tag{108}$$

$$\frac{\partial}{\partial \eta} \left( \mathcal{U}^{\prime} \ln \frac{6e|(\mathcal{U}^{\prime})^2 - \frac{1}{4} \mathcal{U}^4|}{a^4 \lambda^4} \right) + \frac{1}{2} \mathcal{U}^3 \ln \frac{6e|(\mathcal{U}^{\prime})^2 - \frac{1}{4} \mathcal{U}^4|}{a^4 \lambda^4} = 0. \tag{109}$$

It should be noted that a particular exact solution to Equation (109) can be obtained if the logarithm evaluates to zero at all times: that is, if |*Q*| = 1 for

$$Q \equiv 6\epsilon \left[ (\mathcal{U}^{\prime})^2 - \frac{1}{4} \mathcal{U}^4 \right] a^{-4} \left( \mathcal{\zeta} \Lambda\_{\text{QCD}} \right)^{-4} . \tag{110}$$

This may be solved for the two special cases *Q* = ±1, and the solutions are shown in Figure 13. The homogeneous background *U* = *U*(*t*) displays quasi-periodic singularities in physical time in both cases. It should be stressed that the exact compensation of the CE and CM gluon condensate contributions to the QCD ground-state energy density, as discussed earlier, is realized in particular if the two components *Q* = ±1 co-exist in the universe. The cancellation happens over macroscopic distances as the average of the background

vanishes in the large-time limit, and importantly, this occurs without any fine tuning. Crucially, the cancellation will arise due to the time-attractor nature of the contributions coming from the two minima; a property that is demonstrated in the following.

**Figure 13.** The homogeneous QCD condensate amplitude oscillations. The homogeneous component *U*(*t*) displays quasi-periodic singularities in the physical time *t* = *a* d*η*, plotted here in units of the characteristic time scale Λ−1 QCD. To the left, the CE vacuum solution of Equation (110) with *Q* = 1 is shown, and to the right, the CM ditto with *Q* = −1 is displayed. *ξ* = 4.0 has been used along with initial conditions *U* = 0 and *U* = 0, respectively. These results are compatible with those of Ref. [66] up to the scaling of the figure on the right-hand side.

The conformal integral of Equation (108) is

$$\begin{split} \frac{3}{\kappa} \frac{(a')^2}{a^4} &= \bar{\epsilon} + T\_0^{0,II}, \\ T\_0^{0,II} &= \frac{3b}{64\pi^2 a^4} \left\{ \left[ (\mathcal{U}')^2 + \frac{1}{4} \mathcal{U}^4 \right] \ln \frac{6e|(\mathcal{U}')^2 - \frac{1}{4} \mathcal{U}^4|}{a^4 \lambda^4} + (\mathcal{U}')^2 - \frac{1}{4} \mathcal{U}^4 \right\}, \end{split} \tag{111}$$

which may be verified by differentiating this expression with respect to *η* and making use of Equation (109).

Equations (108) and (111) can now be combined in order to find a solution for the scale factor *a*, the trace of the EMT *Tμ μ* and the total energy density *T*0 0. This is possible since the latter equation incorporates the constraint of Equation (109). Solving this set of equations provides the benefit of obtaining solutions for observable quantities that must necessarily be smooth functions in time. Hence, the quasi-periodic singularities of *U*(*t*) may be avoided. Since *t* = d*η a*(*η*), Equations (108)–(111) may be recast in terms of the physical time as

$$\frac{\partial}{\partial \kappa} \left[ \frac{\ddot{a}}{a} + \frac{\dot{a}^2}{a^2} \right] = 4\ddot{\varepsilon} + T^{\mu,II}\_{\mu} \equiv T^{\mu}\_{\ \mu}(t) \,, \qquad \frac{3}{\kappa} \frac{\dot{a}^2}{a^2} = \ddot{\varepsilon} + T^{0,II}\_0 \equiv T^0\_{\ 0}(t) \,, \tag{112}$$

where the energy density of the gluon condensate and the trace in the one-loop effective YM theory read

$$T^{\mu,II}\_{\ \mu} = \frac{3b}{16\pi^2 a^4} \left[ a^2 \mathcal{U}^2 - \frac{1}{4} \mathcal{U}^4 \right],\tag{113}$$

$$T\_0^{0,ll} = \frac{3b}{64\pi^2 a^4} \left\{ \left[ a^2 \dot{\mathcal{U}}^2 + \frac{1}{4} \mathcal{U}^4 \right] \ln \frac{6\epsilon |a^2 \dot{\mathcal{U}}^2 - \frac{1}{4} \mathcal{U}^4|}{a^4 (\frac{\pi}{6} \Lambda\_{\text{QCD}})^4} + a^2 \dot{\mathcal{U}}^2 - \frac{1}{4} \mathcal{U}^4 \right\}. \tag{114}$$

In order to eliminate the explicit dependence on *U*(*t*) and, hence, the obstructing singularities in the two equations, [66] introduced a universal analytic function *g*(*t*) that parameterizes the relation between the trace of the EMT and the total energy density. The defining equation of this function is

$$T^{\mu,II}\_{\ \mu} - \mathbb{C} = (\mathcal{g}(t) + 1) \left[ T^{0,II}\_0 - \frac{\mathbb{C}}{4} \right],\tag{115}$$

$$\mathbb{C} \equiv -4\epsilon\_{\text{top}}^{\text{QCD}} = \frac{3b}{16\pi^2} \frac{(\not\!\!\!\!\/\Lambda\_{\text{QCD}})^4}{6\pi}.$$

Using *g*(*t*), Equation (112) can be written entirely in terms of continuous functions:

$$\frac{6}{\kappa} \left[ \frac{\ddot{a}}{a} + \frac{\dot{a}^2}{a^2} \right] = 4\epsilon\infty + \left( g(t) + 1 \right) \left[ T\_0^{0,II} - \frac{C}{4} \right],\tag{116}$$

$$\frac{3}{\kappa} \frac{\dot{a}^2}{a^2} = \epsilon\_{\text{CC}} - \frac{\mathcal{C}}{4} + T\_0^{0,II} \,. \tag{117}$$

Note here that *T*0,*U* 0 = *T*0,*U* 0 (*U*, *U*˙ , *a*). After excluding *T*0,*U* 0 above, the resulting equation for the scale factor that is left to be solved is

$$
\hbar \frac{\dot{\mathfrak{d}}}{a} + 3 \left( 1 - \mathfrak{g}(t) \right) \frac{\dot{\mathfrak{d}}^2}{a^2} + \kappa \epsilon\_{\rm CC} \left( \mathfrak{g}(t) - 3 \right) = 0. \tag{118}
$$

The general solution is6

$$\begin{split} a(t) &= \\ a^\* \exp\left[ \sqrt{\frac{\kappa \varepsilon\_{\rm CC}}{3}} \int\_{t0}^t \frac{1 + \sqrt{\frac{\varepsilon\_{\rm CC}}{\varepsilon\_0}} + \left( 1 - \sqrt{\frac{\varepsilon\_{\rm CC}}{\varepsilon\_0}} \right) \exp\left\{ \sqrt{\frac{\kappa \varepsilon\_{\rm CC}}{3}} \left( -3(t' - t\_0) + \int\_{t\_0}^t g(\tau) \, d\tau \right) \right\}}{1 + \sqrt{\frac{\varepsilon\_{\rm CC}}{\varepsilon\_0}} - \left( 1 - \sqrt{\frac{\varepsilon\_{\rm CC}}{\varepsilon\_0}} \right) \exp\left\{ \sqrt{\frac{\kappa \varepsilon\_{\rm CC}}{3}} \left( -3(t' - t\_0) + \int\_{t\_0}^{t'} g(\tau) \, d\tau \right) \right\}} \right] \end{split} \tag{119}$$

in terms of the initial values of the scale factor *a*∗ ≡ *a*(*t* = *t*0) and the total energy density 0 ≡ *T*0 0(*t* = *t*0), respectively. Note that CC 0.

The total energy density, *T*0 0(*t*), and the trace of the EMT, *Tμ μ*(*t*), both explicitly defined in Equation (112), can be found by insertion of the solution above into Equation (117) together with a manipulation of Equation (115). The result is

$$\frac{T\_0^0(t)}{^\text{€}\_\text{CC}} = \left[ \frac{1 + \sqrt{\frac{\epsilon\_\text{CC}}{c\_0}} + \left( 1 - \sqrt{\frac{\epsilon\_\text{CC}}{c\_0}} \right) \exp\left\{ \sqrt{\frac{\kappa c\_\text{CC}}{3}} \left( -3(t - t\_0) + \int\_{t\_0}^t \mathbf{g}(\tau) \, d\tau \right) \right\}}{1 + \sqrt{\frac{\epsilon\_\text{CC}}{c\_0}} - \left( 1 - \sqrt{\frac{\epsilon\_\text{CC}}{c\_0}} \right) \exp\left\{ \sqrt{\frac{\kappa c\_\text{CC}}{3}} \left( -3(t - t\_0) + \int\_{t\_0}^t \mathbf{g}(\tau) \, d\tau \right) \right\} \right]^2. \tag{120}$$

$$\frac{T^{\mu}\_{\ \mu}(t)}{^{\text{c}}} = 4 + \frac{4\left(g(t) + 1\right)\left(1 - \frac{\varepsilon\_{\text{C}}}{c\_{0}}\right) \exp\left\{\sqrt{\frac{\text{x}\varepsilon\_{\text{C}}}{3}} \left(-3(t - t\_{0}) + \int\_{t\_{0}}^{t} \mathbf{g}(\tau) \,\mathrm{d}\tau\right)\right\}}{\left[1 + \sqrt{\frac{\varepsilon\_{\text{C}}}{c\_{0}}} - \left(1 - \sqrt{\frac{\varepsilon\_{\text{C}}}{c\_{0}}}\right) \exp\left\{\sqrt{\frac{\text{x}\varepsilon\_{\text{C}}}{3}} \left(-3(t - t\_{0}) + \int\_{t\_{0}}^{t} \mathbf{g}(\tau) \,\mathrm{d}\tau\right)\right\}\right]^{2}}.\tag{121}$$

It shall be pointed out here that the above solutions for the scale factor, the energy density, and the trace of the EMT do not rely on any approximations but are the general solutions that can be obtained from Equation (112). These cosmological observables may therefore be studied on the full range from *t*0 to *t*, provided that *g*(*t*) is known.

For practical analyses, the auxiliary function *g* may be studied in the vicinity of the exact, large-time cancellation point where *Q*(*t*) ∼ 1. This is done by introducing an expansion of the YM energy density around the asymptotic value of the exact solution, where *T*0,*U*∗ 0 = *C*/4, such that

$$T\_0^{0,II}(t) \simeq \mathbb{C}/4 + \delta \epsilon(t), \qquad \delta \mathfrak{e} \ll \mathbb{C} \,. \tag{122}$$

Depending on the relation between the expansion parameter *δ* and the remaining scale CC, the time derivatives *a*˙ and *T*˙ 0,*U* 0 (*t*) take two different asymptotic forms. Firstly, in the case of large *δ*(*t*)  CC, these are

$$\begin{aligned} \dot{a} &\simeq \sqrt{\frac{\kappa}{3}} a \sqrt{\delta \varepsilon}, \\\dot{T}\_0^{0,II} &\simeq \sqrt{\frac{\kappa}{3}} (g(t) - 3) \left(\delta \varepsilon\right)^{3/2}, \end{aligned} \tag{123}$$

when keeping only the leading terms in *δ*(*t*) *C*. Secondly, in the opposite case when *δ*(*t*) CC, the same quantities instead become

$$\begin{split} \dot{a} &\simeq \sqrt{\frac{\kappa \varepsilon\_{\text{CC}}}{3}} a \left( 1 + \sum\_{n=1}^{\infty} \binom{\frac{1}{2}}{n} \left( \frac{\delta \varepsilon}{\epsilon\_{\text{CC}}} \right)^{n} \right), \\ \mathcal{I}\_{0}^{0,II} &\simeq \sqrt{\frac{\kappa \varepsilon\_{\text{CC}}}{3}} (g(t) - 3) \delta \varepsilon \left( 1 + \sum\_{n=1}^{\infty} \binom{\frac{1}{2}}{n} \left( \frac{\delta \varepsilon}{\epsilon\_{\text{CC}}} \right)^{n} \right). \end{split} \tag{124}$$

The difference between Equations (123) and (124) introduces only a very small correction to the period of *g*(*t*), and this correction can safely be neglected, as will be shown later in this section.

It is now possible to study *g* in the vicinity of the asymptote (Equation (122)) and in the two limits of *δ*(*t*) relative to CC. First, solve Equations (113) and (114) for *U*, *U*˙ in terms of *T*0,*U* 0 , *Tμ*,*U μ* . Then, insert the expansion of the energy density from Equation (122) in the resulting expressions as well as in Equation (115). Explicit expressions for *U*, *U*˙ in terms of the EMT components and the scale factor allows for the formulation of

$$
\partial\_t \mathcal{U}(t) - \mathcal{U} \equiv 0 \,. \tag{125}
$$

Explicit computation of the first term results in a differential equation for *g*(*t*) when the expansions in Equation (123) are inserted. The final form of such a differential equation is

$$\dot{\mathbf{g}}^4 - \frac{8\left(\frac{\mathbf{g}}{\hbar}\mathbf{A}\_{\text{QCD}}\right)^4}{3\varepsilon} \left(1 - \mathbf{g}^2\right)^3 = 0.\tag{126}$$

Its implicit analytic solution can be found for the inverse function *t*(*g*) over half a period of oscillation of *g*(*t*) as

$$t(\mathbf{g}) = -\frac{(6e)^{1/4}}{2\xi \Lambda\_{\text{QCD}}} \left[ \,\_2F\_1\left(\frac{1}{2}, \frac{3}{4}, \frac{3}{2}; \mathbf{g}^2\right) \mathbf{g} - k \right], \qquad 0 < t(\mathbf{g}) < T\_{\mathcal{S}}/2 \,. \tag{127}$$

The constant *k* ≈ 2.622 is defined through the above equation, and the initial condition *g*(*t*0) = 1 is adopted for simplicity. These conditions determine *g* = *g*(*t*) as a periodic quasi-harmonic function with unit amplitude. The period of oscillation may be found from Equation (126) as

$$T\_{\mathcal{S}} = \frac{2(6e)^{1/4}}{\,^{\mathcal{S}}\Lambda\_{\text{QCD}}} \int\_0^1 \frac{\text{dg}}{(1-\mathcal{g}^2)^{3/4}} = \frac{2k(6e)^{1/4}}{\,^{\mathcal{S}}\Lambda\_{\text{QCD}}}.\tag{128}$$

Using instead the expansions in Equation (124), the above calculation can be repeated for the case of *δ*(*t*) CC. The analogue of Equation (126), now with an additional term proportional to *α* 1, is

$$\begin{aligned} &\frac{\text{dg}}{\text{d}\tau} \pm 2(1-\text{g}^2)^{3/4} - a(1-\text{g}^2) = 0, \qquad \tau = t\frac{\text{g}^2\Lambda\_{\text{QCD}}}{(6e)^{3/4}},\\ &T\_{\text{g}}' = T\_{\text{\textdegree}} \left[1 + \frac{\text{\textdegree\textdegree}}{4k^2}a^2\right],\\ &a = \frac{2(6e)^{1/4}}{\text{\textdegree\textdegree}} \sqrt{\frac{\kappa\epsilon\_{\text{QCD}}}{\text{\textdegree}}} \approx 4.8 \times 10^{-24}, \qquad \text{for } \Lambda\_{\text{QCD}} \sim 210 \text{ MeV}. \end{aligned} \tag{129}$$

Hence, the additional term may safely be neglected.

Note that the function *g*(*t*) satisfies the following integral constraint

$$\int\_0^t \mathrm{d}\tau \, g(\tau) = \pm \frac{(6\varkappa)^{1/4}}{\int\_{\mathbb{S}} \Lambda\_{\mathrm{QCD}}} (1 - g^2)^{1/4}, \qquad \frac{nT\_{\mathbb{S}}}{2} < t < \frac{(n+1)T\_{\mathbb{S}}}{2},\tag{130}$$

where the upper sign corresponds to even *n* and the lower corresponds to odd *n*. It should be noted that this constraint can be used in Equations (119)–(121) in order to explicitly express the general solutions for *a*(*t*), *T*0 0(*t*) and *Tμ μ*(*t*) in terms of *g*. A very good analytic approximation to the exact *g*(*t*) may be constructed when keeping only the first two non-vanishing terms of the harmonic Fourier expansion s.t.

$$\log(t) \approx A \cos\left(\frac{2\pi t}{T\_{\mathcal{S}}}\right) + (1+A)\cos\left(\frac{6\pi t}{T\_{\mathcal{S}}}\right). \tag{131}$$

The amplitude *A* is found through

$$A = \frac{2}{k} \int\_0^1 \mathrm{d}y \, \frac{\mathcal{g}}{(1 - \mathcal{g}^2)^{3/4}} \cos \left(\frac{\pi}{2k} \int\_{\mathcal{S}}^1 \frac{\mathrm{d}x}{(1 - x^2)^{3/4}}\right) \approx 1.14\ldots \tag{132}$$

A comparison between this approximation and the exact solution *g*(*t*) is provided in Figure 14 from which it is clear that the approximation is indeed capturing the universal function to a very good accuracy. For a particular approximation to *g*(*t*) discussed above, the physical observables have been plotted in Figure 15 for illustration.

**Figure 14.** The time dependence of the quasi-harmonic universal function *g* = *g*(*t*). The exact solution in Equation (127) (solid line) has been extrapolated from the solution over a single period *Tg*/2. A harmonic approximation in Equation (131) (dashed line) captures the behavior of *g*(*t*) well.

**Figure 15.** Solutions for the total energy density *T*0 0(*t*) (**left**), the trace of the total QCD EMT *Tμ μ*(*t*) (**middle**) and the scale factor *a*(*t*) (**right**). The asymptotic values for which *Q* → 1 are indicated by horizontal lines in the left and middle panels, respectively. The initial conditions have been chosen as *U*0 = 0, *U*˙ 0 = (*ξ*ΛQCD)2/ √3*e*, *Q*0 > 1, *ξ* = 4.0 ΛQCD = 332 MeV and *κ* = 10−7 MeV−2. The energy density and the trace are plotted in dimensionless units, rescaled by Λ4 QCD and for illustrative purposes, CC was set to ∼ 0.5 % of ¯. These results are compatible with the qualitative picture in Ref. [66].

#### *4.5. SU(N) and the Functional RG Approach*

So far, we have addressed gluodynamics resorting to an all-order effective perturbative approach. Nonetheless, we can extend our investigation so as to include nonperturbative all-order analyses, resorting to the Functional Renormalization Group (FRG) approach [357,366–375]. This latter is a Wilsonian momentum-shell-wise integration method for the path integral, which was developed to delve into the dynamics of interacting quantum field theories and statistical systems in a non-perturbative way when couplings cannot be dealt with using perturbative techniques. A regulator function *Rk* is taken into account so as to suppress quantum fluctuations at momenta lower than some physical scale, i.e., at an IR cut-off scale *k*. This scale is in principle different to the one defined in the subsections above, namely *μ*0. The former denotes the scale above which all quantum fluctuations are integrated out, while the latter captures instead the one-loop renormalization scale and may further be extended to all-loop accuracy. The regulator function has been discussed by, e.g., Gies in Ref. [359]. A scale-dependent effective action that flows with *k* is then recovered, i.e., Γ*k*, which encodes quantum fluctuations effects at momenta larger than the IR cut-off *k*. Varying *k* then allows to smoothly interpolate among the microscopic/short-scale action and the full quantum effective action Γ*k*→0.

This elucidates why this procedure looks to be tailored ad hoc for cosmological applications, at the IR scale. The Wetterich equation for a non-zero background that fixes the running dependence on the cut-off scale in the FRG approach reads [357,366]

$$
\partial\_t \Gamma\_k = \frac{1}{2} \text{STr} \left( \Gamma\_k^{(2)} + R\_k \right)^{-1} \partial\_t R\_k \, \, \, \tag{133}
$$

where Γ(2) *k* , a matrix in the field space, denotes the second functional variation of the effective (running) action with respect to the field content of the theory7. Above, STr is the supertrace, including a summation over all field components and discrete indices, as well as all the eigenvalues of the Laplacian in the kinetic term. The FRG equation retains a dependence on the full (field-dependent) non-perturbative regularized propagator [367], namely, (Γ(2) *k* + *Rk*)−1.

The FRG approach has been then adapted to YM SU(*N*) theories [358–360]. Specifically, in Ref. [360], a numerical extrapolation among low and high energy scales for the full propagators was deployed in order to derive the gluon condensate. Refining these results in Ref. [376], the FRG approach was extended, within a cosmological setting, to the SU(2) case.

In Ref. [376], resorting to approximations that are necessary to solve the FRG equation, which is otherwise too complicated to provide analytic results, the authors considered

replacing Γ*k* in Equation (133) with the bare action *S*, allowing for integration on both sides of the FRG equation, namely,

$$\Gamma\_k = -\int \mathcal{L}\_{\text{eff}} = \int \text{d}k \, \frac{1}{2} \text{STr}(S^{(2)} + \mathcal{R}\_k)^{-1} \partial\_t \mathcal{R}\_k = \frac{1}{2} \text{STr} \, \log(S^{(2)} + \mathcal{R}\_k) + \text{const.} \tag{134}$$

Setting the bare action to the standard expression *S* = 1 4 *dx Fa μνFa μν*, which here corresponds to the UV limit of the effective theory, the integration constant can be fixed by requiring that Γ*k* vanishes for a vanishing field strength.

The inversion of the regularized propagator then requires the use of a harmonic gauge fixing, entering the action *S*gf and depending on the *α*-parameter, with associated Faddeev-Popov ghosts, in *S*gh, namely,

$$S\_{\rm gf} = \frac{1}{2\alpha} \int \mathrm{d}x \, D^{\mu} d^{a}\_{\nu} D^{\nu} a^{a}\_{\mu} \, \qquad \qquad \qquad S\_{\rm gh} = \int \mathrm{d}x \, D\_{\mu} \mathfrak{c}\_{\nu} D^{\mu} \mathfrak{c}^{\nu} \, \, \, \, \tag{135}$$

where the background methods have been deployed, and barred quantities are calculated with respect to background fields. Then, in the Landau gauge where *α* → 0, the super-trace recasts along the transverse sector as

$$\frac{1}{2}\text{STr}\,\log(S^{(2)} + R\_k) = \frac{1}{2}\text{Tr}\_{\text{trans}}\log\left[\bar{D}\_T^{\mu\nu} + \mathcal{R}\_k(\bar{D}\_T^{\mu\nu})\right] - \frac{1}{2}\text{Tr}\_{\text{gh}}\log\left[\bar{D}\_{\text{gh}}^{\mu\nu} + \mathcal{R}\_k(\bar{D}\_{\text{gh}}^{\mu\nu})\right],\tag{136}$$

where the differential operators can be expressed in terms of the SU(*N*) structure constant *f abc* and the YM coupling constant *g*YM as

$$
\vec{D}\_{\rm T}^{\mu\nu} = \vec{\Box} \,\delta^{\rm cb} \delta^{\mu\nu} + \mathcal{g}\_{\rm YM} \vec{F}^{a} \vec{\mu}^{\mu\nu} f^{\rm abc}, \qquad \vec{D}\_{\rm gh}^{\mu\nu} = \eta^{\mu\nu} \vec{\Box}. \tag{137}
$$

In order to disentangle the emergence of a condensate as a solution to the FRG equation, one restricts the consideration to the case of SU(2). Although this could seem to be limiting in the wider theoretical perspective, nonetheless, the restriction to SU(2) shall be considered as a selection of a subgroup of SU(*N*).

Bearing the discussion above in mind, one may select a straightforward expression for the regulator function in terms of a cut-off scale, *Rk*(D) = *k*2. The super-trace STr may then be evaluated using the Schwinger formula

$$\ln A \equiv \int\_0^\infty \frac{ds}{s} \left[ \mathcal{e}^{-sA} - \mathcal{e}^{-s} \right], \qquad A > 0,\tag{138}$$

and opting for a convenient choice of the background, the self-dual one, as investigated in Ref. [360]. These choices allow for the evaluation of the super-trace as a sum over the eigenvalues of the kinematic terms. In general, the eigenvalues of D¯ T are not known for an arbitrary background. However, those were found in the case of a self-dual background considered in Ref. [360], and a general discussion on the trace technology used in this approach can be found, e.g., in Refs. [358,359]. Applying the results of Ref. [360] in the case of a self-dual background dominated by the color-magnetic field *B*, Ref. [376] obtained an expression for the regularized effective action in the following form:

$$L\_{\rm eff} = \frac{g\_{\rm YM}^2 \theta}{4\pi^2} \int\_0^\infty \frac{\text{ds}}{\text{s}} \left[\varepsilon^{-As} - \varepsilon^{-s}\right] \left(\frac{1}{4\sinh^2(s)} + 1 - \frac{1}{4\text{s}^2}\right), \quad A = \sqrt{\frac{k^4}{g\_{\rm YM}^2 \theta}}, \quad \theta \equiv B^2 > 0. \tag{139}$$

This expression can be recast in terms of the order parameter in the CM domain, J = 2*g*2 YM*θ* > 0, which was introduced earlier and then analytically continued to the CE branch J < 0 as follows [66],

$$\mathcal{L}\_{\text{eff}} = \frac{2\mathcal{J}}{16\pi^2} \int\_0^\infty \frac{\text{ds}}{s} \left[ \varepsilon^{-sA\left[\mathcal{J}\right]} - \varepsilon^{-s} \right] \left[ \frac{1}{4\sinh^2 s} - \frac{1}{4s^2} + 1 \right], \qquad A = \sqrt{\frac{\lambda^4}{\mathcal{J}\tanh(\mathcal{J}/\lambda^4 \varepsilon)}}\tag{140}$$

for *ε* 1. Indeed, due to

$$\left. \mathcal{J} \tanh \left( \frac{\mathcal{J}}{\varepsilon \lambda^4} \right) \right|\_{\varepsilon \to 0} \to \left| \mathcal{J} \right| , \tag{141}$$

the transition between Equation (139) and Equation (140) becomes apparent. Performing the explicit expansion of the rightmost parenthesis above, the form matches that of the one-loop effective Lagrangian in Equation (91) for SU(2).

The results of Appendix D, valid at one-loop order, may now be compared to the allloop order results. In Figure 16 (left), we show a direct comparison of the one-loop with the all-loop effective Lagrangian for J > 0 (CE branch only). A unique non-trivial minimum is found in the non-perturbative region, 0 < J ∗ < *λ*4, and it is therefore identified with the CE condensate [66] whose values for one-loop and all-loop cases differ at a permille level and thus express a remarkable consistency of the one-loop approximation. The allloop running coupling may be straightforwardly extracted from the all-loop effective Lagrangian as

$$\left(\left(\bar{g}^2\right)^{-1} = \frac{2}{4\pi^2} \int\_0^\infty \frac{\mathrm{d}s}{\mathrm{s}} \left[e^{-sA\left[\mathcal{J}\right]} - e^{-s}\right] \left[\frac{1}{4\sinh^2 s} - \frac{1}{4s^2} + 1\right],\tag{142}$$

and this is shown in Figure 16 (middle) for the one-loop and all-loop cases. To find an expression for the *β*-function, one may study the RG equation given by Equation (86). Using

$$\frac{\mathrm{d}\bar{\mathcal{g}}^2}{\mathrm{d}\mathcal{J}} = -\left(\mathcal{g}^2\right)^2 \frac{\mathrm{d}}{\mathrm{d}\mathcal{J}} \left(\frac{1}{\mathcal{\bar{g}}^2}\right). \tag{143}$$

We may express *β* as

$$\begin{split} \frac{\beta}{\mathcal{S}^2} &= -2\mathcal{J} \frac{\mathbf{d}}{\mathbf{d}\mathcal{J}} \Big( \frac{1}{\mathcal{S}^2} \Big) \\ &= -2\mathcal{J} \frac{2}{4\pi^2} \left( -\frac{\mathbf{d}A}{\mathbf{d}\mathcal{J}} \right) \int\_0^\infty \mathrm{d}s \, e^{-sA[\mathcal{J}]} \Big[ \frac{1}{4\sinh^2 s} - \frac{1}{4s^2} + 1 \Big], \end{split} \tag{144}$$

where

$$\frac{\mathrm{d}A}{\mathrm{d}\mathcal{J}} = -\frac{1}{2\mathcal{J}}\sqrt{\frac{\lambda^4}{\mathcal{J}\,\mathrm{tanh}(\mathcal{J}/\lambda^4\varepsilon)}}\left(1 + \frac{\mathcal{J}}{\lambda^2\varepsilon}\frac{1-\mathrm{tanh}^2(\mathcal{J}/\lambda^2\varepsilon)}{\mathrm{tanh}(\mathcal{J}/\lambda^2\varepsilon)}\right). \tag{145}$$

The *β* over the running coupling at one-loop order compared to the corresponding all-loop quantity for SU(2) is displayed in Figure 16 (right). It is clear from the figures that the one-loop approximation accurately captures the main features of the pure YM theory as the differences from the corresponding all-loop quantities are negligible.

**Figure 16.** (**Left**) All- and one-loop effective Lagrangian of SU(2) as dependent on J /*λ*4. Plotted here for the CE branch with J > 0, it is clear that the one-loop result captures the main features of pure YM theory. In the non-perturbative regime shown (J < *λ*4), the two minima in the vicinity of J∗ > 0 coincide. Here, = 0.01 has been used. (**Middle**) The inverse running coupling *g*¯−2 of SU(2). The (dashed) one-loop result captures well the overall behavior of the running coupling in comparison to the all-loop coupling (solid). Here, = 0.01 has been used. (**Right**) *β* over the running coupling at one-loop order compared to the all-loop quantity for SU(2). The ratio of the *β*-function to the coupling coincides with the all-loop result at one-loop level for small J , while a discrepancy is seen away from zero. Here, = 10−5 has been used.

As outlined in Ref. [376], when considering the formation of the non-perturbative CE condensate as the YM system rolls down toward the minimum of the effective Lagrangian depicted in Figure 16 (left), the equation that dictates the evolution of this system in cosmological time is found from the continuity equation [376] (see also Ref. [377] for a more recent discussion)

$$
\rho\_{\rm YM} + 3\frac{d}{a}(\rho\_{\rm YM} + p\_{\rm YM}) = 0 \,. \tag{146}
$$

Accounting for only homogeneous CE YM fields (i.e., for *θ* = *E*2) in the EMT, for simplicity, the expressions of the energy and pressure densities of the YM system *ρ*YM, *p*YM can be written as functionals of the effective action [376]:

$$p\_{\rm YM} = -\mathcal{L}\_{\rm eff}(\theta) + 2\theta \mathcal{L}\_{\rm eff}'(\theta) \,, \qquad p\_{\rm YM} = \mathcal{L}\_{\rm eff}(\theta) - \frac{2}{3}\theta \mathcal{L}\_{\rm eff}'(\theta) \, , \tag{147}$$

where prime denotes the functional variation with respect to *θ*. In this case, Equation (146) transforms to

$$
\dot{\theta} \left( \mathcal{L}'\_{\rm eff} + 2\theta \mathcal{L}'\_{\rm eff} \right) + 4 \frac{\dot{a}}{a} \theta \mathcal{L}'\_{\rm eff} = 0 \,, \tag{148}
$$

in comoving coordinates. Equation (148) may be integrated, given that Leff(*θ*) is sufficiently well behaved, such that

$$
\sqrt{\theta} \mathcal{L}'\_{\rm eff}(\theta) = \mathcal{C} \, a^{-2} \,, \tag{149}
$$

where C denotes a coefficient of proportionality that is fixed by the initial conditions. Equation (149) may be solved by inversion, providing the dynamics of a SU(2) YM condensate in cosmological time. Note that Equation (147) is only valid for the dominant electric-field configurations. An analogical analysis of magnetic-field quantum configurations and their cosmological evolution in the effective Lagrangian approach at the one-loop level has been performed in Ref. [377]. Note, generic YM configurations contain both magnetic and electric components. We refer the reader to Section 4.4 above for a thorough discussion of those more generic configurations and their cosmological (real-time) dynamics.

#### **5. Cosmological Implications of Gauge-Fields Driven Inflation**

Gauge-fields driven inflation has been considered in several models, starting form the paper by Ford [378], in which a hypercharge cosmological inflationary scenario was envisaged. The theory that was taken into account included, besides the Einstein-Hilbert action, the action for a massive hypercharge field, individuated by the Lagrangian

$$L = \frac{1}{4} F\_{\mu\nu} F^{\mu\nu} + V(A^a A\_a) \,. \tag{150}$$

Anisotropies that are naturally present in the Maxwell tensor impose to consider a Bianchi type-I metric of the form

$$ds^2 = dt^2 - a^2(t)(dx^2 + dy^2) - b^2(t)dz^2.\tag{151}$$

Provided the expression for the energy-momentum tensor of the massive vector field

$$T\_{\mu\nu} = F\_{\mu\beta} F\_{\nu}^{\beta} - \frac{1}{4} \mathcal{g}\_{\mu\nu} F\_{a\beta} F^{a\beta} - \mathcal{g}\_{\mu\nu} V + 2V' A\_{\mu} A\_{\nu} \tag{152}$$

with prime denoting the functional derivative in *AαAα*, the Einstein equations were recast as

$$2\frac{d\vec{b}}{ab} + \frac{\vec{a}^2}{a^2} = 8\pi\epsilon \, , \qquad \qquad 2\frac{\vec{a}}{a} + \frac{\vec{a}^2}{a^2} = -8\pi p\_{z'} \, , \tag{153}$$

with as the energy density and *pz* as the pressure density component along the *z*-axis. The conservation law, derived assuming that the pressure density components along the *x* and *y* axes equal *px* = *py*, encodes

$$
\dot{\epsilon} + \left(2\frac{\dot{a}}{a} + \frac{\dot{b}}{b}\right)\epsilon + 2\frac{\dot{a}}{a}p\_x + \frac{\dot{b}}{b}p\_z = 0 \,, \tag{154}
$$

and finally, the only non-vanishing component of the vector field fulfills the equation

$$
\ddot{A}\_z + \left[2\frac{\dot{a}}{a} - \frac{\dot{b}}{b}\right] \dot{A}\_z + 2V'A\_z = 0 \,. \tag{155}
$$

The solutions to this system of equations, Equations (153)–(155), immediately outlined the impossibility to suppress anisotropies over the dynamical evolution of the background and the vector field. Thus, already with the seminal attempt of Ref. [378], it was evident that the over-production of anisotropies at the cosmological level would have provided a main issue, eventually ruling out this class of models.

A possible way to overcome the abundance of anisotropies, severely constrained by the CMB observations by WMAP and Planck, was imagined by Golovnev, Mukhanov and Vanchurin [379], who considered a stochastic distribution of *N* 1012 vector fields, randomly spanning the space directions. Each vector field was assumed to be massive and regulated by the action

$$S = \int d^4x \sqrt{-g} \left( -\frac{R}{16\pi} - \frac{1}{4} F\_{\mu\nu} F^{\mu\nu} + \frac{1}{2} (m^2 + \frac{R}{6}) A\_{\mu} A^{\mu} \right), \tag{156}$$

where *m* denotes the mass of the hypercharge vector field considered, such that *Fμν* = ∇*μAν* − ∇*νAμ* = *∂μAν* − *∂νAμ*.

The equations of motion for the gauge field, which in a covariant fashion recast

$$\frac{1}{\sqrt{-\mathcal{g}}} \frac{\partial}{\partial \mathbf{x}^{\mu}} \left( \sqrt{-\mathcal{g}} F^{\mu \nu} \right) + \left( m^2 + \frac{R}{6} \right) A^{\nu} = 0 \,, \tag{157}$$

split into a temporal "0" component, implying *A*0 = 0, and into space components, which finally provide the equation fro the "field strength" *Bi* ≡ *Ai*/*a*, i.e.,

$$
\ddot{B}\_i + 3HB\_i + m^2B\_i = 0 \,\,,
\tag{158}
$$

having introduced comoving coordinates such that *ds*2 = *dt*2 − *a*2(*t*)*dx*2, and these were denoted with dot derivative with respect to the cosmological time, so that *H* = *a*˙/*a*. For a homogeneous vector field, it was then noticed that the energy-momentum tensor can be expressed by

$$\begin{aligned} T^0\_{\ j} &= -\frac{1}{2} \left( \dot{B}^2\_k + m^2 \dot{B}^2\_k \right), \\ T^i\_{\ j} &= -\left[ -\frac{5}{6} \left( \dot{B}^2\_k - m^2 \dot{B}^2\_k \right) - \frac{2}{3} H \dot{B}\_k B\_k - \frac{1}{3} (\dot{H} + 3H^2) \dot{B}^2\_k \right] \delta^i\_j + \dot{B}^i \dot{B}\_j + H (\dot{B}^i B\_j \quad \text{(159)}\\ &+ \quad B\_j B^i) + (\dot{H} + 3H^2 - m^2) \dot{B}^i B\_j, \end{aligned} \tag{159}$$

where the summation over the index *k* is meant. Summing up now the contributions of a triplet of mutually orthogonal fields *B*(*a*) *i* with the same magnitude |*B*|, the total energy momentum tensor is averaged to the quantity

$$T\_0^0 = \epsilon = \frac{3}{2} \left(\mathcal{B}\_k^2 - m^2 B\_k^2\right), \qquad \qquad T\_j^i = -p\delta\_j^i = -\frac{3}{2} \left(\mathcal{B}\_k^2 - m^2 B\_k^2\right)\delta\_j^i,\tag{160}$$

where *Bk* satisfies

$$
\mathcal{B}\_i + 3H\mathcal{B}\_i + m^2 \mathcal{B}\_i = 0 \,. \tag{161}
$$

These latter relations for the energy-momentum tensor can be proved by the fact that for *B*(*a*) *i* , it holds:

$$\sum\_{i} B\_{i}^{(a)} B\_{(b)\ i} = |B|^{2} \delta\_{(b)}^{(a)} \to \sum\_{a} B\_{j}^{(a)} B^{(a)} \, ^{i} = |B|^{2} \delta\_{j}^{i} \,. \tag{162}$$

Summing up over a large amount of *N* triple of fields, the energy-momentum tensor finally acquires the expression

$$T\_0^0 = \mathfrak{e} \simeq \frac{N}{2} \left( \dot{B}\_k^2 + m^2 B\_k^2 \right) \qquad \qquad T\_{\vec{\eta}} \propto \sum\_{a=1}^N B\_i^{(a)} B\_{\vec{\eta}}^{(a)} \simeq \frac{N}{3} B^2 \delta\_{\vec{\eta}}^i + O(1) \sqrt{N} B^2. \tag{163}$$

Within this scenario, anisotropies are then proven to fall off as 1/ √*N*, hence justifying consistency with the experimental data.

Finally, accounting for an inflationary slow-roll phase, with *B*˙ *i* 0, one finally finds for the first Friedmann equation,

$$H^2 = \frac{8\pi}{3}\epsilon \simeq \frac{4\pi}{3}Nm^2B^2. \tag{164}$$

A different perspective was suggested by Maleknejad and Sheikh-Jabbari, who imagined in Refs. [380–382] that it could become relevant at the level of cosmic perturbations, while at the same time driving the inflationary background evolution of the universe. Isotropy could be guaranteed here by aligning the internal indices (of the adjoint representation) of a SU(2) subgroup of SU(*N*) to the space directions, namely, for the space component of the connection *A*

$$A^a\_i = \psi(t)\epsilon^a\_i = \psi(t)a(t)\delta^a\_i \, , \qquad \qquad \mathfrak{e}^a = a(t)\delta^a\_i \, , \tag{165}$$

where *ψ*(*t*) is a scalar field, *ea i* denotes the triads, which recombine into the space metric *hij* = *ea i ea j* and are expressed here in the comoving coordinates. For compactness of notation, the authors reshuffled *φ*(*t*) = *ψ*(*t*)*a*(*t*).

Specifically, the authors considered a YM action provided with the square of a Pontriagin term, i.e.,

$$S = \int d^4x \sqrt{-g} \left[ -\frac{R}{2} - \frac{1}{4} F^a\_{\mu\nu} F^{\mu\nu}\_a + \frac{\kappa}{384} \left( \varepsilon^{\mu\nu\lambda\sigma} F^a\_{\mu\nu} F^a\_{\lambda\sigma} \right)^2 \right],\tag{166}$$

having set 8*πG* = 1, denoted with *εμνλσ* the totally antisymmetric tensor, and chosen the real parameter *κ* to be positive.

Labeling with a subscript "YM" the contributions arising from the YM terms, and with *κ* the ones related to the Pontriagin terms, it was found that

$$
\epsilon\_{\rm YM} = \frac{3}{2} \left( \frac{\dot{\phi}^2}{a^2} + \frac{g^2 \phi^4}{a^4} \right), \qquad \qquad \qquad \epsilon\_\kappa = \frac{3}{2} \frac{\kappa g^2 \phi^4 \dot{\phi}^2}{a^6}, \tag{167}
$$

where *g* is the YM coupling constant entering the definition of the field strength of the SU(2) connection

$$F^{a}\_{\mu\nu} = \partial\_{\mu}A^{a}\_{\nu} - \partial\_{\nu}A^{a}\_{\mu} - \g{\epsilon}^{abc}A^{b}\_{\mu}A^{c}\_{\nu} \,. \tag{168}$$

and where = YM + *κ* and *p* = 1 3 YM − *κ*. Having provided these definitions, it was possible to show the emergence of a slow-roll phase of inflation, which was dominated by the *κ*-term contribution *κ* over the YM contribution YM, namely *κ*  YM. Furthermore, within the scenario introduced in Refs. [380–382], it was conceivable that sizable effects could be turned on, so as to source primordial magnetic fields at the cosmic perturbation level, with specific possible observational features on the CMB and primordial cosmic magnetic fields outlined.

A novel framework with respect to the ones so far discussed was elaborated by Alexander, Marcianò and Spergel in Ref. [383], in which the authors considered the interaction a model of inflation with as a new ingredient the interaction of an Abelian gauge field with a fermionic charge. This scenario then dramatically differs from only adopting gauge fields in order to generate a realistic inflationary epoch. As a by-product of this approach, researchers considered the possibility of generating the a net-lepton asymmetry. The Sakharov conditions are realized in the model presented in Ref. [383] during the inflationary epoch, due to a dynamical inter-change of the gauge field fluctuations into the lepton asymmetry of the universe.

The action the authors moved from involved, as well as a U(1) hypercharge field *Aμ*, also a massive scalar field *θ*, interacting with *Aμ* through a Chern-Simons term, i.e.,

$$\begin{split} S\_{-}&=S\_{D}+\int d^{4}x\sqrt{-g}\left[\frac{M\_{\rm PL}^{2}R}{8\pi}-\frac{1}{2}\partial\_{\mu}\theta\partial^{\mu}\theta-\frac{1}{4}F\_{a\beta}F^{a\beta}+\frac{\theta}{4M\_{\ast}}F\_{a\beta}\bar{F}^{a\beta}\right],\\ S\_{D}&=\int d^{4}x\sqrt{-g}\left(\imath\bar{\Psi}\gamma^{\mu}\nabla\_{\mu}\Psi+\mathfrak{c}\mathfrak{c}+M\bar{\Psi}\Psi+\eta\bar{\Psi}\gamma^{\mu}\psi A\_{\mu}\right),\end{split} \tag{169}$$

where *M*PL denotes the Planck mass, *M*∗ is the mass-scale of the pseudo-scalar decay constant, regulating with *theta* the CP-violating Chern-Simons such as term, *Fμν* = *∂*[*μAν*] is the field strength of the U(1) connection *A*, *F*˜*αβ* = *εαβλσF*˜ *λσ* is the gravitational Hodge-dual, and *γμ* = *e μ I γI* , with *γI* Dirac matrices, *e μ I* inverse tetrad and internal indices *I* = 0, . . . 3.

The model, which is then based on the interaction between a homogeneous and isotropic configuration of a U(1) gauge field and a fermionic charge density J0, relies on the regulated fermionic charge density as generated from a Bunch-Davies vacuum state, using the procedure outlined by Koksma and Prokopec in Ref. [384]. In conformal coordinates, this was found to redshift as 1/*a*(*η*). Then, within the scenario of Ref. [383], the time-like component of the hypercharge gauge field was found to be sourced by the fermionic charge, consistently with a growth in the gauge field proportional to the scale factor, namely,

$$A\_0(\eta) \sim a(\eta)\,. \tag{170}$$

This motivates results with inflation dominated by the energy density stored within the interaction among the gauge field and the fermionic charge, namely *A*0J0, which is approximately constant over the inflationary epoch. The appealing feature of Ref. [383] stands in the possibility to obtain an epoch of cosmic inflation involving the physical description fields already existing in nature, specifically the time-like U(1) gauge field interacting with a fermionic charge density. Nonetheless, the role of the scalar field cannot be underestimated, retaining a certain relevance in producing baryogenesis, and providing a graceful exit from inflation. Indeed, the mechanism that accounted for the graceful exit is strictly interconnected to the one advocated for reproducing the baryogenesis, with the right baryon asymmetry index. The Chern-Simons term, through the coupling to the pseudo-scalar field, converts gauge field fluctuations into lepton number, while the rapid oscillation of the pseudo-scalar field near its minimum allows achieving thermalization of the gauge field and thus to end inflation. The relevance of the coupling between scalar modes, there interpreted as axions, was further investigated in Ref. [385].

An improvement of the scheme first addressed in Ref. [383] was provided in Ref. [386], where the authors analyzed the consistency of the model via the Stückelberg mechanism [387]. This provided an incorporation of the longitudinal scalar DoFs into the hypercharge field, which is now massive. This could be thought again as a further stepforward, toward the realization of an inflationary mechanism relying on the YM dynamics, the hypercharge sector being eventually recognized as an Abelian subgroup of the SU(*N*) gauge sector.

The action of the theory was then considered to be

$$S = \int d^4x\sqrt{-g}\left[\frac{M\_{\rm PL}^2R}{8\pi} - \frac{1}{4}G\_{a\beta}G^{a\beta} - \frac{1}{2}m^2\mathcal{C}\_{\mu}\mathcal{C}^{\mu} + \mathcal{C}\_{\mu}\mathcal{I}^{\mu} + L\_{\rm D}\right]$$

$$L\_{\rm D} = -i\Psi\gamma^{\mu}\nabla\_{\mu}\Psi + c.c. + M\Psi\psi,\tag{171}$$

having introduced the massive Stückelberg field *Cμ* = *Aμ* − 1 *m ∂μθ* and its field strength *Gμν* = *∂*[*μCν*] = *∂*[*μAν*] = *Fμν*, and the fermionic vector current J *μ* = *qψγ*¯ *μψ*. Gauge invariance is ensured in this framework by the transformations

$$A\_{\mu} \to A\_{\mu}^{\prime} = A\_{\mu} + \partial\_{\mu} \Lambda \,, \qquad \qquad \theta \to \theta^{\prime} = \theta + m\Lambda \,. \tag{172}$$

Within this framework, the authors could prove the existence and the stability of dynamical attractor solutions for the cosmological inflation epoch, which is again driven by the coupling among the fermions and a (massive) gauge field. Numerical analyses then showed that stability is attained for a large basin of the initial conditions, making this inflationary scenario almost independent on these latter: inflation arises without fine tuning and without the need of postulating any effective potential or any non-standard coupling.

An alternative scenario featuring a coupling between an axion-like field and an SU(2) gauge field is known as the chromo-natural inflation [388]. The rotationally invariant homogeneous condensate of the gauge field satisfies an attractor solution that enables it to drive cosmic inflation for the axion decay constant having a natural value at a sub-Planckian scale. Interestingly enough, this scenario features a possibility for termination of inflation as soon as the axion potential vanishes, simultaneously providing a small tensor-to-scalar perturbation ratio.

An inflationary scenario, taking into account non-trivial topological features deployed in Refs. [389–391], was extended in Ref. [392] so as to account, over the universe expansion, for a sector of a strongly coupled QCD-like gauge theory. The idea at the base of investigations in Refs. [389–391] is to perform a periodic (between the Σ and Σ surfaces) path integral over Euclidean geometries,

$$e^{-\Gamma} = \int\_{\mathbb{S},\mathfrak{\phi}|\_{\Sigma} = \mathfrak{g}, \mathfrak{\phi}\_{\Sigma'}} D[\mathfrak{g}, \mathfrak{\phi}] e^{-S\_{E,\mathfrak{F}}\mathfrak{\phi}},\tag{173}$$

with *g* and *φ*, respectively, representing the metric and matter field and *SE* denoting the Euclidean action, so as to extract the 'density matrix' of the universe,

$$\rho[\boldsymbol{\varphi}, \boldsymbol{\varphi}'] = e^{\Gamma} \int\_{\mathcal{S}, \Phi \vert\_{\Sigma, \Sigma'}} D[\boldsymbol{\mathcal{g}}, \boldsymbol{\Phi}] \boldsymbol{\varepsilon}^{-S\_E[\boldsymbol{\mathcal{g}}, \Phi]} \; \; \; \tag{174}$$

which describes a microcanonical ensemble, *ϕ* denoting field configurations that encode both gravitational and matter variables.

The uniform distribution over the Euclidean spaces actually corresponds, over Lorentzian spaces, to a distribution that is peaked about complex saddle points of the path integral. The latter can be then represented by cosmological instantons, entailing a bounded range values for the cosmological constant.

On the other hand, inflationary cosmologies can be engendered by the very same instantonic solutions [389–391]. The low energy of the accelerated expansion can be then attained at its late stage, resorting to the dynamical evolution of extra dimensions specifically postulated in string theory framework [391]. This results in a bounded range for the very early (inflationary) cosmological constant, which provides a constraint on the available landscape of the string vacua. Finally, the same mechanism can be advocated to give rise to a possible DE candidate, accounting for the quasi-equilibrium decay of the microcanonical state of the universe. Within this scenario, Barvinsky and Zhitnitsky promoted a new picture for the emergence of an inflationary spacetime [392], resorting to considerations developed in Refs. [14,393,394] on the generation, in a strongly coupled QCD-like gauge theory, of the vacuum energy from non-trivial topological features.

The limits of the usual semi-classical expansion were overcome by the dominant contribution of the numerous conformal modes. Integration over the modes then provides the quantum effective action of the conformal field theory ΓCFT[*gμν*], which can be calculated with methods similar to those ones implemented in determining the conformal anomaly. Starting from the FLRW background, accounting for a periodic factor *a*(*τ*)—this is due to the fact that functions of the Euclidean time are supported to the circle S1—and finally using a local conformal transformation to the static Einstein universe and the very same well-known trace anomaly, one finds

$$g\_{\mu\nu} = \frac{\delta\Gamma\_{\rm CFT}}{\delta g\_{\mu\nu}} = \frac{1}{4(4\pi)^2} g^{1/2} \left( \alpha \square R + \beta E \gamma \mathcal{C}\_{\mu\nu\rho\sigma}^2 \right),\tag{175}$$

where we introduced the Gauss-Bonnet term *E* = *R*2 *μναγ* − 4*R*2 *μν* + *R*2 and the Weyl tensor *Cμνρσ*.

Considering a spacetime with topology S3 × S1, and moving from the expression of the energy of the gauge field holonomy, winding across the compactified coordinate of the length T , Barvinsky and Zhitnitsky found that

$$\rho = \rho\_{\text{vac}}[\mathbb{S}^3 \times \mathbb{S}^1] - \rho\_{\text{vac}}[\mathbb{R}^4] = \frac{\mathbb{E}\_{\mathcal{T}} \Lambda\_{\overline{\text{QCD}}}^{\overline{\text{3}}}}{\mathcal{T}},\tag{176}$$

with Λ*QCD* being the scale of the QCD-like gauge theory, *c*¯T being a dimensionless constant of order one, and similarly, the full period of the proper Euclidean time on these periodic *m*-fold garland instantons is given by the analogous integral,

$$\mathcal{T} = \oint\_{\mathbb{S}\_1} d\mathbf{r} \,\mathrm{s}\tag{177}$$

which in an FLRW metric background reduces to the 2*m*-multiple result

$$
\mathcal{T} = 2m \int\_{a\_-}^{a^+} \frac{da}{a} \, , \tag{178}
$$

where the integral is between the two neighboring turning points of *a*(*τ*) such that *a*˙(*τ*±) = 0.

#### **6. Summary**

In this review, we have made a brief outlook of the current status of confined and de-confined QCD dynamics in the early universe as well as the key methodology for studies QCD in the strongly coupled regime relevant for cosmological evolution. The covered research areas are broadly inter-disciplinary, and our discussion may not be fully exhaustive. Still, we have identified a few quite unexpected and intriguing connections between currently pursued research in particle physics and possible dynamics of the early universe. Such fundamental questions as the gauge-fields-driven inflation, cyclic universe, particle production mechanisms, non-perturbative real-time dynamics of the QCD ground state, a rather challenging problem of dynamical generation of cosmological DE and DM, the structure of the QCD vacuum, the QCD phase transitions and the role of QCD matter in late-time universe evolution are among the key points of this review. Such a wide breadth of topics, with deep roots into QCD or, more generically, quantum YM field theories, exhibits enormous and critical significance of microscopic dynamics of particle physics and confined field theories for understanding of the macroscopic cosmic evolution. The picture is far from its final shape though, and many more pillars of such connections and possible interplay are yet to be established. We believe that our review can be useful for both young researchers and for more senior experts specialized in both particle physics and cosmology research areas.

**Author Contributions:** Conceptualization, A.A., A.M., T.L., R.P. and M.Š.; validation, R.P., T.L. and M.Š.; formal analysis, R.P., T.L. and M.Š.; investigation, R.P., T.L. and M.Š.; resources, R.P. and M.Š.; writing—original draft preparation, A.M., T.L., R.P. and M.Š.; writing—review and editing, A.A., A.M., T.L., R.P. and M.Š.; supervision, R.P.; project administration, R.P.; funding acquisition, R.P. and M.Š. All authors have read and agreed to the published version of the manuscript.

**Funding:** Work of A.A. is supported by the Talent Scientific Research Program of College of Physics, Sichuan University, Grant No.1082204112427 & the Fostering Program in Disciplines Possessing Novel Features for Natural Science of Sichuan University, Grant No. 2020SCUNL209 and 1000 Talent program of Sichuan province 2021. A.M. wishes to acknowledge support by the Shanghai Municipality, through the grant No. KBH1512299, by Fudan University, through the grant No. JJH1512105, the Natural Science Foundation of China, through the grant No. 11875113, and by the Department of Physics at Fudan University, through the grant No. IDH1512092/001. R.P. is supported in part by the Swedish Research Council grant, contract number 2016-05996, as well as by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 668679). M.Š. is partially supported by the grants LTT17018 and LTT18002 of the Ministry of Education of the Czech Republic.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Elements of Relativistic Hydrodynamics**

Equations of relativistic hydrodynamics are based on conservation of the energymomentum and the current

$$
\partial\_{\mu}T^{\mu\nu} = 0 \,, \quad \partial\_{\mu}j\_i^{\mu} = 0 \,, \tag{A1}
$$

where *j μ i* , *i* = *B*, *Q*, *L*, ... denotes the conserved currents corresponding to baryon number *B*, electric charge *Q*, lepton number *L*, etc. Both *Tμν* and *j μ i* can be decomposed into time-like

and space-like components using natural projection operators, the local flow four-velocity *uμ*, and the second-rank tensor perpendicular to it Δ*μν* = *gμν* − *uμuν* [11,25,395,396]:

$$T^{\mu\nu} = \varepsilon\mu^{\mu}\mu^{\nu} - p\Delta^{\mu\nu} + \mathsf{W}^{\mu}\mu^{\nu} + \mathsf{W}^{\nu}\mu^{\mu} + \pi^{\mu\nu},\tag{A2}$$

$$j\_i^{\mu} \quad = \quad n\_i \mu^{\mu} + V\_i^{\mu} \tag{A3}$$

where = *uμTμνuν* is the energy density, *p* ≡ *ps* + Π = −1 3Δ*μνTμν* is the total (hydrostatic *ps* + bulk Π) pressure, *Wμ* = Δ*μ αTαβuβ* is the energy (or heat) current, *ni* = *uμj μ i* is the charge density, *Vμ i* = Δ*μ ν j ν i* is the charge current, and *πμν* = *Tμν* is the shear stress tensor. The angular brackets in the definition of the shear stress tensor *πμν* stand for the following operation:

$$
\langle A^{\mu \nu} \rangle = \left[ \frac{1}{2} (\Delta^{\mu}\_{\ a} \Delta^{\nu}\_{\ \beta} + \Delta^{\mu}\_{\ \beta} \Delta^{\nu}\_{\ \ a}) - \frac{1}{3} \Delta^{\mu \nu} \Delta\_{a \beta} \right] A^{a \beta} \,. \tag{A4}
$$

To further simplify our discussion, we restrict ourselves in the following to only the one conserved charge, the baryon number *B*, and denote the corresponding baryon current as *j μ* ≡ *j μ B*. The various terms appearing in the decompositions (A2) and (A3) can then be grouped into ideal and dissipative parts as follows

$$T^{\mu\nu}\_{\ \nu} = \left. T^{\mu\nu}\_{\ \dot{\alpha}} + T^{\mu\nu}\_{\ \dot{\alpha}} = \left[ \epsilon \mu^{\mu} \mu^{\nu} - p\_s \Delta^{\mu\nu} \right]\_{\ \dot{\alpha}} + \left[ -\Pi \Delta^{\mu\nu} + \mathcal{W}^{\mu} \mu^{\nu} + \mathcal{W}^{\nu} \mu^{\mu} + \pi^{\mu\nu} \right]\_{\ \dot{\alpha}s} \tag{A5}$$

$$j^{\mu} \quad = \quad j\_{\rm id}^{\mu} + N\_{\rm dis}^{\mu} \quad = [nu^{\mu}]\_{\rm id} + [V^{\mu}]\_{\rm dis} \,. \tag{A6}$$

Neglecting the dissipative parts, the energy-momentum conservation and the current conservation (A1) define ideal hydrodynamics. In this case (and for a single conserved charge), a solution of the hydrodynamical Equation (A1) for a given initial condition describes the spacetime evolution of the six variables—three state variables (*x*), *p*(*x*), *n*(*x*), and three space components of the flow velocity *uμ*. However, since (A1) constitutes only five independent equations, the sixth equation relating *p* and , the EoS *p*(), has to be added by hand in order to solve them.

Two definitions of flow can be found in the literature, see e.g., Refs. [11,395]; one related to the flow of conserved charge (Eckart):

$$
\mu\_E^{\mu} = \frac{j^{\mu}}{\sqrt{j\_{\nu} j^{\nu}}} \, \, \, \, \tag{A7}
$$

the other related to the flow of energy (Landau):

$$u\_L^{\mu} = \frac{T^{\mu}\_{\ \nu} u\_L^{\nu}}{\sqrt{\mu\_L^{\alpha} T\_{\alpha}^{\beta} T\_{\beta \gamma} u\_L^{\gamma}}} = \frac{1}{\varepsilon} T^{\mu}\_{\ \nu} u\_L^{\nu}. \tag{A8}$$

Let us note that *Wμ* = 0 (*Vμ* = 0) in the Landau (Eckart) frame. In the case of vanishing dissipative currents, both definitions represent a common flow. The Landau definition is more suitable when describing the evolution of matter at zero chemical potential, such as in the case of the mid-rapidity particle production in ultra-relativistic HIC at the LHC and at the top RHIC energy, or in the early universe. In this case, all momentum density is due to the flow of energy density, *uμTμν id* = *uν* and *uμTμν dis* = 0, i.e., the heat conduction effects can be neglected.

#### **Appendix B. Hydrodynamical Description of Dissipative Effects**

In its modern formulation, relativistic fluid dynamics provides an effective description of a system that is in local thermal equilibrium, and it can be derived from the underlying kinetic description through Taylor expansion of the entropy four-current *Sμ* = *suμ* in gradients of the local thermodynamic variables [25]. In zeroth order in gradients, one obtains ideal fluid dynamics

$$
\partial\_{\mu}S^{\mu} = \partial\_{\mu}(su^{\mu}) = u^{\mu}\partial\_{\mu}s + s\partial\_{\mu}u^{\mu} = 0 \,,\tag{A9}
$$

and the evolution of the scale factor of the universe is driven solely by the entropy conservation *s*(*t*)*a*3(*t*) = const. The higher orders describe effects due to irreversible thermodynamic processes such as the frictional energy dissipation between the fluid elements that are in relative motion or their heat exchange with its surroundings on its way to approach thermal equilibrium with the whole fluid.

When solving the hydrodynamic equations with the dissipative terms, it is customary to introduce the following two phenomenological definitions (so-called *constitutive equations*) for the shear stress tensor *πμν* and the bulk pressure Π appearing in Equation (A5) [395],

$$
\pi^{\mu\nu} = 2\eta \left< \nabla^{\mu} u^{\nu} \right> , \quad \Pi = -\mathcal{J}\partial\_{\mu} u^{\mu} = -\mathcal{J}\nabla\_{\mu} u^{\mu} \, , \tag{A10}
$$

where the angular brackets ... are defined in Equation (A4) and ∇*μ* = Δ*μν∂ν*. Neglecting the charge current *Vμ* in Equation (A6), the first-order expansion of *Sμ* is completely determined by the *shear viscosity η* and *bulk viscosity ζ* coefficients [395]:

$$T\partial\_{\mu}S^{\mu} = \pi\_{\mu\nu}\langle\nabla^{\mu}u^{\nu}\rangle - \Pi\partial\_{\mu}u^{\mu} = \frac{\pi\_{\mu\nu}\pi^{\mu\nu}}{\eta} + \frac{\Pi^2}{\overline{\zeta}} = 2\eta\langle\nabla^{\mu}u^{\nu}\rangle^2 + \zeta(-\partial\_{\mu}u^{\mu})^2 \ge 0. \tag{A11}$$

A well-known example of the flow involving both coefficients *η* and *ζ* is provided by boost-invariant one-dimensional expansion with the velocity in the *z* direction, *vz*, proportional to *z* co-ordinate [397]

$$
\mu\_{\rm Bj}^{\mu} = \frac{\mathbf{x}^{\mu}}{\mathbf{\bar{\tau}}} = \frac{\mathbf{t}}{\mathbf{\bar{\tau}}} \left( 1, 0, 0, \frac{\mathbf{z}}{\mathbf{t}} \right), \ \mathbf{\bar{\tau}} = \sqrt{t^2 - z^2}. \tag{A12}
$$

After inserting this solution into the constitutive Equation (A10), we arrive at the equation of motion [395]:

$$\frac{d\varepsilon}{d\tau} = -\frac{\epsilon + p\_s}{\tau} \left( 1 - \frac{4}{3\tau T} \frac{\eta}{s} - \frac{1}{\tau T} \frac{\zeta}{s} \right). \tag{A13}$$

The last two terms on the right-hand side in Equation (A13) describe a compression of the energy density due to viscous corrections. Two dimensionless coefficients in the viscous correction, *η*/*s* and *ζ*/*s*, where *s* is the entropy density, reflect the intrinsic properties of the fluids. It is worth mentioning that neglecting *η* and *ζ* in Equation (A13), i.e., for the ideal fluid EoS with *ps* = 1 3 , one obtains the celebrated Bjorken solution of ideal hydrodynamics [397]

$$\frac{d\epsilon}{d\tau} = -\frac{\epsilon + p\_s}{\tau} = -\frac{4}{3}\frac{\epsilon}{\tau} \qquad \Rightarrow \qquad \epsilon = \tau^{-4/3},\tag{A14}$$

frequently used to discuss the salient features of the ultra-relativistic HICs.

The one-dimensional character of the Bjorken flow (A12) makes it possible to replace the *z* co-ordinate with the radial one *r*. Radial flow in the transverse direction, i.e., when *r*⊥ = *x*2 + *y*2, was studied in Ref. [398]. For this case and the constant sound velocity *cs*, analytic solutions of relativistic viscous hydrodynamics describing expanding fireballs were developed in Ref. [399]. In a three-dimensional case, a new class of exact fireball solutions of relativistic dissipative hydrodynamics for arbitrary shear and bulk viscosities, as well as for other dissipative coefficients, was studied in Ref. [400]. The common property of these solutions is the presence of the relativistic Hubble flow.

However, the analogy between the solutions describing HICs and expansion of the early universe must not be pushed too far, since in the latter case form of the energymomentum tensor *Tμν* and particle four-current *jμ* of the matter, cf. Equations (A5) and (A6) is strongly constrained by the symmetries of the FLRW metric (38). In particular, due to the local momentum isotropy, the term ∇*μuν* appearing in the viscous shear-stress tensor *πμν*, cf. Equation (A10), vanishes [401]. Consequently, the term proportional to *η* in Equation (A13) disappears. The shear viscosity *η* also disappears in the theories with scalar perturbations of the metric tensor *gμν* [84,270]. There, the fluctuations of energy density destroy the homogeneity but not the isotropy of the early universe FLRW metrics. Hence, in the following discussion, we will consider mainly the bulk viscosity—the property of expanding matter arising typically in mixtures. They can be either of different species, as in a radiative fluid, or of the same species but with different energies, as in a Maxwell–Boltzmann gas. In each of these instances, the bulk viscosity provides the internal 'friction' that sets in due to the different cooling rates in the expanding mixture [401].

The relativistic Navier-Stokes description given by Equation (A10) accounts only for terms that are *linear* in velocity gradient. This leads, unfortunately, to severe problems. In particular, when the thermodynamic force ∇*μuν* or ∇*μuμ* is suddenly switched off/on, the corresponding thermodynamic flux *πμν* or Π which is a purely local function of the velocity gradient also instantaneously vanishes/appears [402]. The linear proportionality between dissipative fluxes and forces causes an instantaneous (acausal) influence on the dissipative currents, leading to numerical instabilities [403]. The solution of this problem requires the inclusion of terms that are second order in gradients [404]. The resulting equations for the dissipative fluxes *πμν* and Π then become the relaxation-type equations [395,405]. The latter encode the time delay between the appearance of thermodynamic gradients that drive the system out of local equilibrium and the associated build-up of dissipative flows in response to these gradients, thereby restoring causality [405]. Accounting for non-zero relaxation times at all stages of the evolution constrains departures from local equilibrium, thereby both stabilizing the theory and improving its quantitative precision.