*Editorial* **Shedding Light to the Dark Sides of the Universe: Cosmology from Strong Interactions**

**Roman Pasechnik 1,\* and Michal Šumbera <sup>2</sup>**


**\*** Correspondence: roman.pasechnik@thep.lu.se

The basic aim of this Special Issue was to reflect upon the modern status of research on strong interactions and their implications in Cosmology. In addtion to a few important original research articles, our Special Issue comprises a collection of reviews on existing implications of the physics of strongly coupled and non-perturbative phenomena such as those in Quantum Chromodynamics (QCD) to the universe evolution in different epochs. Such implications concern, in particular, cosmological aspects of quark–gluon plasma and phase transition dynamics, the physics of neutron stars, QCD equations of state (EoS), as well as the origins of Dark Matter (DM) and Dark Energy (DE) from both dark sectors and QCD, inflationary cosmology, etc.

A broad review [1] by the Guest Editors and their collaborators provides a comprehensive outlook of the key research results in the field of confined and de-confined QCD dynamics and their implications in physics of the early Universe. In addition, it briefly covers the basic methodology for studies of quantum field theories in the strongly coupled regime on the non-stationary background of the expanding Universe. A few rather important connections between currently pursued research in particle physics and possible dynamics of the early Universe have been identified and elaborated upon, potentially addressing such fundamental questions as particle production mechanisms in the early Universe, the origins of cosmic acceleration, and the non-perturbative real-time dynamics of the QCD ground state, among others. Given the broadly inter-disciplinary coverage of a variety of different and challenging topics, this review represents an important, but not exhaustive, reference for frontier research at an intersection of particle physics and cosmology.

The review [2] by Vitaly Beylin, Maxim Khlopov, Vladimir Kuksa, and Nikolay Volchanskiy discusses some of the basic cosmological effects of strongly coupled New Physics focusing on the possible nontrivial role of strong interactions. One particular example is the presence of new stable colored particles, such as exotic quarks, which could give rise to a composite DM candidate. In addition, an overview of new stable DM composite candidates was given in the context of QCD-like interactions in various scenarios of New Physics, broadly referred to as Techni (or Hyper)-Color. A particular emphasis was given on possible interactions of new stable particles with those in the Standard Model, with interesting implications for cosmic-ray and high-energy neutrino astrophysics and for the phenomenology of stable fractionally charged particles.

Another review [3] by Ralf Hofmann overviews possible cosmological consequences of a scenario when the SU(2) gauge principle governs the Cosmic Microwave Background (CMB) instead of U(1) one attempting to address the existing tensions between local and global cosmological measurements in the framework of ΛCDM. These consequences concern possible changes in the radiation and dark sector components of the Universe, as well as in the structure formation assisted by the de-percolation of condensed ultralight axion configurations with the Peccei–Quinn scale close to the Planck mass. It has been demonstrated that a compatibility of the axionic field profiles with typical DM galactic

**Citation:** Pasechnik, R.; Šumbera, M. Shedding Light to the Dark Sides of the Universe: Cosmology from Strong Interactions. *Universe* **2022**, *8*, 545. https://doi.org/10.3390/ universe8100545

Received: 9 October 2022 Accepted: 14 October 2022 Published: 19 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

halos emerges if cosmological Yang–Mills fields confining at much higher energy scales than that of SU(2)CMB are responsible for the generation of the axion mass.

An extensive review [4] by G. Fiorella Burgio and Isaac Vidaña considers possible correlations between astrophysical observables of neutron stars (NS) and the properties of atomic nuclei. This review focuses in particular on the major role of the EoS of nuclear matter in defining the main characteristics of NS, including the existing experimental constraints. Such aspects of the NS phenomenology as the NS tidal deformability and its correlations with the stellar radius, the stiffness of the symmetry energy, and the neutronskin thickness have been discussed.

The article [5] by Petr Jizba, Lesław Rachwał, Stefano G. Giaccari, and Jaroslav K ˇnap addresses the issue of a dynamical breakdown of scale invariance in quantum Weyl gravity (QWG), together with the related cosmological implications. In their study, which is based on the approach of functional renormalization group, the authors investigate the infrared (IR) physics of the quantum Weyl gravity and find it to be surprisingly rich and interesting, in particular in connection with the structure of its IR fixed point and the corresponding cosmological consequences. In a cosmology that is implied by the broken phase of QWG, they are able to map the broken-phase's effective action on a two-field hybrid inflationary model that, in its low-energy phase, approaches the Starobinsky *f*(*R*) model with a gravi-cosmological constant that has a negative sign in comparison to the usual matter-induced cosmological constant. The implications of this finding for cosmic inflation are also discussed.

In the article [6] by Roland Kirschner and George Savvidy, it has been considered a new possibility that inside hadrons inhabitated by standard SU(3)*<sup>c</sup>* partons—quarks and gluons—there are additional partons–tensorgluons, which can carry a part of the proton momentum. Since the gluons do not couple directly to the photon in deep-inelastic scattering measurements, their density inside the hadrons is one of the least constrained functions. Here comes the opportunity for the tensorgluons: although their existence does not predict a new hadronic state, it leads to a modification of the parton distribution functions of a proton. Moreover, because tensorgluons have a larger spin than ordinary gluons, they can influence the spin structure of the nucleon.

The study [7] by Andrea Addazi, Stephon Alexander and Antonino Marcianò addresses an important fundamental issue of origin of the late-time cosmological acceleration due to the possible existence of an addition dark QCD-like matter sector. Using the arguments of strong dynamics, such as the formation of dark gluon and dark quark condensates breaking the chiral symmetry in the dark sector, the authors draw a conclusion that the interaction energy between the dark condensates may cause late-time cosmic acceleration, reproducing the observable effect of the cosmological constant.

Last but not least, the article [8] by Dmitri N. Voskresensky is devoted to evolution of quasiperiodic structures in a non-ideal hydrodynamic description of phase transitions. It starts with a general introduction to first- and second-order phase transitions. The latter could have taken place either in the early universe, in the course of heavy-ion collisions and supernova explosions, and also in proto-neutron stars, in cold compact stars, and in the condensed matter at terrestrial conditions. The author presents some novel solutions of non-ideal hydrodynamics describing the evolution of quasiperiodic structures that are formed in the course of the phase transitions. The most important result of this work concerns the finding that viscosity and thermal conductivity are the driving forces of the first-order liquid–gas and quark–hadron phase transitions to the state characterized by the zeroth wave number and by the instability occurring for temperatures below the isothermal spinodal region.

**Funding:** This Guest Editors' activity received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


### *Review* **Cosmology from Strong Interactions**

#### **Andrea Addazi 1,2, Torbjörn Lundberg 3, Antonino Marcianò 2,4, Roman Pasechnik 3,\* and Michal Šumbera <sup>5</sup>**


**Abstract:** The wealth of theoretical and phenomenological information about Quantum Chromodynamics at short and long distances collected so far in major collider measurements has profound implications in cosmology. We provide a brief discussion on the major implications of the strongly coupled dynamics of quarks and gluons as well as on effects due to their collective motion on the physics of the early universe and in astrophysics.

**Keywords:** QCD in the early universe; phase transitions; hydrodynamical evolution; equation of state of super-dense matter; classical Yang-Mills fields; Dark Energy; Dark Matter; gluon condensate; effective Yang-Mills action; cosmic inflation

**PACS:** 98.80.Qc; 98.80.Jk; 98.80.Cq; 98.80.Es

#### **1. Introduction and Historical Perspective**

The strongly coupled dynamics of quarks and gluons has many important implications in particle physics, astrophysics, and cosmology [1–16]. The fundamental theory of strong interactions, known as Quantum Chromodynamics (QCD), provides a successful description of a variety of observables in high-energy hadronic collisions [17], hadronic masses [18], and, to a lesser extent, also of the properties of phases of the QCD matter [19,20]. While QCD is successful in the interpretation of short-distance phenomena (i.e., in the weakly coupled regime), a long-standing theoretical problem is a dynamical description of the color confinement phenomenon. The latter appears in the infrared (strongly coupled) regime of QCD and still remains the major unsolved problem of the Standard Model (SM) of particle physics [21,22].

Due to confinement, color-charged particles do not exist as free states at large spacetime separations. They are instead bound together into colorless collective excitations that evolve into a gas of hadrons. No exact dynamical transition in spacetime between the fundamental (parton) and the composite (hadron) states of QCD is known to date despite the wealth of phenomenological information available from particle and heavy-ion collision experiments. Therefore, one usually resorts to a heuristic description using the concept of quark-hadron duality [23,24] together with effective field theoretical (EFT) approaches [25]; this is used also in the framework of thermal field theory (for recent reviews of the latter, see also Refs. [26–28]). On the theory side, effective (typically, static or equilibrium) approaches, such as lattice QCD (LQCD) [19,21], are commonly being exploited while very little has been done on first-principle real-time evolution of QCD states [8].

The term "quark matter" was first used in 1970 by Itoh [29] in the context of neutron stars. Even before then, in 1965, Ivanenko and Kurdgelaidze [30] considered a star made of quarks. Since the mechanism for quark confinement was unknown at that time, they had to

**Citation:** Addazi, A.; Lundberg, T.; Marcianò, A.; Pasechnik, R.; Šumbera, M. Cosmology from Strong Interactions. *Universe* **2022**, *8*, 451. https://doi.org/10.3390/ universe8090451

Academic Editor: Antonino Del Popolo

Received: 27 April 2022 Accepted: 20 August 2022 Published: 29 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

assume that the quark masses are much larger than the masses of ordinary baryons. A few years later, however, following the realization that QCD exhibits asymptotic freedom [31,32], several authors have suggested that the transition from a hadronic phase to a one dominated by quarks and gluons may be relevant to describe the state of matter in the early universe or inside the neutron stars with a possibility to re-create such a condition also in the laboratory by colliding heavy ions [3,5,33–38].

The terms "hadronic plasma" [35] and "quark-gluon plasma" (or QGP) [36] were coined by Shuryak to describe a hypothetical state of matter existing at temperatures of order 100 MeV. The corresponding phenomena were expected to occur at a characteristic energy density close to 1 GeV/fm3. This makes a good analogy with a classical gaseous plasma in which electrically neutral gas at high enough temperatures turns into a statistical system of mobile charged particles [39]. While for such plasma the particle interactions obey the U(1)e*m* gauge symmetry of Quantum Electrodynamics (QED), in the former QCD case, the interactions between plasma constituents are driven by their SU(3)c color charges. For an exhaustive collection of key references tracing the development of theoretical ideas on the QGP up to 1990, see e.g., Ref. [2]. For a summary of later developments, see more recent reviews [8,10,40].

Let us note that contrary to initial oversimplified expectations [2], strongly interacting multi-particle systems feature numerous emergent phenomena that are difficult to predict from the underlying QCD theory, just like in condensed matter and atomic systems where the interactions are controlled by QED. In addition to the hot QGP phase, several additional phases of QCD matter were predicted to exist [15,41]. In particular, the long-range attraction between the quarks in the color anti-triplet (**3¯**) channel was predicted to lead to the color superconductivity (CSC) phase with condensation of <sup>1</sup>*S*<sup>0</sup> Cooper pairs [42,43]. This result was anticipated, though using a different reasoning, already in 1969 by Ivanenko and Kurdgelaidze [44], who predicted that the superconducting quark phase may be relevant for the super-dense star interiors. At high baryon density, an interesting symmetry breaking pattern SU(3)c× SU(3)L× SU(3)R× U(1)B → SU(3)c+L+R× Z(2) leading to the formation of quark Cooper pairs was found in QCD with three massless quark flavors (i.e., under an assumption that *mu* = *md* = *ms* = 0) [41,45]. This breaking of color and flavor symmetries down to the diagonal subgroup SU(3)c+L+<sup>R</sup> implies a simultaneous rotation of color and flavor called the color-flavor locking (CFL). It is expected that CSC and CFL phases might play important role in the equation of state (EoS) of neutron stars [46].

Another interesting phase of QCD matter called quarkyonic matter situated in the QCD phase diagram between the chirally restored and the confined phases was proposed in Ref. [47]. The quarkyonic matter is expected to exist at densities parametrically large compared to <sup>∼</sup>Λ<sup>4</sup> QCD when the number of colors *Nc* is large. Since gluons are in the adjoint representation of SU(3)*c*, their effects are scaled as <sup>∼</sup>*N*<sup>2</sup> *<sup>c</sup>* , and so, they dominate all quarkinduced ∼*Nc* effects. This provides the binding of gluons into quark-free states, so-called glueballs, and so, the quarkyonic matter has only *Nc* degrees of freedom (DoFs). This form of matter is expected to play some role in the structure of neutron stars [48]. The existence of another peculiar form of hadronic matter—the pion condensate—was suggested by Migdal already in the 1970s [49].

The rich phase structure of QCD at nonzero temperature and baryon chemical potential was recently reaffirmed by the proposed existence of phases with spatial modulations; see [50] and references therein. Their moat-shaped energy spectrum with a minimum of the energy over a sphere at nonzero momentum leads to a characteristic peak. In heavy-ion collisions at low energy, these new QCD phases are expected to leave their imprints in particle spectra and their correlations. Their cosmological implications are so far unexplored.

Our current knowledge of the QCD phase diagram is illustrated in Figure 1. Comparing this diagram to the phase diagram of water, see e.g., Ref. [8], one notices that (at least, theoretically) the complexity of the former approaches the latter.

**Figure 1.** The schematic phase diagram for QCD matter in terms of the temperature *T* and net baryon density *n* normalized to the cold nuclei baryon density *no*. From https://nica.jinr.ru/physics.php accessed on 27 April 2022 (see also Ref. [51]).

Experimental study of the QCD phase diagram at high temperatures, see Figure 2, dates back to the CERN SPS fixed-target program with the lead ion beams in 1995–2000 and covers the domain of the baryon chemical potential *μ*<sup>B</sup> =200 − 400 MeV [8]. With the advent of a first heavy-ion collider in 2000, the investigation of the *μ*<sup>B</sup> 0 region soon led to a discovery of the strongly interacting quark-gluon plasma (sQGP) at RHIC in 2005 [52–55]. The existence of this new phase of hot and strongly interacting QCD matter was five years later confirmed at order-of-magnitude higher energies of the LHC at CERN [8,40].

**Figure 2.** A schematic QCD phase diagram in the thermodynamic parameter space spanned by the temperature *T* and the baryonic chemical potential *μB*. The corresponding (center-of-mass) collision energy ranges for different accelerator facilities, especially the RHIC Beam Energy Scan (BES II) program, are indicated in the figure. Adapted from Ref. [56] (see also Ref. [57]).

Starting from 2010, it became possible to explore systematically the phase structure of hot and dense matter at nonzero baryon density and, in particular, to search systematically for the critical endpoint (CEP) of the QCD phase diagram. The CEP, a postulated secondorder phase transition point, is an expected endpoint of a line of the first-order phase transitions (FOPTs) that separates the low-temperature, low-density hadronic phase from a low-temperature, large-baryon number density QGP phase. Similarly to the watersteam transition where at the critical point, one finds bubbles of steam and drops of water intermixed at all length-scales from macroscopic, visible sizes down to atomic scales (with drops and bubbles near micron scale causing the strong light scattering called "critical opalescence" [58]), several interesting phenomena are also expected to occur near the CEP of the QCD phase diagram [57,59–62]. The search for the CEP is conducted by the STAR collaboration at RHIC within its Beam Energy Scan (BES) program at the energies indicated in Figure 1.

Current experimental and theoretical studies of the QCD phase diagram thus cover a wide region in temperature and baryon chemical potential (*T*, *μ*B), particularly, at small *μ*<sup>B</sup> 0 [19,20,63,64] and large *μ*<sup>B</sup> 100 − 600 MeV [20,41,57,60], see Figure 2. The red and black full circles denote the critical endpoints of the chiral and nuclear liquid-gas phase transitions, respectively. The (dashed) freeze-out curve indicates where the hadro-chemical equilibrium is attained at the final stage of the collision. The nuclear matter ground-state at *T* = 0 and *μ*<sup>B</sup> = 0.93 GeV and the approximate position of the QCD critical point at *μ*<sup>B</sup> ∼ 0.4 GeV are also indicated. The dashed line is the chiral pseudo-critical line associated with the crossover transition at low temperatures.

The hot and dense QCD matter is considered to be a dominant ingredient of the early universe evolution in its first few microseconds. Physics of heavy-ion collisions (HIC), therefore, provides necessary means for theoretical understanding of the cosmological processes at those time scales. In HIC theory, an important progress has been made when relativistic viscous fluid dynamics was formulated starting from the first principles in an EFT framework, which was based entirely on the knowledge of symmetries and longlived degrees of freedom, see e.g., Ref. [25] and Appendix B of this review. However, for proper understanding of the cosmological evolution, at least in a vicinity of the QCD phase transition epoch, the precise dynamical information on color-field media at finite temperatures is mandatory. Ongoing precision tests of QCD under extreme conditions, in particular those at the Large Hadron Collider (LHC) at CERN and the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory (BNL), are currently pushing the energy, temperature and density frontiers, opening up new largely unexplored possibilities for understanding also the cosmological QCD phase transition. There is strong hope that the growing amount of data and phenomenological concepts will eventually boost theoretical developments on infrared and finite-temperature dynamics of QCD. The latter is particularly relevant for understanding the real-time evolution of its ground state and the associated phase transitions as well as hadronization processes relevant for the dynamics of the early universe.

The QCD dimensional transmutation mechanism, breaking the conformal symmetry of the classical QCD action, has deep implications for the early universe evolution. Indeed, from higher to lower primordial plasma temperatures, QCD crosses a phase transition to a chiral symmetry-breaking ground state related to the color confinement phase. Thus, an attractive possibility is that the QCD vacuum energy may provide a source of universe acceleration and Dark Energy (DE) [65,66]. For the current status of this problem, see also Ref. [67] and references therein.

At high temperatures above the confinement scale ΛQCD, i.e., during the first microseconds after the Big Bang, the thermal bath in the early universe was dominated by the primordial QGP [11,12,68,69]. When the temperature of the universe decreases down to ΛQCD, the QGP dissolves out through collective hadronization phenomena. It is worth remarking that the QGP formation can be highly favoured under the very high-density conditions, where the matter chemical potential starts to be comparable to the QCD critical

scale. Indeed, this can happen in high-density objects such as in the core of neutron stars as an effect of the gravitational potential [9,70,71]. Nevertheless, the issue of whether the QGP exists inside the neutron stars is still highly controversial and is under intense debate in the literature. In fact, the critical phenomena of QCD dynamics and hadronization and their connections to the QCD ground-state evolution in real time have paramount consequences on the whole cosmological scene, which may also shed light on the late-time cosmic acceleration [65], the formation of Dark Matter (DM) [72], primordial black holes (PBHs) [73] and could be imprinted in the spectra of primordial density perturbations and gravitational waves (GWs) [74].

The main aim of our review is to provide a new critical sight on our current picture of quantum Yang-Mills (YM) field theories in the strong-coupling regime in a dynamical (i.e., non-stationary) spacetime background and in cosmology, in connection to the empirical knowledge that comes from particle physics measurements and cosmological data. In the following, unless otherwise noted, we will mainly exploit the standard natural units *h*¯ = *c* = *kB* = 1, where *kB* is Boltzmann constant, *c* is speed of light in the vacuum and *h*¯ = *h*/(2*π*), with *h* is the Planck constant.

The review is organized as follows. In Section 2, we discuss the nowadays standard scenario of the phase transitions in the early universe, making a connection to the production of primordial black holes and to super-dense weakly interacting saturated QCD matter. We also discuss possible applications of the axion dynamics to the early universe and close with the possible role of non-perturbative QCD ground-state cosmological evolution. For completeness, we also mention the possible role of the phase transitions in grand unified theories of particle interactions. In Section 3, we first introduce the basic notions of the hydrodynamical description of an expanding universe. There, we discuss simple models with constant speed of sound and then move on to a more complicated equation of states for the early universe. We also present current progress in the description of the dissipative effects in relativistic hydrodynamics. The section is finalized by an overview of the problematics regarding the Cosmological Constant and the Vacuum Catastrophe. Section 4 is devoted to a brief discussion of the real-time dynamics of the ground state in an effective action approach to quantum YM theories. We first discuss the YM ground state as time crystal; then, we develop an effective action approach providing the equation of state of the quantum ground state of the universe. The section is closed with a discussion of cosmological attractors—the solutions of the YM-Einstein equations using the Renormalization Group (RG) methods. Section 5 provides an overview of basic concepts of cosmic inflation models driven by YM dynamics in the early universe. Finally, a short summary is given in Section 6.

#### **2. The Phase Transitions in the Early Universe**

#### *2.1. The Phase Transitions in the Standard Model*

In the SM of elementary particle interactions, the dynamics of fireball expansion is based on the asymptotic freedom property of underlying non-Abelian gauge theories [31,32]. QCD is a quantum non-Abelian field theory, an important part of the SM, that describes the fundamental interactions between colored quarks and gluons. The generalization of classical electrodynamics to non-Abelian gauge theories was first studied and exemplified in SU(2) by Yang and Mills in 1954 [75]. The classical Lagrangian density of an SU(*N*) gauge theory reads,

$$\mathcal{L}\_{\rm cl} = -\frac{1}{4} F^{a}\_{\mu\nu} F^{a \; \mu\nu} \, , \tag{1}$$

in terms of the field strength tensor defined in terms of YM fields *A<sup>a</sup> <sup>ν</sup>* as *F<sup>a</sup> μν* = *∂μA<sup>a</sup> <sup>ν</sup>* <sup>−</sup> *∂νA<sup>a</sup> <sup>μ</sup>*+ *g*YM *f abcAb μA<sup>c</sup> <sup>ν</sup>*, where *f abc* are the structure constants of the SU(*N*) group. Throughout, *a*, *b*, ... denote internal indices of SU(*N*) in the adjoint representation. Here, the parameter *g*YM is known as the YM coupling constant. Gauge theories based on SU(*N*) are known as YM theories, and they became the target of a wider interest prompted by the discovery that massless particles may acquire a mass and a longitudinal polarization through spontaneous

symmetry breaking (or Higgs) mechanism of a massless YM theory [76–78]. The latter is a vital part of the SM framework realizing the classical Higgs mechanism of EW symmetry breaking that has found an excellent confirmation through the discovery of the Higgs boson [79,80].

In the framework of the quantum field theoretical approach, the YM field fluctuations are quantized around a given ground state through the first quantization procedure à la QED. However, due to self-interactions of the YM quanta, manifested via the term ∝*A*<sup>2</sup> in the field strength tensor, the quantum YM ground state acquires, in general, a very non-trivial structure. This structure is well understood in the weakly coupled (perturbative) regime of the theory, implying *g*YM 1, which is the case of the EW theory or in the UV regime of QCD with the so-called asymptotic freedom of color charges. The strongly coupled (non-perturbative) regime, in which *g*YM - 1 that is realized in particular in the infrared limit of QCD, corresponds to the color confinement phenomenon and has remained the subject of active research over the last few decades. In recent years, significant progress has been made in understanding of the quantum YM ground state in SU(*N*) gauge theories at finite temperatures, see e.g., Refs. [81,82].

As a result of a series of cosmological phase transitions that occurred in early universe during the first few microseconds after the Big Bang, new vacuum subsystems associated with breaking of the fundamental symmetries were formed. In the early universe, the SM predicts that cooling proceeds as a series of two phase transitions associated with the various spontaneous symmetry breakings of the corresponding gauge symmetries [12,68,83–86]. One at the temperature *T*EW *<sup>c</sup>* = 160 GeV [87] is responsible for spontaneous breaking of the EW symmetry providing masses to the elementary particles; see the left panel on Figure 3. Due to the large value of its critical temperature *Tc*, it is not amenable to experimental study under the laboratory conditions. The second and the only one accessible in the laboratory, QGP-to-hadronic matter phase transition happening at *T*QCD *<sup>c</sup>* ≈ 160 MeV [88], is related to the spontaneous breaking of the chiral symmetry and manifesting itself in the massless quark limit of the QCD Lagrangian. Since both phase transitions are considered to be analytic crossovers, the bulk motion of the corresponding plasmas did not depart from thermal equilibrium. Therefore, such transitions, if realized in nature, are not expected to generate cosmological relics [86,89,90] or to be helpful for a baryogenesis mechanism.

The QCD phase transition has occurred at characteristic temperatures of above 200 MeV that correspond to a cosmological time-scale of above 10−<sup>5</sup> s and the Hubble length-scale of approximately 10 km. The nature of the QCD phase transition is still a matter of intense debates in the literature [74,91–103], with results derived so far heavily relying either on lattice field theory methods applied to QCD [92,93], i.e., lattice QCD, or on holographic analyses of QCD at early cosmology [103]. For a thorough review on various aspects of the QCD transition epoch, see e.g., Refs. [12,104–108].

It is undeniable that an abrupt QCD phase transition occurring reversibly in the early universe would lead to a promising cosmological scenario, according to which a large part of the quark excess would be condensed into invisible quark nuggets—a possible explanation for DM only relying on QCD. As Witten suggested in Ref. [5], this would happen only if quark matter retains an energy per baryon which is less than 938 MeV: then, neutron stars might generate a quark matter component for cosmic rays, and detectable gravitational radiation could be produced during the QCD phase transition. Conversely, several recent studies drew a different conclusion, pointing toward the realization of second-order or crossover phase transition scenarios [101,102].

In a hot and dense QCD matter, the *u*, *d* and to some extent, depending on the temperature, also the *s* quarks become nearly massless, and the QCD Lagrangian acquires an approximate chiral symmetry SU(*NF*)L × SU(*NF*)R, with the number of massless quark flavors *NF* = 2 (*u*, *d*) or 3 (*u*, *d*,*s*). At low *T* < *T*QCD *<sup>c</sup>* , the QCD vacuum becomes unstable, and this symmetry is spontaneously broken by *q*−*q*¯ pairing. The corresponding order parameter *qq*¯ =−(245 MeV)3, known as the chiral quark condensate, gives rise to masses of light hadrons as well as to constituent masses of *u*, *d* quarks and to some extent also to

the *s* quark; see the left panel of Figure 3. Recent lattice calculations with *mu* = *md* = 0 and strange quark having its physical mass reveal that the chiral phase transition occurs at *T*QCD *<sup>χ</sup>* =132+<sup>3</sup> <sup>−</sup><sup>6</sup> MeV [109].

**Figure 3.** (**Left**) Quark masses in the QCD vacuum and the Higgs vacuum. A large fraction of the light quark (*u*, *d*,*s*) mass is due the chiral (*χc*) symmetry breaking in the QCD vacuum, with numerical values from Ref. [110] (see also Ref. [111]). (**Right**) The effective number of relativistic DoFs *g*eff in the cosmological plasma in the SM as a function of temperature *T*, taking into account interactions between particles, obtained with both perturbative and lattice methods. From Ref. [112].

Let us note that in addition to the above scenario when the thermal history of the universe proceeds by a sequence of phase transitions from a more symmetric to a less symmetric state of matter, there is also a possibility of the reverse evolution when part of the zero-temperature unbroken gauge group of the SM or other gauge theory might have been broken in the early universe by thermal effects. As first noted by Weinberg [113] in the context of an *O*(*n*) × *O*(*n*) gauge theory, with decreasing *T*, one may encounter a transition to a state of higher symmetry *O*(*n*) × *O*(*n* − 1) → *O*(*n*) × *O*(*n*). Within the minimal extensions of the SM containing an additional color triplet scalar field, the scenario in which the early universe underwent an epoch when SU(3)*<sup>c</sup>* was spontaneously broken but later restored was analyzed in Ref. [114]. The attractiveness of such a multi-step phase transition scenario stems from the fact that it may generate the observed baryon asymmetry of the universe [115].

To describe the evolution of energy density (*T*) and entropy density *s*(*T*) of the early universe, it is customary to normalize both quantities to their values 0(*T*) and *s*0(*T*) corresponding to an ideal massless Bose gas with a single degree of freedom (DoF) [112,116]

$$g\_{\rm eff}(T) \equiv \frac{\epsilon(T)}{\epsilon\_0(T)}, \quad \epsilon\_0(T) = \frac{\pi^2}{30} T^4 \,. \tag{2}$$

$$h\_{\rm eff}(T) \equiv \frac{\mathbf{s}(T)}{\mathbf{s}\_0(T)}, \quad \mathbf{s}\_0(T) = \frac{2\pi^2}{45} T^3,\tag{3}$$

and call *g*eff(*T*) and *h*eff(*T*) the effective numbers of DoFs in energy and entropy, respectively. For the particular case of a non-interacting gas consisting of *N*<sup>F</sup> Dirac fermions, *N*<sup>V</sup> massive vectors, *N*V0 massless vectors and *N*<sup>S</sup> neutral scalars, the two functions are identical and read

$$\mathcal{g}\_{\rm eff}(T) = h\_{\rm eff}(T) = \frac{7}{8} 4\text{N}\_{\rm F} + 3\text{N}\_{\rm V} + 2\text{N}\_{\rm V0} + \text{N}\_{\rm S} \tag{4}$$

where the prefactors account for the DoF of each of the considered particles.

It is worth mentioning that in a generic case of interacting (non-ideal) gas *g*eff(*T*) and *h*eff(*T*), they are not identical and depend on temperature. Using the relationship

*p* = *sT* − between the pressure *p*, the entropy density *s*, the energy density and the generalized EoS parameter *w*,

$$
\boldsymbol{\nu} = \boldsymbol{w} \boldsymbol{\epsilon} \,. \tag{5}
$$

and the speed of sound *cs* can be expressed in terms of the effective DoF measures *g*eff(*T*) and *h*eff(*T*),

$$w(T) = \frac{sT}{\epsilon} - 1 = \frac{4h\_{\rm eff}(T)}{\Im g\_{\rm eff}(T)} - 1,\tag{6}$$

$$c\_s^2(T) = \frac{dp}{d\epsilon} = T\frac{ds}{d\epsilon} + s\frac{dT}{d\epsilon} - 1 = \frac{4}{3} \left[ \frac{4h\_{\rm eff}(T) + Th\_{\rm eff}^{'}(T)}{4g\_{\rm eff}(T) + Tg\_{\rm eff}^{'}(T)} \right] - 1,\tag{7}$$

where the prime indicates differentiation with respect to temperature *T* [116]. It is worth mentioning that the causality condition between the speed of sound and the speed of light *cs* ≤ 1 induces the inequality

$$\frac{h\_{\rm eff}(T)}{g\_{\rm eff}(T)} \le \frac{3}{2} \,\,\,\,\,\tag{8}$$

with an upper bound saturated for *w* = 1 corresponding to the case of absolutely stiff fluid.

Taking into account all permissible interactions in the SM, one can calculate either directly [117] or on a lattice the temperature dependence of (*T*) and *s*(*T*) and extract corresponding DoFs. This is illustrated on the right panel of Figure 3 showing the temperature dependence of *g*eff(*T*) in the SM. For realistic values of both *g*eff(*T*) and *h*eff(*T*) for a wide temperature interval from 10 keV to 10 TeV, see Ref. [117].

The simultaneous presence of the EW and the QCD matter in thermal equilibrium is one of the remarkable differences between the QGP produced in accelerator experiments and the primordial QGP in the early universe [13,118,119]. To find out which form of matter prevails, let us use Equation (2) and compare the number of DoFs in an ideal massless gas consisting only of EW or QCD matter. Including only the particles which at the temperature *T T*EW can be considered as massless, we obtain *g*EW eff <sup>=</sup> <sup>7</sup> <sup>8</sup> (12 + 6) + 2 = 17.75 DoFs for the EW case. The first term in the brackets corresponds to charged leptons *e*, *μ*, *τ*, and the second term corresponds to neutrinos *νe*, *νμ*, *ντ*. The last term corresponds to photons. For non-interacting or weakly interacting QCD matter, Equation (4) reduces to

$$\mathbf{g}\_{\rm eff}^{\rm QCD} = 2 \times 8 + \frac{7}{8} (3 \times N\_F \times 2 \times 2) \,\,,\tag{9}$$

where the first term accounts for the two spin and *N*<sup>2</sup> *<sup>c</sup>* − 1 = 8 color DoFs of the gluons and the second term accounts for *Nc* = 3 colors, *NF* flavors, two spin and two particleantiparticle DoFs of the quarks. Including only the quarks with the mass *mi*/*T* 0, *i* = (*u*, *d*), *s*, *c*, *b*, we obtain successively *g*QCD eff = (37, 47.5, 56, 68.5) DoFs. In thermal equilibrium at the temperature *T T*EW and for *NF* active quark flavors, the QGP contains a factor of *g*QCD eff /*g*EW eff 2 − 4 more energy and pressure than those for the EW matter. For temperatures *<sup>T</sup> <sup>T</sup>*EW *<sup>c</sup>* , deep inside the EW era, all six quarks *u*, *d*,*s*, *c*, *b*, *t* can be considered to be massless, cf. the left panel of Figure 3, and *g*QCD eff = 79. At the same time, the EW matter acquires *g*EW eff <sup>=</sup> <sup>7</sup> <sup>8</sup> (12 + 6) + 8 + 4 = 26.75 DoFs, where 8 = 2 × (3 + 1) are the DoFs of massless gauge bosons, *W*±, *W*0, *B*0, and the last term is due to the Higgs scalar doublet. For this case, the QGP has a factor of *g*QCD eff /*g*EW eff 3 larger energy density and pressure than those of the EW matter. Hence, we conclude that the QGP was the most dense form of matter filling the early universe during both the QCD and EW epochs.

#### *2.2. Creation of Primordial Black Holes during the Phase Transitions*

According to inflationary theories, initially, very small inhomogeneities in the matter distribution were produced by the end of the exponential expansion regime. Such inhomogeneities filling the early universe are described by the metric perturbations *δgμν* which can be decomposed into three irreducible pieces—scalar, vector and tensor ones, see, e.g., Refs. [84,120]. While the scalar part is induced by energy density fluctuations *δ*, the vector and tensor perturbations are related to the rotational motion of the fluid and to the gravitational waves, respectively [84]. Given the scope of this review, in the following, we focus only on one spectacular phenomenon related to metrics fluctuations in the early universe—the matter collapse into primordial black holes (PBHs).

Whereas their existence was proposed already a half-century ago first by Zeldovich and Novikov [121] and later by Hawking [122], it was the detection of gravitational waves from mergers of tens of solar mass *M* black hole binaries [123] which has led to a surge of current interest in the PBHs as a Cold Dark Matter (CDM) candidate [73,124–126]. It can be shown that the creation of PBHs due to the gravitational collapse of hot and dense matter occurs for the density contrast *δ* = *δ*/ exceeding the critical threshold *δc*(*w*[*T*]) ≈ 0.3 − 0.45, which generally depends on the EoS parameter *w* [73,125]. Such large values of *δ* can be generated, e.g., during a period of inflation in the very early universe [126] or during an intermediate period dominated by long-lived massive particles (for recent work, see e.g., Ref. [127] and references therein) or when the universe in the course of the phase transition passes a local minimum in the pressure-to-energy density ratio *w* = *p*/ [125].

For the PBHs forming from Gaussian inhomogeneities with root-mean-square amplitude *δ*rms, the present CDM fraction for PBHs with a mass around *M* is found as [125,128]

$$f\_{\rm PBH}(M) \approx 2.4\beta(M)\sqrt{\frac{M\_{\rm eq}}{M}}, \quad \beta(M) = \frac{2}{\pi} \int\_{x}^{\infty} e^{-y^2} dy, \quad x = \frac{\delta\_c(w[T(M)])}{\sqrt{2}\delta\_{\rm rms}(M)}, \tag{10}$$

where *β*(*M*) is the fraction of horizon patches undergoing collapse to PBHs when the temperature of the universe is *<sup>T</sup>*, *<sup>M</sup>*eq <sup>=</sup> 2.8 <sup>×</sup> <sup>10</sup>17*<sup>M</sup>* is the horizon mass at matter-radiation equality, and the numerical factor comes from the ratio of measured baryon Ω<sup>b</sup> and CDM ΩCDM abundances. In Equation (10), we have explicitly taken into account dependence of the critical overdensity *δ<sup>c</sup>* on the EoS parameter *w*(*T*). The temperature depends on the PBH mass *<sup>M</sup>* as *<sup>T</sup>* <sup>≈</sup> <sup>200</sup>√*<sup>M</sup>* /*<sup>M</sup>* MeV. Note that the parameter *<sup>δ</sup>*rms(*M*) appearing in Equation (10) can be always adjusted to counterbalance the theoretical uncertainties in the value of *δ<sup>c</sup>* so that the current PBH DM fraction is preserved [128]. It is worth mentioning that in the scenarios where PBHs are formed during inflation, their abundance is larger than the Gaussian result by orders of magnitude, but also the mass function has a more pronounced tail at larger masses [126].

In fact, there are a plethora of other mechanisms for PBHs formation (including besides the already mentioned FOPTs, bubble collisions, and the collapse of cosmic strings, necklaces, domain walls, non-standard vacua, etc., see e.g., the recent reviews [73,125]). In the following, in conformity with the topic of our review, we will concentrate on the softest point (SP) mechanism of creation of the PBHs discussed in Ref. [128]. Its virtue stems from the fact that by tracing the origin of PBHs to the SM phase transitions, it is capable of explaining the provenance of part, if not all, of the CDM in the universe [124].

The SP, corresponding to a local minimum in the pressure-to-energy density ratio *w* = *p*/ as given in Equation (6), gives rise to elongation of the expansion time of the hot and dense matter. In HICs, the interest in locating the SP was started by the recognition that the longest-lived fireball might provide a clear signature of the QGP-to-hadron phase transition [129]. Shortly after that, the formation of horizon-size PBHs due to a substantial reduction of pressure during adiabatic collapse in the course of the QCD transition was analyzed in the context of the early universe in Refs. [130,131]. Even though the previously used assumption of the first-order character of the phase transition was later on replaced by a crossover scenario, the lattice calculations have found a local minimum in *w* = 0.145(4) at *T* = 159(5) MeV [132].

To become acquainted with the influence of the SPs on the cooling of the universe during its radiation-dominated era, let us follow Ref. [128] and inspect the behavior of the function *g*eff(*T*) shown on the right panel of Figure 3. Let us focus on the temperatures when a part of the radiation matter transforms into a non-relativistic matter. Starting from *T* ≈ 200 GeV downwards, this happens first to the top quark at *T* ≈ *m*<sup>t</sup> = 172 GeV, which is followed by the Higgs boson at 125 GeV and the *Z* and *W* bosons at 92 and 81 GeV, respectively. The fact that these particles become non-relativistic at nearly the same time of universe expansion induces a significant drop in the number of relativistic DoFs, from *g*eff = 106.75 down to *g*eff = 86.75. Further changes at the *b*,*c*-quark and *τ*-lepton thresholds are too small to be noticeable. Hence, further on, *g*eff remains approximately constant until the QCD transition at around 160 MeV. The number of relativistic DoFs then falls abruptly to *g* = 17.25. A little later, pions become non-relativistic, and then muons, yielding *g*eff = 10.75. Thereafter, *g*eff remains constant until *e*+*e*<sup>−</sup> annihilation and neutrino decoupling at around 1 MeV, when it drops down to *g*eff = 3.36 [128].

Provided that total entropy is conserved, an abrupt reduction of *g*eff(*T*) leads to a sudden drop in the speed of sound *cs*(*T*), cf. Equation (7), and hence to a drop of pressure, *p* = *w*(*T*), cf. Equation (6). The effect is clearly visible on the left panel of Figure 4 showing the four periods in thermal history of the universe when *w*(*T*) reaches its local minimum. After each period, *w* returns back to its relativistic value of 1/3, but each sudden drop modifies the probability of gravitational collapse of any large curvature fluctuations present at that time [128].

**Figure 4.** (**Left**) EoS parameter *w* as a function of temperature *T*. The gray vertical lines correspond to the masses of the electron, pion, proton/neutron, *W*, *Z* bosons and top quark, respectively. The gray dashed horizontal line indicates value of *w* = 1/3. Adapted from Ref. [128]. (**Right**) The mass spectrum of PBHs *f*PBH(*M*) in solar mass units *M* . The gray vertical lines correspond to the EW and QCD phase transitions and *e*+*e*<sup>−</sup> annihilation. The vertical colored lines indicate the masses of the three LIGO-Virgo events. Gray curves are constraints from microlensing (M), ultra-faint dwarf galaxies and Eridanus II (E), X-ray/radio counts (X), and halo wide binaries (W). The accretion constraint (A) is shown dashed, as it relies on uncertain astrophysical assumptions. Adapted from Ref. [128].

Consider one cooling period *T*<sup>1</sup> < *T* < *T*<sup>2</sup> with *w*(*T*) < *w*(*T*1,2) = 1/3 centered around the local minimum *w*(*T*SP) and define the quantity,

$$
\Delta h\_{\rm eff}(T) \equiv \mathcal{g}\_{\rm eff}(T) - h\_{\rm eff}(T) \,, \tag{11}
$$

measuring departure from the *w* = 1/3 case; see Equation (6). At the endpoints Δ*h*eff(*T*2) = Δ*h*eff(*T*1) = 0 but for *w*(*T*) < 1/3, it is always positive Δ*h*eff(*T*) > 0, cf. Equation (7). Hence, the initial drop in the entropy DoFs, *h*eff(*T*), always precedes the jump in the energy density DoFs, *g*eff(*T*). This leads to the following "coarse-grained" scenario for the PBH formation: the reduction in *h*eff(*T*) occurring for *T*<sup>2</sup> > *T* > *T*SP is followed by a fall in *g*eff(*T*) for *T*<sup>1</sup> > *T* > *T*SP. An excess in entropy ∼ Δ*h*eff(*T*SP) lost by the radiation during its cooling is dumped into the collapsing matter, emerging eventually in the form of PBHs—the matter with the largest entropy density in the universe [133]. This may explain why even at the present stage of the universe evolution, there is by a huge factor far more entropy in

supermassive black holes (BHs) in galactic centers than in all other sources of entropy put together [134].

Assuming that the amplitude of the primordial curvature fluctuations is approximately scale-invariant [128], one obtains from Equation (10) the mass spectrum of PBHs *f*PBH(*M*) shown on the right panel of Figure 4. The peaks at *<sup>M</sup>* <sup>10</sup>−6, 2, 30 and 106*<sup>M</sup>* correspond to the EW and QCD phase transitions, to pions and muons becoming non-relativistic and to *e*+*e*<sup>−</sup> annihilation, respectively. The latter may also provide seeds for the supermassive BHs' formation in galactic nuclei. The largest contribution to *f*PBH(*M*) comes from the PBHs formed at the QCD transition epoch and that would naturally have the Chandrasekhar mass (1.4 *M* ) [128]. Moreover, the peak in the range 1–10 *Mrm* could explain the LIGO/Virgo observations [123]. The latter favor mergers with low effective spins as expected for PBHs, but it is hard to explain BHs of stellar origin [135].

The simple analytical models that describe the dynamical process of gravitational collapse which may be relevant for PBH formation were studied in Ref. [136]. It is also worth noting that the gravitational collapse of large inhomogeneities during the quarkhadron transition epoch may also explain the baryon asymmetry of the universe [137]. The asymmetry can be generated in local hot spots through the violent process of PBH formation at the QCD transition triggered by a sudden drop in the radiation pressure and the presence of large amplitude curvature fluctuations caused by the axion field—the subject to be discussed in Section 2.6.

#### *2.3. Perturbative and Strongly Coupled Regimes of QCD*

An important contribution to the effective number of relativistic DoFs, *g*eff, comes from the hadron-to-QGP phase transition—see a big jump in the interval 102 *T* 103 MeV in Figure 3 (right panel). At higher temperatures, in the QGP region, the strength of the interactions between the quarks and gluons is set by the QCD coupling *α*S(*Q*) which at the one-loop order of perturbation theory takes the form,

$$a\_{\rm S}(Q) \simeq \frac{2\pi}{b\_0 \ln(Q/\Lambda\_{\rm QCD})}, \qquad b\_0 = 11 - \frac{2}{3} N\_{\rm F} \tag{12}$$

where *Q* is the momentum transferred during the interaction, ΛQCD 200 MeV is the characteristic QCD scale, and *NF* is the number of active quark flavors. The logarithmic decrease of *α*S(*Q*) with increasing *Q*, i.e., with decreasing distance among the quarks and gluons, is due to the fact that, in contrast to the photon in QED, the force carriers in QCD, the gluons, have color charge. Their exchanges in higher-order processes involving both the quarks and the gluons occur more frequently with increasing *Q* and lead to a color charge spread (or anti-screening). Indeed, the gluon multiplicity increases at low momentum fractions corresponding to the limit of large energies. Dilution of the initial color charge is responsible for the weakening of *α*<sup>S</sup> at small distances - <sup>Λ</sup>−<sup>1</sup> QCD, i.e., when the quark experiences a large momentum transfer *Q*, see Equation (12). This effect known as asymptotic freedom [31,32,138] is illustrated on the left panel of Figure 5 where the values of the *α*S(*Q*) extracted from proton-(anti)proton and lepton-proton collisions are shown [139]. In agreement with Equation (12), a slow logarithmic decrease from *α*S(*Q*min =5 GeV) = 0.22 to *α*S(*Q*max =1500 GeV) = 0.08 is observed.

Before proceeding further, let us recall that quantum field theory (QFT) at finite temperature *T* is often considered to be equivalent to Euclidean QFT in a space which is periodic, with period 1/*T* along the "imaginary time" axis (for a recent review of this subject, see e.g., Refs. [26,140] and references therein). Thus, in order to formulate the theory at *T* > 0 using its variant at *T* = 0, one should replace zero components of all 4-momenta *k<sup>μ</sup>* in the Euclidean integrals by the discrete Matsubara frequencies—2*πnT* for bosons and (2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)*π<sup>T</sup>* for fermions, and sum over *<sup>n</sup>* <sup>∈</sup> <sup>Z</sup> instead of integrating over *<sup>k</sup>μ*. Consequently, the average momentum transferred during the interactions in the hot medium *Q* can be related to the temperature as *Q* = 2*πT*. In particular, the maximum value of the momentum

transferred *Q*max =1500 GeV, which has been so far measured in *pp* collisions at the LHC, roughly corresponds to the "temperature" *<sup>T</sup>*max <sup>=</sup> *<sup>Q</sup>*max/(2*π*) 240 GeV *<sup>T</sup>*EW *<sup>c</sup>* .

**Figure 5.** (**Left**) The QCD coupling *α*S(*Q*) as a function of the momentum transfer scale *Q* obtained by using the MSTW2008 NLO PDF set [141]. Adapted from Ref. [139]. (**Right**) The partonic phase diagram showing evolution of the partons' density and size as a function of their rapidity *Y* = ln(1/*x*) and the logarithm of momentum transfer squared ln *Q*2. Adapted from Ref. [142].

The asymptotic freedom formula given by Equation (12) is based on the applicability of QCD perturbation theory to the processes at high momentum transfers, as it is a well-known fact that the perturbative expansion is an example of asymptotic series. It looses its validity with decreasing *Q* when the perturbative approximation breaks down. Interestingly, by performing the matching of the fundamental theory onto the effective chiral Lagrangian at the infrared scale *Q* 4*π f<sup>π</sup>* 1 GeV, at which the ranges of validity of perturbative QCD and chiral perturbation theory (describing interactions among low-momentum hadrons) descriptions meet, one can infer the information about the behavior of *α*S(*Q*) at large distances. Such a matching implies that the QCD coupling in the infrared region is "frozen" at *α*SIR 0.56 [143], incidentally at twice the upper scale value of *α*S(*Q*) shown on the left panel of Figure 5.

Let us note that evolution of the strong coupling parameter *α*S(*Q*) described at the leading order by Equation (12) is valid only for DoFs that dominate the thermodynamical evolution, i.e., for the partons (quarks and gluons) with momenta of order *T* = *Q*/2*π*, and it does not apply to the long wavelength non-perturbative modes residing at the length scales of - > *T*−1. Those modes are occupied by a liquid in which neighboring "unit cells" are tightly coupled to each other [144]. This strongly coupled regime [10] makes the QGP behave as the ideal fluid [145,146]. The fluidity of the QGP was first established in the collisions of ultra-relativistic nuclei at RHIC [53–55] and later confirmed at higher energies of the LHC [147–149]. The most prominent signals of the strong interaction in the deconfined bulk manifest in a collective flow of matter [8,146] and in a spectacular phenomenon of suppression of very energetic partons passing through the QGP medium [8,150,151]. Direct evidence for the non-perturbative character of deconfined matter comes from the low-momentum spectra of direct photons measured in Au+Au collisions at RHIC. The temperatures obtained from the spectra *T* 220 MeV [152] point to the initial temperatures *T*ini 300–600 MeV at early times of *τ*<sup>0</sup> = 0.6–0.15 fm/*c*, which are way below the perturbative regime of QCD.

By definition, plasma is a state of matter in which charged particles interact via long-range (massless) gauge fields [153]. This distinguishes it from neutral gases, liquids, or solids in which the inter-particle interaction is of short range. So, plasmas themselves can be gases, liquids, or solids, depending on the value of the plasma parameter Γ, which

is the ratio of interaction energy to kinetic energy of the particles forming the plasma [39]. For a classical plasma of *N* particles with charge *Ze* occupying a volume *V*,

$$
\Gamma \equiv \frac{(Z\varepsilon)^2}{ak\_B T}, \qquad a(T) = \left(\frac{3V}{4\pi N}\right)^{1/3} \approx 0.62 n(T)^{-1/3} \,\text{s}.\tag{13}
$$

where *n*(*T*) = *N*/*V* is the temperature-dependent particle number density. While most plasmas are ideal with Γ < 10−3, a strongly interacting plasma has Γ - 1. For plasma with *<sup>Z</sup>* 1 at the temperature *<sup>T</sup>* 106 <sup>K</sup> 100 eV, the number density *<sup>n</sup>* must be as high as <sup>10</sup><sup>26</sup> cm−<sup>3</sup> to make <sup>Γ</sup> 1 [39]. Ion plasma in a white dwarf has <sup>Γ</sup> = 10–200, in transient plasmas produced in explosive shock tubes, the values of Γ = 1–5 are found [39]. A more down-to-earth example is table salt (NaCl), which can be considered as a crystalline plasma made of permanently charged ions Na<sup>+</sup> and Cl<sup>−</sup> [153]. At temperatures of *<sup>T</sup>* <sup>≈</sup> 103 K, still too small to ionize non-valence electrons, the table salt transforms into a molten salt, which is a liquid plasma with Γ ≈ 60.

Generalization of Equation (13) to the QGP case was suggested in Ref. [154]

$$
\Gamma \simeq 2 \frac{\mathcal{C}\_{\mathfrak{q}, \mathfrak{g}} a\_{\mathbb{S}}}{aT}, \qquad \mathcal{C}\_{\mathfrak{q}} = \frac{N\_{\mathfrak{c}}^2 - 1}{2N\_{\mathfrak{c}}} = \frac{4}{3}, \qquad \mathcal{C}\_{\mathfrak{g}} = N\_{\mathfrak{c}} = 3,\tag{14}
$$

where the strong coupling *α*<sup>S</sup> is a slowly varying function of temperature, *Cq* and *Cg* are the Casimir invariants of fundamental and adjoint irreducible representations of the color SU(3)*<sup>c</sup>* group corresponding to quarks and gluons, respectively, and *a* = *a*(*T*) is the average inter-parton distance at a given temperature *T* as follows from Equation (13). The factor 2 in Equation (14) takes into account the equal importance of chromoelectric (CE) and chromomagnetic (CM) interactions in ultra-relativistic systems. For ideal massless QCD gas with *NF* active quarks and *dF* = *g*QCD eff degrees of freedom, see Equation (9), the particle number density reads

$$m = d\_F \frac{\zeta(3)}{\pi^2} T^3 \approx d\_F \left(\frac{T}{2}\right)^3, \qquad d\_F = 2 \times 8 + \frac{7}{8} (3 \times N\_F \times 2 \times 2) \,. \tag{15}$$

From Equations (13) and (15), it follows that *<sup>a</sup>* 1.24*d*−1/3 *<sup>F</sup> <sup>T</sup>*−1, and so, the term *aT*, appearing in the denominator of Equation (14), depends on *T* through the temperaturedependent number of active quark flavours *NF*(*T*) only. Consequently, for an ideal massless QCD gas, the temperature-dependence of the plasma coupling parameter Γ is driven by *α*S*d*1/3 *<sup>F</sup>* . For the QGP created in HICs at RHIC *T* ≈ 200 MeV and *α*<sup>S</sup> = 0.3–0.5 with *NF* = 2, Equations (15) and (14) yield *dF* = 37 and Γ 2–8 well inside the strongly coupled regime. At much higher temperatures, say, at *T T*EW, with *α*<sup>S</sup> = 0.08 and *NF* = 5, we obtain *dF* = 52.5 and Γ 0.5–1.5—the value located in the vicinity of the strongly coupled regime. At even higher temperatures, the number of active quark DoFs saturates at *NF* = 6 and the evolution of Γ(*T*) becomes solely driven by the (logarithmically decreasing) QCD coupling *α*S(*T*), cf. (12). Let us note that the ideal gas approximation serves only as a lower estimate of Γ because it ignores the interactions in the partonic liquid. The latter will slow down the temperature dependence of the average inter-parton distance *a*, thus weakening the strong coupling parameter dependence on *T*.

A more in-depth approach to strongly coupled non-Abelian plasmas [155] was expected to come from the gauge/string duality [156]—a correspondence between *d*-dimensional conformal QFT and (*d* + 1)-dimensional string or gravity theory. In these theories, the graviton needs not live in the same spacetime as the QFT, but due to the holographic principle, the description of gravity within a volume of spacetime can be thought of as encoded on a lower-dimensional boundary to the region in the formalism of conformal field theory [157,158]. However, the inherently conformal character of the gauge QFT used in the duality with antide Sitter gravity (AdS) is at variance with QCD where the scale invariance is broken by the confinement scale (for a recent review of confinement dynamics, see Ref. [22]), causing

the running of the coupling. This limits the applicability of gauge/string duality [156] to temperatures *<sup>T</sup> <sup>T</sup>*QCD *<sup>c</sup>* and hence to weak couplings.

#### *2.4. QCD at High Parton Densities and Saturation*

To proceed further, a different type of analysis of the QCD dynamics of partonic matter is needed. There are two independent paths along which the density of partons can evolve, and these are illustrated in the right panel of Figure 5. The two together form the basis of our current understanding of high-energy scattering in QCD.

The first path follows the development of partonic cascade in variable *Q*. For partons that occupy a transverse area 1/*Q*2, the increase of *Q* and hence of the temperature *T* ∼ *Q* leads to dilution of their density. The process is controlled by the Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) [159–161] equations describing the evolution of partonic density as a function of evolution variable ln(*Q*2/Λ<sup>2</sup> QCD) [17,162].

The second path follows the development of a parton shower in variable *x* = *k*+/*P*+, which is a fraction of the light cone momentum<sup>1</sup> *P*<sup>+</sup> of the parent parton, which has radiated a parton emerging with the light-cone momentum *<sup>k</sup>*+. In the *<sup>x</sup><sup>⊥</sup>* plane transverse to the direction of the fast-moving primary parton, the partonic cascade initiated by the primary parton can be visualized as a Brownian motion-like trajectory developing from *x* = 1 toward *x* → 0. The corresponding Gribov diffusion process is controlled by the so-called evolution parameter *Y* = ln(1/*x*), leading to a difference in rapidity between the primary and radiated partons, with a coefficient being the diffusion constant proportional to *α*S. Its evolution in the *Y* variable is described by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation that is complementary to the DGLAP evolution realized in ln *Q*<sup>2</sup> (see standard textbooks, e.g., Refs. [17,162] and references therein).

At fixed *Q*, the radiated partons (mostly, soft gluons with *x* 1) are typically of the same size. When a parton-parton interaction cross-section <sup>∼</sup>*α*S/*Q*<sup>2</sup> multiplied by *xGA*(*x*, *Q*2)—the probability to find at fixed *Q* a parton carrying a fraction *x* of the parent parton momentum—becomes comparable to the geometrical cross section *πR*<sup>2</sup> *<sup>A</sup>* of the object *A* occupied by the gluons, the partons "overlap". The repulsive interactions among the gluons ensure that their occupation number *fg* (the number of gluons with a given *x* multiplied by the area each gluon fills up divided by the transverse size of the object) saturates at *fg* ∼ 1/*α*S. Note that this is a very generic behavior—the same density scaling as the inverse interaction strength *α*−<sup>1</sup> is characteristic of a number of condensation phenomena such as the Higgs condensate, see, e.g., Ref. [163], or superconductivity [164].

The phenomenon of saturation [165] is thus important for gluons with transverse momenta *k*<sup>⊥</sup> ≤ *Qs* [166,167], where

$$Q\_s^2(\mathbf{x}) = \frac{a\_\mathcal{S}(Q\_s)}{2(N\_c^2 - 1)} \frac{\mathbf{x} G\_A(\mathbf{x}, Q\_s^2)}{\pi R\_A^2} \sim \frac{1}{\mathbf{x}^\lambda} \tag{16}$$

is the *x*-dependent saturation scale representing a fixed point of the parton density evolution in *x* or, equivalently, the emergent "close packing" scale [167]—see the right panel of Figure 5 where the saturation line ln *Q*<sup>2</sup> *<sup>s</sup>*(*Y*) = *λY*, *Y* = ln *x* is also displayed. Such gluons form a highly coherent configuration called Color Glass Condensate (CGC) [142,167,168], or glasma [169], which due to the high occupation number *fg* has properties of QCD in the classical regime [166].

At high temperatures, one usually expects that quantum effects become less important. To show that, for the CGC, we follow Ref. [166] and write the gluonic part of the QCD action in terms of the gauge field potential <sup>A</sup>*<sup>a</sup> <sup>μ</sup>* and field strength <sup>F</sup>*<sup>a</sup> μν* which are obtained

by rescaling of their equivalents *A<sup>a</sup> <sup>μ</sup>* and *F<sup>a</sup> μν* used in a more traditional approach when the coupling constant *gs* multiplies the interaction terms in the Lagrangian,

$$A^a\_{\mu} \to \mathcal{A}^a\_{\mu} \equiv \mathcal{g}\_s A^a\_{\mu}, \quad F^a\_{\mu\nu} \to \mathcal{g}\_s F^a\_{\mu\nu} \equiv \mathcal{F}^a\_{\mu\nu} = \partial\_{\mu} \mathcal{A}^a\_{\nu} - \partial\_{\nu} \mathcal{A}^a\_{\mu} + f^{abc} \mathcal{A}^b\_{\mu} \mathcal{A}^c\_{\nu}, \tag{17}$$
 
$$\square \quad \square \quad f\_{\square \quad \square \quad \square \quad \square}, \quad \square \quad f\_{\square \quad \square \quad a \square \quad \square}, \quad \square$$

$$S\_{\mathcal{S}} = -\frac{1}{4} \int F^{a}\_{\mu\nu} F^{\mu\nu,a} d^4 \mathbf{x} = -\frac{1}{4\mathcal{S}\_s^2} \int \mathcal{F}^{\mu\nu,a} \mathcal{F}^a\_{\mu\nu} d^4 \mathbf{x} \,, \tag{18}$$

where *<sup>f</sup> abc* with (*a*, *<sup>b</sup>*, *<sup>c</sup>*) ∈ {1, ... , 8} are the SU(3) group structure constants. For a classical configuration of gluon fields, by definition, <sup>F</sup>*<sup>a</sup> μν* does not depend on the coupling, and the action is large, *Sg h*¯. The number of quanta in such a configuration is then

$$f\_{\mathcal{S}} \sim \frac{S\_{\mathcal{S}}}{\hbar} \sim \frac{1}{\hbar \mathcal{g}\_s^2} \rho\_4 V\_{4\ \prime} \tag{19}$$

where we rewrote Equation (18) as a product of four-dimensional gluon condensate density *<sup>ρ</sup>*<sup>4</sup> ∼ F*μν*,*a*F*<sup>a</sup> μν* and spacetime volume *V*4. The number of quanta *fg* in such a configuration depends only on the product of the Planck constant *h*¯ and the strong coupling squared *g*<sup>2</sup> *<sup>s</sup>* = 4*πα*S. The classical limit *h*¯ → 0 is indistinguishable from the weak coupling limit *g*<sup>2</sup> *<sup>s</sup>* → 0 [166]. Thus, the weak coupling limit of small *α*<sup>S</sup> corresponds to the semi–classical regime.

An equivalent argument employs the fact that the path integral formulation of the quantum theory in Minkowski space sums over all field configurations weighted with exp(−*iSg*/¯*h*). Since *<sup>g</sup>*<sup>2</sup> *<sup>s</sup>* appears in the exponential in the same place as *h*¯, cf. Equation (18), this already suggests that for *g*<sup>2</sup> *<sup>s</sup>* → 0, the path integral is dominated by the classical configurations. Such configurations are believed to describe the matter inside incident nuclei during the initial stage of relativistic HICs at RHIC and the LHC [142,168].

Although in its original formulation, the QCD saturation is used for partons with fixed *Q* in case of the macroscopic bodies, it is more relevant to consider the partons at a fixed temperature *T* = *Q*/2*π*, avoiding at the same time the quantum entanglement problem [170], which is inevitably present in the description of microscopic objects. In this generalized setting, an object *A* filled with gluons may represent not only a fast-moving proton or nucleus but also the interior of the expanding early universe. Moreover, as follows from Equation (16), the saturation phenomenon is not necessary related to the growth of the gluon density at small *x*. For a big fast-moving domain of space filled with deconfined quarks and gluons, with radius *R*  1 fm, and hence also for the fast-expanding early universe itself characterized by the Hubble horizon *LH* ≫ 1 fm, the saturation limit can be reached even at *x* 1 [162] 2. In the extreme case, when the gluonic part of QCD matter completely decouples from the QCD fermionic fields and forms the vacuum condensate, the first term in Equation (9) can be neglected, and we obtain *g*QCD eff /*g*EW eff 8/3, making the CGC a prevailing form of matter during both the QCD and EW epochs.

Thus, in the periods of cosmological evolution when *T*  ΛQCD, including a very hot QCD era, EW era and beyond, it is perfectly conceivable that the universe was dominated by the fully saturated gluonic matter with occupation number *fg* <sup>∼</sup> *<sup>α</sup>*−<sup>1</sup> <sup>S</sup> . If during its subsequent cooling, the universe followed a trajectory in the ln(1/*x*), ln *Q*<sup>2</sup> plane staying still above the ln *Q*<sup>2</sup> *<sup>s</sup>* = *λY* line, see the right panel of Figure 5, the CGC phase would be a prevailing form of matter down to the temperatures *T* ≈ (2 − 5) × ΛQCD. For lower temperatures, the glasma is expected to fragment into a strongly interacting QGP.

The issue of emergence of classical behavior in the cosmological history has drawn recently a great deal of attention because of its conceptual as well as practical importance, see e.g., Ref. [171] and references therein. Although the origin of the observed anisotropies in the cosmic microwave background (CMB) radiation is traced back to vacuum fluctuations of quantum fields in the very early universe [84,120,172], there is a general expectation that the main characteristics of the universe can be described in classical terms even in its early history [171,173,174]. This is consistent with the fact that the initial conditions of the Hot Big Bang were determined by cosmic inflation driven by the so-called inflaton

field [84,86,175,176]. Apart from the quantum fluctuations of this field, its very emergence may well be a quantum phenomenon, e.g., the axion condensation (for a more thorough discussion, see below).

Let us note that the word "glass" appearing in the acronym CGC is used in condense matter physics to describe a non-equilibrium, disordered state of matter acting like solids on short time scales but liquids on long time scales [177]. In the glasma case, there are two scales present—the light cone time *τ*wee of low *x* (or wee) partons and the light cone time *τ*valence of primary (or valence) partons. For partons of transverse momentum *k*⊥, we obtain

$$\tau\_{\text{wee}} = \frac{1}{k^-} = \frac{2k^+}{k\_\perp^2} = \frac{2\varkappa P^+}{k\_\perp^2} \ll \frac{2P^+}{k\_\perp^2} \approx \tau\_{\text{valence}}.\tag{20}$$

suggesting that the valence parton modes are static over the times scales over which the wee modes are probed [167]. It is quite tempting to identify *τ*valence with the quantum break-time discussed in Refs. [178,179] defined as the time-scale after which true quantum evolution of the parton densities departs from the classical mean field evolution.

Glasses are formed when liquids are cooled too fast to form the crystalline equilibrium state. The fast cooling leads to an enormous number of possible configurations *N*gl(*T*) into which the glasses can freeze and consequently to their large entropy *S* = ln *N*gl(*T*), which does not vanish, even at *T* = 0, see e.g., Refs. [177,180]. In case of the CGC, the fast cooling is expected to take place in the Grand Unified Theory (GUT) era (see Section 2.5) when about one-third of all gauge bosons are gluons. By the end of that period, the excess of effective entropy DoF *h*eff is almost completely absorbed by the saturated gluonic matter.

The gluon condensation into (many domains of) the saturated phase was also facilitated by the fact that the wee partons "see" the color charge of other gluons over very large distances given by their transverse wavelength *<sup>λ</sup>*wee <sup>∼</sup> 1/*k*<sup>+</sup> <sup>=</sup> 1/(*xP*+). Since the glasma domains were formed in separate and completely different gluonic configurations, the saturated gluonic matter occupying the early universe had a substantial excess in the entropy DoF, *h*QCD eff , over the effective number of DoF in energy, *<sup>g</sup>*QCD eff . Consequently, the value of the generalized EoS parameter *w* was higher than that of the ideal massless gas; see Equation (6).

One possibility of how the EoS of a non-equilibrium matter comprising weakly interacting gluons can be approximated by the ideal massless gas of quasi-particles with *w*(*T*) > 1/3 follows from Equation (48) discussed later in Section 3.1. The glasma with *w*(*T*) ≈ 1/*D* may be looked upon as either a two-dimensional sheet (*D* = 2), a onedimensional string (*D* = 1) or, more generally, a fractal with the Hausdorff dimension 1 ≤ *D* < 3. Recent investigations of the dynamics of expanding glasma show that the spatial asymmetry introduced by the initial geometry is effectively transmitted to the azimuthal distribution of the gluon momentum field, even at very early times [181,182].

#### *2.5. The Running Couplings of the Standard Model and Their Unification*

The importance of QCD interactions in the EW era of the universe evolution can be also seen by comparing the corresponding couplings—the strong *α*S, electromagnetic *α*EM and weak *α*<sup>W</sup> ones. This is illustrated on the left panel of Figure 6 showing the RG flow in the scale *μ* = *Q* of the electromagnetic, weak and strong coupling parameters above *Q* = 100 GeV. Note that due to the fact that the gauge group of SM interactions SU(3)C × SU(2)L × U(1)Y is not a simple Lie group, the theory has not one but three coupling parameters, which are often denoted as *α*<sup>1</sup> ≡ *α*EM, *α*<sup>2</sup> ≡ *α*<sup>W</sup> and *α*<sup>3</sup> ≡ *α*<sup>S</sup> [183].

To see how the values of the coupling parameters influence the state of early universe, let us compare the collision time among its constituents *tc* ∼ 1/(*σnv*) (*σ* is the effective cross-section, *n* is the particle number density and *v* is their relative velocity) with the characteristic time-scale of the universe expansion *tH* ∼ 1/*H* [84]. Let us first restrict ourselves to the temperatures *T* - *T*EW when all particles of the SM are ultra-relativistic and the gauge bosons are massless. Then, the cross-sections for strong and EW interactions have a similar energy dependence *<sup>σ</sup>* <sup>∼</sup> *<sup>α</sup>*2/*T*2, where *<sup>α</sup>* <sup>10</sup>−<sup>1</sup> <sup>−</sup> <sup>10</sup>−<sup>2</sup> are the corresponding

dimensionless running couplings *α*1,2,3 varying logarithmically with *T*; see Figure 6. Taking into account Equations (39) and (42) for *<sup>n</sup>* <sup>∼</sup> *<sup>T</sup>*3, we find

$$t\_c \sim \frac{1}{\sigma nv} \sim \frac{1}{a^2 T} \ll t\_H \sim \frac{1}{H} \sim \frac{1}{\sqrt{\epsilon}} \sim \frac{1}{T^2} \,. \tag{21}$$

Thus, for the temperatures 10<sup>15</sup> <sup>−</sup> <sup>10</sup><sup>17</sup> GeV - *T* - *T*EW, the local equilibrium in the fluid is established before expansion of the universe becomes relevant.

At the lower temperatures down to *T*QCD *<sup>c</sup>* ≈ 160 MeV, the strong interactions prevail over the EW ones, and the local equilibrium is controlled by the interactions among the quarks and gluons in the strongly interacting QGP liquid. Even though the number density of the medium may be smaller3, the big effective cross-section among the particles forming the medium guarantees that the condition for local equilibrium *tc tH* is satisfied. Thus, over that whole range of the temperatures, the early universe is in the local equilibrium. It develops along the maximal possible entropy path, making it amenable to the hydrodynamical description.

In spite of its success, the SM cannot be the ultimate theory of particle physics. Such long-standing problems as the absence of a suitable DM candidate, no explanation of the observed baryon asymmetry in the universe, as well as various hierarchy problems in the underlined mass spectra (such as the flavor problem, the neutrino mass problem and the Higgs hierarchy problem) call for various bottom-up extensions of the SM framework as well as continuous attempts to derive the SM structure from a top-down perspective.

As was earlier discussed, e.g., in Ref. [184], the Higgs boson quartic coupling in the SM turns negative at scales -1010 GeV, rendering the vacuum state of the theory unstable at high energies. The current theoretical developments and experimental measurements suggest that the metastability of the Higgs vacuum is favored. This means a vacuum decay may occur with possibly catastrophic consequences for cosmology, since there are many catalysts that could trigger such a decay in the early universe. For a comprehensive review on cosmological implications of the Higgs vacuum metastability, see, e.g., Ref. [185].

Incidentally, almost immediately after the discovery of the asymptotic freedom, it was suggested [186] that at very high energies, the three gauge interactions of the SM are merged into a single force. The model, a first example of the Grand Unified Theory (GUT), is based on the smallest simple Lie group which contains the SM gauge groups SU(5) ⊃ SU*c*(3) × SU(2)L × U(1)Y. Among its 24 gauge bosons, there are in addition to eight gluons of QCD and four EW gauge bosons *W*±, *Z* and *γ* also 12 new ones called *X* and *Y*. Their emission or absorption makes it possible to transform a lepton into a quark or vice versa. Hence, the SU(5) GUT does not conserve baryon and lepton numbers separately, making it the first theory providing an explicit mechanism for the proton decay *<sup>p</sup>* <sup>→</sup> *<sup>e</sup>*+*π*0, with the half-time *<sup>τ</sup><sup>p</sup> <sup>M</sup>*<sup>4</sup> *X*/*m*<sup>5</sup> *<sup>p</sup>* <sup>10</sup><sup>30</sup> <sup>−</sup> 1031 years, where *MX* is the mass of SU(5) gauge boson at the scale of Grand Unification and *mp* is the mass of the proton. Although it was later found to disagree with experimental lower limit of *<sup>τ</sup><sup>p</sup>* <sup>≥</sup> 8.2 <sup>×</sup> <sup>10</sup><sup>33</sup> years [187] SU(5), unification is still considered an important example and a reference point of GUT model-building.

The basic property of the SU(5) theory and its later GUT successors [188–192] is that by virtue of the unification into a single (simple Lie) gauge group at very high energies, strict unification of the SM gauge couplings must take place. This is hinted at but not really achieved in the SM; see the left panel of Figure 6. First, *<sup>α</sup>*<sup>1</sup> and *<sup>α</sup>*<sup>2</sup> cross each other at *<sup>μ</sup>* <sup>∼</sup> 1013 GeV; then, *<sup>α</sup>*<sup>1</sup> crosses *<sup>α</sup>*<sup>3</sup> at *<sup>μ</sup>* <sup>∼</sup> <sup>10</sup><sup>14</sup> GeV, and finally, *<sup>α</sup>*<sup>2</sup> and *<sup>α</sup>*<sup>3</sup> cross at *<sup>μ</sup>* <sup>∼</sup> 1017 GeV, providing a hint of unification close to the Planck scale. Thus, as ultraviolet (UV) completions of the SM describing the physics at very high energies, GUTs are sometimes connected to theories of gravity such as string theory.

**Figure 6.** (**Left**) RG flow of the inverse SM gauge couplings *α*−<sup>1</sup> *<sup>a</sup>* as functions of the renormalization scale parameter *μ*. Index *a* = 1, 2, 3 stands for QED (*a* = 1), weak (*a* = 2) and QCD (*a* = 3) couplings. (**Right**) RG flow of the MSSM gauge couplings. Adapted from Ref. [192].

Let us now try to speculate on how the Grand Unification might have influenced the dynamics of the early universe at the temperatures when gluons decouple from quarks, forming a condensate. In SU(5) GUT, with increasing temperatures *<sup>T</sup> <sup>T</sup>*EW *<sup>c</sup>* , the gluon exchange between quarks (antiquarks) becomes overshadowed by the exchange of EW massless gauge bosons *W*±, *W*0, *B*<sup>0</sup> to be finally superseded by the GUT super-heavy *X* and *Y* gauge bosons. The latter also facilitate the conversion of quarks into leptons and vice versa. Since the transformations occur in the thermal and chemical equilibrium, the number of quarks and leptons remains constant on average.

One of the fundamental issues with the GUT models, which remains a challenge today, is the large hierarchy between the mass scale of Grand Unification and the EW scale. The latter is a source of large loop corrections to the Higgs mass. The problem is usually solved by extending the GUT with supersymmetry, the hypothetical symmetry between fermions and bosons (for a recent review on the concepts of supersymmetric GUTs, see e.g., Ref. [192]), which is broken at lower energies. Some of its minimal realizations, such as the Minimal Supersymmetric Standard Model (MSSM) [193], predict the unification of all three gauge couplings at the same scale; see the right panel of Figure 6.

Notwithstanding the drawbacks of GUTs, their cosmological signatures look quite promising. With the critical temperature *T*GUT *<sup>c</sup>* of the GUT phase transition approximately of the same order of magnitude as the particle masses at that temperature, the phase transitions that take place in GUTs at *T* - 10<sup>14</sup> GeV, as a rule, prove to be FOPTs [194]. Such transitions proceed via bubble nucleation [195]. While an isolated spherical bubble may produce GWs through sound waves in the plasma and magneto-hydrodynamics turbulence effects (see for example Refs. [196–199] and references therein), the process of bubble collision contributes to the GWs spectrum in the quadrupole approximation [86,200]. This contrasts with the SM phase transitions where the crossovers do not lead to a strong enhancement over the primordial GW spectrum. Moreover, the FOPTs also generate a primordial magnetic field during the turbulence phase of the plasma and bubble collision [86,201], and in some instances, they may generate topological defects such as domain walls and strings [195,202].

Let us add that an FOPT in the EW sector though precluded in the SM is possible in many of its scalar sector extensions [203]. In the most exotic scenario, a very peculiar history of the universe may occur: a first-order QCD phase transition (with six massless quarks) triggers an EW FOPT, which is eventually followed by a low-scale reheating of the universe where hadrons (likely) deconfine again, before a final, conventional crossover QCD transition to the current vacuum [204].

#### *2.6. Axions*

A promising avenue connecting the strong interactions with physics beyond the SM having at the same time far-reaching consequences in cosmology is provided by hypothetical ultra-light particles—the axions [183,205–208]. The QCD, unlike the EW interaction, is symmetric under time reversal and hence under a combined charge conjugation C and parity P operation CP. In principle, one can add to the QCD Lagrangian the term:

$$\mathcal{L}\_{\mathbb{Q}} = \theta \frac{\mathcal{g}\_s^2}{32\pi^2} F^{a}\_{\mu\nu} F^{\mu\nu}\_a \, , \tag{22}$$

where *F<sup>a</sup> μν* is the gluon field strength tensor (see Equation (17)), and *<sup>F</sup>*˜*μν <sup>a</sup>* = <sup>1</sup> <sup>2</sup> *μναβF<sup>a</sup> αβ* is its dual. For a non-zero value of the parameter *θ*, called the vacuum angle [209], the strong coupling permits violation of the CP symmetry. However, since L*<sup>Q</sup>* can be written as a total derivative of *<sup>K</sup>μ*, the Chern–Simons current, <sup>L</sup>*<sup>Q</sup>* <sup>=</sup> *∂μKμ*, the new term does not produce any effects in perturbation theory and is therefore usually neglected. Nevertheless, classical configurations, topological in nature, do exist, one example being the instantons [210], for which this term cannot be ignored. For instance, in the semi-classical dilute instanton gas approximation (DIGA) [211], the QCD vacuum energy density <sup>0</sup> depends on *θ* as

$$\varepsilon\_0(\theta) = -2\mathcal{C}\varepsilon^{-S\_{\rm inst}}\cos\theta \, \prime \, \, \, S\_{\rm inst} = \frac{8\pi^2}{\mathcal{S}\_s^2} \, . \tag{23}$$

Here, *C* is a positive constant and *S*inst is the QCD instanton action [207,208], cf. Equation (18). Moreover, since L*<sup>Q</sup>* in Equation (22) preserves the charge conjugation C, it contributes directly to the neutron electric dipole moment *dn* <sup>≈</sup> *emq*/*m*<sup>2</sup> *<sup>n</sup>θ*, where *e* is the proton charge, *mq* denotes the mass of *u*, *d* quarks, and *mn* is the neutron mass. Current measurements [212] provide an upper bound on the CP-violation parameter, <sup>|</sup>*θ*<sup>|</sup> <sup>10</sup>−10.

Within the SM, the smallness of the *θ* parameter becomes a true fine-tuning problem. Since *θ* could acquire an O(1) contribution from the observed CP-violation in the EW sector (via the common quark mass phase, arg det(*Mq*), where *Mq* is the quark mass matrix), its not obvious why it becomes cancelled to a high precision by the (unrelated) gluon term [207]. To solve this problem, the SM is augmented with an extra pseudo-scalar particle called axion *A*, whose only non-derivative coupling is to the CP-violating topological gluon density *F<sup>a</sup> μνF*˜*μν <sup>a</sup>* that is suppressed by a large scale *fA*. With *θ* → *θ* + *φ*(*x*)/ *fA*, where *φ* is the angular DoF of spin-zero complex field,

$$\varphi = |\varphi|e^{i\theta} = |\varphi|e^{i\phi/f\_A} \, , \tag{24}$$

the minimum of the vacuum energy occurs when the coefficient *θ* + *φ*/ *fA* in front of *F<sup>a</sup> μνF*˜*μν a* vanishes.

It is worth noting that interactions of the scalar field with ordinary matter is controlled by the factor *∂μϕ*/ *fA*. Thus, even at its originally suggested value of *fA* ∼ 250 GeV at the EW breaking scale [213], the axions interact so weakly that they emerge without attenuation from reactor cores or stellar interiors. Present astrophysical constraints push *fA* to substantially higher values, which are somewhere between a few 10<sup>8</sup> GeV and a few 10<sup>17</sup> GeV.

In cosmology, since their introduction, the axion-like particles were considered to be potentially important candidates if not for all than at least for the main component of the DM—a form of matter accounting for about one-quarter of its total energy density [72]. In order to fulfill their mission, the axions must contribute a non-negligible amount to the energy density of the universe and should have not been in thermal equilibrium with the cosmological plasma at any time in the history of the universe. This, together with the smallness of axion mass, implies their large occupation numbers [207]—the situation already encountered when we have discussed the properties of saturated gluon matter, cf. Equation (19). This implies that the axions over the whole history of the universe can be modeled by solving the classical field equations of a scalar condensate [214].

Starting from the Peccei–Quinn (PQ) scalar field *ϕ* introduced in Equation (24), the Lagrangian invariant under the global U(1)PQ transformation reads [215]:

$$\mathcal{L} = \frac{1}{2} |\partial\_{\mu} \varphi|^{2} - V\_{\text{eff}}(\varphi, T) \,, \quad V\_{\text{eff}}(\varphi, T) = \frac{\lambda}{4} (|\varphi|^{2} - f\_{A}^{2})^{2} + \frac{\lambda}{6} T^{2} |\varphi|^{2} \,. \tag{25}$$

Focusing on the early universe, the evolution of the field *ϕ* passes the following milestones. At high temperatures *<sup>T</sup> Tc* <sup>=</sup> <sup>√</sup><sup>3</sup> *fA*, the effective potential *<sup>V</sup>*eff(*ϕ*, *<sup>T</sup>*) depicted on the left panel of Figure 7 has the U(1)PQ symmetric minimum at *ϕ* = 0. With increasing time *t*, the universe cools, and the vacuum with *ϕ* = 0 becomes unstable. Due to the misalignment mechanism, the field starts to roll down from *ϕ* = 0, and the potential becomes tilted. At *T fA*, the PQ symmetry is spontaneously broken—the field acquires the vacuum expectation value *ϕ* = *fA*. Then, the axion—the Nambu-Goldstone boson of the spontaneously broken U(1)PQ symmetry—becomes a massless angular DoF at the minimum of the potential. The detailed results depend on whether the PQ phase transition occurs before or after inflation. While in the former case, only one *θ*<sup>0</sup> angle contributes (all other values are inflated away), in a post-inflationary scenario, the initial value of the angle *θ* takes all values in the interval −*π*, *π*. Eventually, when the universe cools to the temperatures of a few GeV, the axion obtains a mass through the QCD non-perturbative instanton effect known as the axial anomaly [183,216,217].

**Figure 7.** (**Left**) Potential of the PQ scalar *V*(*ϕ*) at different temperatures *T Tc* (pink) and *T Tc* (violet). The yellow circles show the positions of the minimum. Adapted from Ref. [215]. (**Right**) Continuum limit of *χ*top(*T*) from LQCD. The inserted sub-figure shows the behavior around the QCD phase transition temperature. Adapted from Ref. [89].

At the leading order in *f* <sup>−</sup><sup>1</sup> *<sup>A</sup>* , the axion mass *mA*(*T*) at some temperature *T* can be extracted from the QCD-generating functional Z(*θ*) in the presence of a theta term [218],

$$m\_A^2(T) = \frac{\delta^2}{\delta \phi^2} \ln Z\left(\frac{\phi}{f\_A}\right)\Big|\_{\phi=0} = \frac{1}{f\_A^2} \frac{d^2}{d\theta^2} \ln Z(\theta)\Big|\_{\theta=0} = \frac{\chi\_{\text{top}}(T)}{f\_A^2},\tag{26}$$

where *χ*top(*T*) is the QCD topological susceptibility. This quantity is typically computed using the lattice methods developed by several groups, see e.g., Refs. [89,219–228], but also in the framework of analytical approaches [218,229]. Its temperature dependence can be extracted either from the DIGA or from the LQCD calculations; see the right panel of Figure 7. The lattice simulations performed, for instance, in Ref. [89] have revealed that for *<sup>T</sup>* <sup>&</sup>gt; *<sup>T</sup>*<sup>∗</sup> <sup>=</sup> 150 MeV, the susceptibility falls as *<sup>χ</sup>*top(*T*) <sup>∼</sup> *<sup>T</sup>*−*<sup>b</sup>* with *<sup>b</sup>* <sup>=</sup> 8.16, extending thus the previous DIGA result *<sup>χ</sup>*top(*T*) = *<sup>χ</sup>*(0)*T*−<sup>8</sup> up to *<sup>T</sup>* <sup>≈</sup> 3 GeV. At the temperatures of *T* = 100 − 140 MeV in the vicinity of the QCD chiral phase transition temperature *Tc*, see Section 2.1, *χ*top(*T*) flattens. Further analysis performed in Ref. [89] exploiting the QCD EoS obtained therein has revealed that in the post-inflationary scenario, depending on the fraction of DM consisting of axions, *mA* = 50–1500 μeV. In particular, for the 50% axion content of the DM, one obtains the axion mass scale of *mA* = 50(4) μeV. Since at *T* < *Tc*, the chiral perturbation theory [183] becomes applicable, it is worth making a comparison with the original formula for the axion mass,

$$m\_A^2 = \frac{m\_\text{u} m\_d}{(m\_d + m\_d)^2} \frac{m\_\pi^2 f\_\pi^2}{f\_A^2} \implies m\_A \approx 5.7 \left(\frac{10^{12} \text{GeV}}{f\_A}\right) \mu \text{eV} \,\tag{27}$$

where *f<sup>π</sup>* is the pion decay constant, and *md* and *mu* are the down- and up-quark masses appearing in the QCD Lagrangian [183,208]. Let us note that if *mA* - 20 eV, the axions decay faster than the age of the universe.

For temperatures above the chiral phase transition, the axion potential computed in DIGA reads [211]:

$$V(\phi, T) = \chi\_{\text{top}}(T) \left[ 1 - \cos\left(\frac{\phi}{f\_A}\right) \right]. \tag{28}$$

Expanding Equation (28) around *φ*/ *fA* = 0, we obtain *V*(*φ*) = *m*<sup>2</sup> *<sup>A</sup>φ*2/2 at a finite *<sup>T</sup>*. Assuming the Friedmann-Lemaître-Robertson-Walker (FLRW) metric and classical axion field with this potential, the corresponding unperturbed energy density ¯*<sup>A</sup>* and pressure *p*¯*<sup>A</sup>* due to the axion field read

$$\bar{\epsilon}\_A = \frac{1}{2}\phi^2 + \frac{1}{2}m\_A^2\phi^2 \quad \text{and} \quad \bar{p}\_A = \frac{1}{2}\phi^2 - \frac{1}{2}m\_A^2\phi^2,\tag{29}$$

respectively. Substituting ¯*<sup>A</sup>* and *p*¯*<sup>A</sup>* from Equation (29) into the fluid Equation (39), we obtain

$$
\ddot{\Phi} + 3H\dot{\phi} + m\_A^2 = 0 \,. \tag{30}
$$

At early times *H*(*t*)  *mA*(*t*), we can neglect *mA* in Equation (30) to obtain the solution *φ*(*t*) = *φ*0—the axion field is frozen at a constant value *φ*<sup>0</sup> ∼ *fA*.

With increasing time *t*, eventually, the oscillating term proportional to *m*<sup>2</sup> *<sup>A</sup>*(*t*) ≡ *m*2 *<sup>A</sup>*(*T*(*t*)) in the equation of motion of the axion field (28) begins to contribute. At the time *t*osc defined implicitly as *mA*(*t*osc) ≈ 3*H*(*t*osc), the universe is sufficiently large to host a sizeable fraction of one oscillation period—the axion field starts to oscillate with an amplitude damped by the expansion rate. A solution of Equation (30) then reads [120],

$$\phi(t) = \phi\_1 \left(\frac{a(t\_1)}{a(t)}\right)^{3/2} \cos\left(\int\_0^t m\_A(t)dt + a\right),\tag{31}$$

where *t*<sup>1</sup> is the time at which *H*(*t*1) = *mA*, i.e., when the temperature drops below the QCD chiral phase transition temperature *Tc*, *φ*<sup>1</sup> ∼ *fA* is the constant and *α* is the phase. In particular, for *mA*(*t*) = *mA* and a radiation-dominated universe, one obtains *φ*<sup>1</sup> ≈ 1.44*φ*<sup>0</sup> and *α* = −3*π*/8 [120].

For the initial conditions at the onset of oscillations *<sup>θ</sup><sup>i</sup>* <sup>≡</sup> *<sup>θ</sup>*(*t*osc), ˙ *<sup>θ</sup><sup>i</sup>* <sup>≡</sup> ˙ *θi*(*t*osc), where *θ<sup>i</sup>* is called the initial misalignment angle, we obtain from Equation (30)

$$
\theta\_i = \theta\_{\rm PQ} + \frac{\phi\_{\rm PQ}}{H\_{\rm PQ}} \quad \text{and} \quad \theta\_i = \theta\_{\rm PQ} \left( \frac{H(t\_{\rm osc})}{H\_{\rm PQ}} \right)^{3/2} \tag{32}
$$

where *<sup>θ</sup>*(*t*PQ) ≡ *<sup>θ</sup>*PQ, *<sup>φ</sup>*˙(*t*PQ) ≡ *<sup>φ</sup>*˙ PQ, *H*PQ ≡ *H*(*t*PQ) and *a*PQ ≡ *a*(*t*PQ) are the values at the PQ symmetry breaking time *t*PQ *t*osc. In the second equation, we have also used *<sup>a</sup>* <sup>∼</sup> 1/*<sup>T</sup>* and *<sup>H</sup>* <sup>∼</sup> *<sup>T</sup>*2*G*1/2 [208].

While in the pre-inflationary scenario, inflation selects one patch of the universe within which the spontaneous breaking of the PQ symmetry leads to a homogeneous value of the initial misalignment angle *φi*, in the post-inflationary scenario, the PQ symmetry breaks with *θi*, taking different values in patches that are initially out of causal contact; see the left panel of Figure 7. However, today, they populate the volume enclosed by our Hubble horizon. In the post-inflationary scenario, the initial misalignment angle *φ<sup>i</sup>* takes all possible values on the unit circle. For a quadratic potential, *V*(*φ*) = *m*<sup>2</sup> *<sup>A</sup>φ*2/2, this is equivalent to an assumption that the initial condition reads *φ<sup>i</sup>* ≡ *φ*2 *<sup>i</sup>* = *π*/ <sup>√</sup>3, where the angle brackets represent the value averaged over −*π*, *π* [208].

Starting from *T*osc, the number of axions in a comoving frame becomes frozen, and their number density evolves as

$$m\_A(T\_{\rm osc}) \approx m\_A(T\_{\rm osc}) f\_A^2 \langle \phi\_i^2 \rangle \,. \tag{33}$$

For isotropic evolution, the ratio of the number density *nA* to the entropy density *s* in the comoving frame is conserved, i.e., *nA*(*T*)/*s*(*T*) = *nA*(*T*osc)/*s*(*T*osc) leading for *T Tc* to the expression for the axion energy density,

$$\epsilon\_A^{\rm miss} = m\_A n\_A (T\_{\rm osc}) \frac{h\_{\rm eff}(T)}{h\_{\rm eff}(T\_{\rm osc})} \left(\frac{T}{T\_{\rm osc}}\right)^3 = m\_A f\_A^2 \langle \phi\_i^2 \rangle \frac{h\_{\rm eff}(T)}{h\_{\rm eff}(T\_{\rm osc})} \left(\frac{T}{T\_{\rm osc}}\right)^3 T\_{\rm c} \,, \tag{34}$$

where the effective number of DoFs of the entropy *h*eff(*T*) is defined in Equation (3). Let us note that in contrast to Ref. [208] in Equation (34), the constancy of the axion mass for *T Tc* is already taken into account.

After the spontaneous U(1)PQ symmetry breaking the axion field, *φ*, being an angular variable, takes values in the interval 0, 2*π fA*, cf. Equation (24). Consequently, the axion potential *V*(*φ*) given by Equation (28) is periodic in *φ* with period Δ*φ* = 2*π fA*/*N*DW. In the other words, *V*(*φ*) has an exact *ZN*DW discrete symmetry. The axion acquires a periodic potential with *N*DW equivalent minima. The Kibble mechanism [195,230] then dictates that, depending on the homotopy group *π*(M) of the manifold M of degenerate vacua, the topological defects—domain walls, strings or monopoles—form each time the symmetry is broken [195,202]. With M = U(1)PQ and *π*(M = *ZN*DW ), the production of axionic strings, which are vortex-like topological defects that form as soon as the symmetry is spontaneously broken, is possible. Those are not important when the PQ symmetry is broken before inflation—they are inflated away—but they play an important role in the post-inflationary scenario. When the Hubble parameter *H* becomes comparable to the axion mass *mA*, the axion starts to roll down to one of the minima. Since the axion field settles into different minima in different places of the universe, domain walls are formed between the different vacua; see the left panel of Figure 7. It is worth mentioning that this phenomenon is similar to ice formation on the surface of a pond or a puddle when the water begins to freeze in many places independently, and the growing plates of ice join up in random fashion, leaving zigzag boundaries between them [195].

As an example, let us consider a planar wall orthogonal to the *z*-axis *φ* = *φ*(*z*). The solution of the classical field equation with potential given by Equation (28) reads [231]:

$$\frac{\phi(z)}{f\_A} = \frac{2\pi k}{N\_{\rm DW}} + \frac{4}{N\_{\rm DW}} \tan^{-1} e^{m\_A z} \,. \tag{35}$$

This configuration interpolates between the two allowed vacua, *φ*/ *fA* = 2*πk*/*N*DW at *z* → −∞ and *φ*/ *fA* = 2*π*(*k* + 1)/*N*DW, which are separated by the wall of thickness 1/*mA*.

Astrophysical signatures of the axion can be broadly divided into its couplings to elementary/composite particles, i.e., photons, electrons, protons or neutrons, and to the macroscopic objects in the universe such as BHs [208]. In the latter case, when the Compton length of the axions becomes of order of the BH size, they form gravitational bound states around it. The phenomenon of superradiance [232,233] causes the axion occupation numbers to grow exponentially, providing a way to extract very efficiently energy and angular momentum from the BH. The presence of axions could be inferred by observations of BH masses and angular momenta. Current measurements exclude the region of 6 × <sup>10</sup><sup>17</sup> GeV <sup>≤</sup> *fA* <sup>≤</sup> <sup>10</sup><sup>19</sup> GeV.

In particle physics signatures, the most important is the decay channel of axion into two photons with the decay width <sup>Γ</sup>A→*γγ* <sup>=</sup> *<sup>g</sup>*<sup>2</sup> A*γγm*<sup>3</sup> *<sup>A</sup>*/(64*π*) in terms of the model-dependent coupling constant *g*A*γγ*. The main uncertainty is due to the electromagnetic and color anomalies of the axial current associated with the axion. The most relevant process induced by the *g*A*γγ* is the Primakoff process—the conversion of thermal photons in the electrostatic field of electrons and nuclei into axions,

$$
\gamma + Z\varepsilon \to A + Z\varepsilon \,. \tag{36}
$$

A strong bound on exotic cooling processes in the sun is provided by the helioseismological considerations [234] giving the following constraint: *<sup>g</sup>*A*γγ* <sup>≤</sup> 4.1 <sup>×</sup> <sup>10</sup>−<sup>10</sup> GeV<sup>−</sup>1. For more details, consult Refs. [205–208,215]. Another option is the inverse Primakoff scattering, which allows solar axions to coherently scatter off the atomic electric field and back-convert into photons in the detector volume, *A* + *Ze* → *γZe*, proceeding through a *t*-channel photon exchange [235].

The axion besides being a well-recognized DM candidate may provide also an interesting explanation of a wider range of phenomena related to the early universe dynamics. As a first example, let us mention the SMASH model—a minimal extension of the SM with additional particle content comprising three sterile right-handed neutrinos *Ni*, *i* = 1, 2, 3, a color triplet *<sup>Q</sup>* and a complex SM-singlet scalar *<sup>σ</sup>*, whose VEV of *<sup>v</sup><sup>σ</sup>* <sup>∼</sup> 1011 GeV breaks the lepton number and the PQ symmetry simultaneously [236]. At low energies, the model reduces to the SM, which is augmented by seesaw-generated neutrino masses and mixing, plus the axion. In this scenario, the inflaton, a scalar field driving cosmic inflation in the very early universe, is a mixture of *σ* and the SM Higgs fields. The reheating of the universe after inflation occurs via a mechanism known as a Higgs portal [237]—by the DM particles, which interact only through their couplings with the Higgs sector of the theory. The model provides a consistent picture of particle physics from the EW scale to the Planck scale *M*PL and of cosmology from inflation until today. In particular, in the SMASH model framework, the PQ symmetry is first broken and then restored non-thermally during preheating for *fA* <sup>=</sup> <sup>4</sup> <sup>×</sup> 1016 GeV.

The second example, the Axion Quark Nugget (AQN) DM model, see, e.g., Ref. [238] and references therein, replaces the commonly accepted baryogenesis scenario with a charge separation process in which the global baryon number of the universe remains zero at all times. Similarly to Witten's idea of stranglets [5], the AQN DM is composed of quarks and anti-quarks but now in a new high-density CSC phase. Initially, nuggets of both matter and antimatter are formed with equal probability as a result of the dynamics of the axion domain walls which at the same time provide the extra pressure needed to stabilize the CSC phase. Later on, due to the global CP violating processes associated with the initial misalignment angle *θ*<sup>0</sup> = 0 during the early formation stage, the populations of the nuggets with the positive and negative baryon number become different. The unobserved antibaryons hidden inside the DM would not participate in nucleosynthesis and, therefore, according to the usual definition would not contribute to the visible matter. However, since antimatter nuggets can interact with regular matter via annihilation leading to electromagnetic radiation, their existence has observational consequences [238].

It is worthwhile mentioning here an interesting generalization of the PQ mechanism where in addition to *θ* angle, also the strong coupling *α*<sup>S</sup> is promoted to a dynamical quantity. The latter evolves through the VEV of a singlet scalar field that mixes with the Higgs field [239,240]. In the resulting cosmic history, the QCD confinement and EW symmetry breaking initially occur simultaneously close to the weak scale.

To conclude this section, as we have already noticed in several examples above, the nonperturbative dynamics of QCD often exhibits very non-trivial and rather unexpected consequences at cosmological scales. An attractive mechanism proposing that the QCD axion may emerge as a composite state has been discussed very recently in Ref. [241]. In particular, it was suggested that Majorana neutrinos, that combine into Cooper pairs, can form collective low-energy degrees of freedom. This motivates the existence of the

QCD axion as a collective excitation of the neutrino condensate. Such a condensate can be produced after the QCD phase transition epoch in a cold coherent state by means of a misalignment mechanism, thus providing an alternative DM candidate. In this case, a QCD anomalous portal provides the necessary means for a tiny mass gap generation by neutrinos. Furthermore, the Cosmological Constant emerges as a result of the spontaneously broken mirror symmetry of the QCD ground state triggered by the quantum gravity effects as suggested by Refs. [65,66] (see also below). Hence, one concludes that QCD may be responsible for the dynamical generation of both the DM and DE components of the universe such that a complete knowledge of the QCD in the infrared regime may be absolutely critical for understanding of the cosmological evolution since the latest QCD transition epoch and for eventual formation of the current state of the universe.

#### **3. Dynamics of the Early Universe**

#### *3.1. Simple Models with Constant Speed of Sound*

In accord with the observations, the Standard Cosmological Model (SCM) postulates that the cosmic matter at scales larger than 100 Mpc is homogeneously and isotropically distributed. Consequently, thermodynamic pressure *p* at early times of its evolution depends on temperature *T* and various chemical potentials *μi*, *i* = *B*, *Q*, *L*, ... corresponding to baryon number *B*, electric charge *Q*, lepton number *L* etc., only via the energy density . Solution of the Einstein equations of general relativity

$$R\_{\mu\nu} - \frac{1}{2}g\_{\mu\nu}(R - 2\Lambda) = -8\pi GT\_{\mu\nu} \tag{37}$$

where *Rμν* is the Ricci tensor, *R* = *Rμνgμν* is the scalar curvature, *gμν* is the metric tensor, and Λ and *G* are the cosmological and gravitational constants; preserving the homogeneity and isotropy of space under its time evolution is a spacetime of constant curvature parameter *k* = {+1, 0, −1}. It is described by a single function—the time-dependent scale factor *a*(*t*) [84,120,242] which connects the Lagrangian (or comoving) coordinates *r* with the physical Euler coordinates *r*ˆ(*t*) = *a*(*t*)*r*. The metric tensor *gμν* in the preferred coordinate system where these symmetries are clearly manifest reads4

$$ds^2 = \mathbf{g}\_{\mu\nu} \mathbf{x}^{\mu} \mathbf{x}^{\nu} = dt^2 - a^2(t) \left[ \frac{dr^2}{1 - kr^2} + r^2(d\theta^2 + \sin^2\theta d\phi^2) \right], \qquad k = \text{const.} \tag{38}$$

Using this metric in Equations (37) and (38) and neglecting the dissipative terms yields the Friedmann equation for the time evolution of *a*(*t*) and the fluid equation for the time evolution of (*t*), respectively (see e.g., Refs. [84,175]),

$$H^2(t) \equiv \left(\frac{\dot{a}}{a}\right)^2 = \frac{8\pi G}{3}\epsilon - \frac{k}{a^2} + \frac{\Lambda}{3}, \qquad \dot{\epsilon} + 3(\epsilon + p)H(t) = 0,\tag{39}$$

where *H*(*t*) is the Hubble parameter. From Equation (39), it follows that the expanding universe is characterized by a natural time-scale *H*−<sup>1</sup> = *a*/*a*˙. Any particle species will remain in thermal equilibrium with the cosmic fluid so long as the mean interaction time *tc* allows rapid adjustment to the falling temperature provided that *tc* < *H*−1.

In the period of the universe evolution when - 1 GeV fm−3, the terms containing constants Λ and *k* in Equation (39) can be safely neglected, transforming the above two equations into a single one describing the time evolution of the energy density [6],

$$-\frac{d\varepsilon}{3\sqrt{\varepsilon}(\varepsilon+p)} = \sqrt{\frac{8\pi G}{3}}dt\,. \tag{40}$$

For the time-independent speed of sound *cs* and for the energy densities negligible compared to the initial density (*t* > 0) (*t* = 0) integration of Equation (40), the calculations yield [13]

$$\epsilon(t) = \frac{1}{6\pi G (1 + c\_s^2)^2 t^2}, \qquad c\_s^2 \equiv \frac{dp}{d\epsilon} = \text{const.} \tag{41}$$

Substituting this equation into the fluid Equation (39), we obtain the expansion rate of the early universe

$$
\dot{a} \sim t^{-\kappa}, \qquad \kappa = \frac{1 + 3c\_s^2}{3(1 + c\_s^2)}. \tag{42}
$$

In particular, for the massless non-interacting gas with *c*<sup>2</sup> *<sup>s</sup>* = 1/3, one obtains *a*˙ ∼ *t* −1/2 and hence *a*(*t*) ∼ *t* 1/2.

It is worth mentioning that in the cosmological literature, see, e.g., Ref. [120], it is customary to express the EoS in terms of the parameter *w*, Equation (5). For the constant speed of sound case, i.e., for *w* = const, the solution of the Friedmann Equation (39) yields

$$
\epsilon \sim a^{-3(1+w)}.\tag{43}
$$

For the non-relativistic matter (dust) with *<sup>p</sup>* <sup>=</sup> 0, we obtain <sup>∼</sup> *<sup>a</sup>*−3, for the massless non-interacting gas, <sup>∼</sup> *<sup>a</sup>*−4, for the EoS *<sup>p</sup>* <sup>=</sup> corresponding to absolutely stiff fluid [243], <sup>∼</sup> *<sup>a</sup>*−6, and for the vacuum energy with the EoS *<sup>p</sup>* <sup>=</sup> <sup>−</sup>, Equation (43) gives <sup>=</sup> const.

Let us now follow Ref. [244] and consider an ideal gas of free particles in *D*-dimensional space. Its particle number density *n*, energy density and pressure *p* expressed in terms of single-particle statistical sum *f*(*E*, *T*, *μ*) of particle with energy *E*, momentum *P* and spin *s* reads (see e.g., Ref. [245]):

$$f(E, T, \mu) = \frac{1}{\exp[(E - \mu)/T] \pm 1}, \qquad n = \gamma \int f(E, T, \mu) d^D P,\tag{44}$$

$$\varepsilon = \gamma \int f(E, T, \mu) E(P) d^D P = \gamma S(D) \int\_0^\infty f(E, T, \mu) E(P) P^{D-1} dP,\tag{45}$$

$$p = -\frac{T}{V} \ln Z = -\gamma T \int \ln f(\mathcal{E}, T, \mu) d^D P = \gamma \frac{S(D)}{D} \int\_0^\infty f(\mathcal{E}, T, \mu) \frac{\partial \mathcal{E}(P)}{\partial P} P^D dP,\tag{46}$$

where *<sup>d</sup>DP* = [*DπD*/2]/[Γ(*D*/2 <sup>+</sup> <sup>1</sup>)]*PD*−1*dP* <sup>≡</sup> *<sup>S</sup>*(*D*)*PD*−1*dP* is a volume element of the *<sup>D</sup>*-dimensional hypersphere and *<sup>γ</sup>* <sup>≡</sup> (2*<sup>s</sup>* <sup>+</sup> <sup>1</sup>)(2*π*)−*D*. When evaluating the first integral in Equation (46), we have performed integration by parts assuming that the particle energy *E*(*P*) is some generic function of *P*. The substitution of Equations (45) and (46) into Equation (5) constrains *E*(*P*) to satisfy the differential equation

$$\frac{P}{D} \frac{\partial E(P)}{\partial P} = wE(P) \, , \tag{47}$$

whose solution for *w* = const reads

$$E(P) = \mathfrak{F}P^{wD},\tag{48}$$

with *ξ* some arbitrary constant. Hence, the medium with constant speed of sound squared *c*2 *<sup>s</sup>* = *w* in ordinary three-dimensional space can be equivalently described as an ideal gas of quasi-particles with energy *E* and momentum *P* satisfying the dispersion relation (48) in *D*-dimensional space [244]. It is worth mentioning that at some instances, the dispersion relation (48) can be satisfied by the real particles. The case of *w* = 0 corresponds to nonrelativistic particles, while *wD* = 1 corresponds to massless particles in *D*-dimensional space with the EoS *p* = /*D* and hence with the sound velocity *cs* = *D*−1/2.

Unfortunately, this is not the case of the absolutely stiff fluid, as first discussed by Zeldovich [243]. Notwithstanding that its EoS *p* = can be used to describe a large variety of systems such as phonon-like excitations in a thin channel (*D* = 1) [244], thin film (*D* = 2) of non-relativistic quasi-particles [244], interiors of neutron stars [246,247], Big Bang nucleosynthesis [248] or warm self-interacting DM component [249]. The cosmology and thermodynamics of the FLRW universe with bulk viscous stiff fluid was studied in Ref. [250]. The fact that the stiff fluid saturates the holographic covariant entropy bound was used in Ref. [251] to describe a cosmology of the very early universe. Last but not least, a stiff perfect fluid is energetically equivalent to a time-like massless scalar field *φ*, see, e.g., Ref. [252]. From its energy-momentum tensor,

$$T^{\phi}\_{\mu\nu} = \partial\_{\mu}\phi \partial\_{\nu}\phi - \frac{1}{2}g\_{\mu\nu}\partial\_{\mu}\phi \partial^{\mu}\phi \ , \qquad \partial\_{\alpha}\phi \partial^{\alpha}\phi > 0 \tag{49}$$

we obtain that *<sup>p</sup>* <sup>≡</sup> *<sup>T</sup><sup>φ</sup> kk* <sup>=</sup> <sup>≡</sup> *<sup>T</sup><sup>φ</sup>* <sup>00</sup> = *∂αφ∂αφ*.

#### *3.2. Equation of State of the Early Universe*

As already discussed in Section 3.1, the only way how to write the EoS compatible with homogeneity and isotropy of the universe is through the pressure expressed as a function of the energy density. Thus, given the barotropic EoS *p*() of the expanding matter, we can for instance use Equation (40) to predict the temporal evolution of the energy density (*t*) and the temperature *T*(*t*). For the ideal gas EoS, Equation (41) yields the relation between the temperature *T* of the early universe and the time *t*sec elapsed from the Big Bang, *T*MeV =O(1)/ <sup>√</sup>*t*sec [84], i.e., the same time-dependence as for the time derivative of the scale factor *a*˙(*t*); see Equation (42) and the discussion below it. For the QGP-to-hadronic matter transition, this leads to *t* QGP sec <sup>∼</sup>10−<sup>5</sup> *<sup>s</sup>*, and for the EW phase transition, this leads to *t* EW sec <sup>∼</sup>10−<sup>11</sup> *<sup>s</sup>*.

More sophisticated EoS can be based either on the phenomenological models or on the microscopic theory. In the latter case, a quantum physics formulation of statistical mechanics in terms of the *S* matrix, which describes the scattering processes taking place in the thermodynamical system of interest, is available—see Ref. [253]. It provides a simple prescription for calculating the grand canonical potential *Z*(*T*, *μ*) of any gaseous system given the free-particle energies and *S*-matrix elements. The application of *S* matrix formulation to study the thermal properties of an interacting gas of hadrons can be found in Ref. [254]. It can be also used to show how the hadron resonance gas model emerges from the *S*-matrix framework [255]. However, this approach is based on the perturbative expansion of the *S*-matrix and is therefore not applicable in the strong coupling regime of QFT. There, one resorts to direct calculation of the grand canonical potential on the lattice.

Our first example of the EoS used in a description of the evolution of the early universe is the Bag Model (BM) EoS [1,256,257] based on the phenomenological description of the mass spectrum of the hadron states [11,258,259] in terms of gas of massless color objects quarks and gluons—moving inside the confining potential—the bag,

$$\varepsilon\_q(T) = \sigma\_q T^4 + \mathcal{B}, \qquad p\_q(T) = \frac{\sigma\_q}{3} T^4 - \mathcal{B}, \qquad p\_q(\epsilon) = \frac{1}{3} (\epsilon - 4\mathcal{B}) \,. \tag{50}$$

In Equation (50), *<sup>σ</sup><sup>q</sup>* = *<sup>π</sup>*<sup>2</sup> <sup>30</sup> *<sup>g</sup>*QCD eff is the Stefan-Boltzmann constant with *<sup>g</sup>*QCD eff given by Equation (9). The BM EoS (50) incorporates color confinement through the bag constant B = bag − vac > 0, indicating the difference between the energy densities of the physical vacuum and the ground state for quarks and gluons in the medium. The latter can be interpreted as the energy needed to create a bubble in the vacuum in which the noninteracting quarks and gluons are confined. While the fit to hadron masses made in the original BM predict <sup>B</sup>1/4 <sup>≈</sup> 140 MeV, the value of <sup>B</sup>1/4 <sup>≈</sup> 220 MeV is frequently quoted in works dealing with the vacuum structure of QCD, see, e.g., Ref. [11].

Let us note that BM EoS represents a bare-bones model of hadron-to-QGP phase transition. At small energy densities, hadrons—the bubbles inside the non-perturbative vacuum—occupy only a small fraction of the total considered volume *V*. An increase of leads to the coalescence of several bubbles into larger ones. For ≥ 4B, the volume *V* is filled with one large bubble, cf. Equation (50), whose surface coincides with the enclosing walls. Hence, there is no longer free surface against the vacuum, and the new de-confined phase of matter is created [256]. The critical temperature *Tc* of the phase transition can be estimated using Gibbs criteria. Equating the BM pressure (50) with the pressure *ph* of the hadron (pion) gas with 3 DoF and hence with *σ<sup>h</sup>* = 3*π*2/30, we obtain

$$p\_{\boldsymbol{q}}(T\_{\boldsymbol{c}}) = \frac{\sigma\_{\boldsymbol{q}}}{3} T\_{\boldsymbol{c}}^{4} - \mathcal{B} = p\_{\boldsymbol{h}}(T\_{\boldsymbol{c}}) = \frac{\sigma\_{\boldsymbol{h}}}{3} T\_{\boldsymbol{c}}^{4}, \qquad T\_{\boldsymbol{c}} = \left(\frac{3\mathcal{B}}{\sigma\_{\boldsymbol{q}} - \sigma\_{\boldsymbol{h}}}\right)^{1/4}.\tag{51}$$

For *NF* <sup>=</sup> 3 active quark flavors *<sup>u</sup>*, *<sup>d</sup>*, and *<sup>c</sup>*, Equation (2) yields *Tc* <sup>≈</sup> 0.67B1/4 and ≈150 MeV. Let us add that Equation (51) describes the first-order phase transition with energy density discontinuity,

$$
\Delta \varepsilon = \varepsilon\_{\mathbb{H}}(T\_{\varepsilon}) - \varepsilon\_{\mathbb{H}}(T\_{\varepsilon}) = 3p\_{\mathbb{H}}(T\_{\varepsilon}) + 4\mathcal{B} - 3p\_{h}(T\_{\varepsilon}) = 4\mathcal{B} \,. \tag{52}
$$

Last but not least, expressing the BM EoS (50) in terms of the dimensionless interaction measure—the trace anomaly describing the thermal contribution to the trace of the energymomentum tensor <sup>T</sup> *μν* <sup>≡</sup> *<sup>T</sup>μν*

$$\Theta \equiv \frac{\mathcal{T}^{\mu\mu}(T)}{T^4} = \frac{\epsilon\_q - 3p\_q}{T^4} = \frac{4\mathcal{B}}{T^4} = \frac{4\sigma\_q \mathcal{B}}{\epsilon\_q - \mathcal{B}},\tag{53}$$

we observe a monotonous weakening of the interaction strength with increasing *<sup>q</sup>* [33], leading ultimately to a Stefan-Boltzmann (SB) value of Θ = 0. At the same time, the entropy density *sq* = (*<sup>q</sup>* + *pq*)/*T* = 4/3*σ*1/4 *<sup>q</sup>* (*<sup>q</sup>* − B)3/4 converges to its SB limit from below. The speed of sound derived from the BM EoS (50) is energy-density-independent and coincides with that of the ideal gas of massless particles *c*<sup>2</sup> *<sup>s</sup>* = *dpq*/*d<sup>q</sup>* = 1/3.

More sophisticated EoS can be constructed, e.g., by adding the term <sup>∼</sup> *<sup>T</sup>*<sup>2</sup> to expressions for (*T*) and *p*(*T*) in Equation (50) with *σ* = *σ<sup>q</sup>*

$$
\varepsilon(T) = \sigma T^4 - CT^2 + \mathcal{B} \,, \qquad p(T) = \frac{\sigma}{3} T^4 - DT^2 - \mathcal{B} \,. \tag{54}
$$

The motivation for such a modified BM EoS comes from the observation [260,261] that in the case of a pure gauge theory up to temperatures a few times the transition temperature *Tc*, the dominant power-like correction to the pQCD high-temperature behavior is <sup>O</sup>(*T*−2) rather than <sup>O</sup>(*T*−4). Moreover, the quadratic thermal terms in the deconfined phase can be also obtained from gauge/string duality [262].

In its original setting with *C* = *D* > 0, Equation (54) represents the LQCD-motivated "fuzzy" BM EoS of Ref. [263]. On the other hand, for *C* = −*D* < 0, it represents a particular case of gas of gluonic quasi-particles EoS with a temperature-dependent bag function <sup>B</sup>(*T*) = <sup>−</sup>*CT*<sup>2</sup> <sup>+</sup> <sup>B</sup>; see e.g., Refs. [264,265]. In the following discussion, we will keep the values of the constants *C* and *D* in Equation (54) unrestricted.

By inverting (*T*) in Equation (54) with respect to temperature squared *T*2,

$$T^2(\varepsilon) = \frac{\mathbb{C} + \sqrt{\mathbb{C}^2 + 4\sigma(\varepsilon - \mathcal{B})}}{2\sigma} > 0 \,, \tag{55}$$

and substituting for *T*2() into Equation (54), we obtain the barotropic form of the EoS (54):

$$p(\epsilon) = \frac{1}{3}(\epsilon - 4\mathcal{B}) - \frac{1}{3}\text{sgn}(A)|A|T^2(\epsilon) \quad A = 3D - \mathcal{C} \,. \tag{56}$$

The corresponding sound velocity squared reads

$$c\_s^2(\varepsilon) = \frac{dp(\varepsilon)}{d\varepsilon} = \frac{1}{3} \left( 1 - \frac{\text{sgn}(A)|A|}{\sqrt{\mathbb{C}^2 + 4\sigma(\varepsilon - \mathcal{B})}} \right). \tag{57}$$

Note that for *A* = 0, both Equations (56) and (57) degenerate to the corresponding result from the BM EoS; see Equation (50). In the context of Equation (54), the situation with *C* = 3*D* = <sup>3</sup> <sup>2</sup> *NFμ*<sup>2</sup> *<sup>B</sup>* describes the EoS of ideal QGP with non-zero baryon chemical potential [266]. It is worth mentioning that by tuning the phase transition temperature to *Tc* = 160 MeV, this EoS predicts a six times higher value of the bag constant B than the original BM [266].

Therefore, in the following, we discuss only the non-trivial case of *A* = 0. First, we check which values of the constants *C* and *D* in Equation (57) are compatible with the condition *c*<sup>2</sup> *<sup>s</sup>* > 0. While for *A* < 0 (i.e., *C* > 0 and −*C* < *D* < *C*/3), Equation (57) is always satisfied for *A* > 0 (*C* ≥ 0 and *D* > *C*/3 or *D* > −*C* > 0), there exists a lower bound <sup>0</sup> on energy density

$$
\varepsilon > \varepsilon\_0 = \frac{A^2 - C^2}{4\sigma} + \mathcal{B} \,. \tag{58}
$$

Thus, in both cases (*A* > 0 or *A* < 0), Equation (55) represents the genuine nontrivial EoS of non-ideal high-density matter each with its own sound velocity approaching for <sup>→</sup> <sup>∞</sup>, either from below or from above, the SB limit of <sup>√</sup>1/3. However, only for *<sup>A</sup>* <sup>&</sup>lt; 0, the second term on the right-hand side of Equation (56) with <sup>−</sup><sup>1</sup> <sup>3</sup> *AT*2() represents independent pressure. It is also worth mentioning that for the latter case (*A* < 0), the trace anomaly

$$\Theta(\epsilon) = \frac{4\mathcal{B}}{T^4(\epsilon)} + \frac{\text{sgn}(A)|A|}{T^2(\epsilon)}\,'\,. \tag{59}$$

with *T*2() defined in Equation (55), acquires a peak at the energy density

$$\epsilon\_p = \frac{2\mathcal{B} \pm \sqrt{2\mathcal{B}|A|}\mathcal{C}}{|A|\sigma} - \mathcal{B} \,. \tag{60}$$

Standard explanation of this phenomenon within the SU(3)*<sup>c</sup>* gauge theory, see, e.g., Ref. [267], relies on the fact that in the region around and just above the critical temperature *Tc* of hadron-to-QGP phase transition, the energy density rises much more rapidly than the pressure, leading to the observed rapid increase of Θ. Since asymptotically /*T*<sup>4</sup> and 3*p*(*T*)/*T*<sup>4</sup> converge to their common Stefan-Boltzmann value of *σ*, see Equation (54), there must be some temperature *Tp* (and hence also some energy density *<sup>p</sup>* = (*Tp*)) at which the growth rates change roles, with the pressure now increasing more rapidly. The further decrease of Θ is in good approximation given by *T*−2, so that *T*2Θ(*T*) becomes approximately constant very soon above *Tc* and up to about 5*Tc* [263].

Note that Equation (40) with the initial condition 0(*t*0) = 104 GeV fm−<sup>3</sup> at *t* QGP sec *<sup>t</sup>*<sup>0</sup> <sup>=</sup> <sup>10</sup>−<sup>9</sup> <sup>s</sup> *<sup>t</sup>* EW sec was used in Ref. [13] to study the sensitivity of dilution and cooling of the early universe to the changes in the primordial QGP EoS. The latter might modify the pattern of the emission of GWs or the generation of baryon number fluctuations. No dramatic changes between different EoS, including Equations (50) and (54) (with *C* = *D*) and others, were found in the whole considered time interval. This finding seems to be supported by the LQCD calculations which show that the transition from primordial QGP matter to hadronic matter proceeded as a continuous crossover [92]. The latter does not introduce any fluctuations on length scales much longer than the natural length scales of QCD <sup>∼</sup> <sup>Λ</sup>−<sup>1</sup> QCD, so it has probably left no imprint in the microseconds-old universe that survived so as to be visible in some way today [7].

A different conclusion was, however, reached in Ref. [74] where a novel mechanism for the production of GWs during the QCD phase transition has been proposed. It was found that while the energy density of the homogeneous gluon condensate is smoothly decaying in cosmological time *t*sec, the pressure of the condensate undergoes a sequence of violent oscillations at the characteristic QCD time scales—the process called relaxation—the basic feature of the QCD transition. Such relaxation processes generate a very specific multi-peaked GWs signature in the domain of radio frequencies. In particular, gravitational echoes of the QCD transition potentially accessible by the FAST [268] and SKA [269] GW telescopes are sourced by the formation of domain walls in the QCD vacuum (also responsible for color confinement in the infrared regime). More details about the dynamical theory of the YM vacuum will be provided below in Section 4.

Recently, the stochastic GWs induced by the scalar perturbations of the metrics were investigated in Ref. [270] as a cosmological probe of the sound speed *cs* during the QCD phase transition. Although the GW propagation itself does not depend on *cs*, the sound speed value affects the dynamics of primordial density. Induced stochastic GWs can thus be an indirect probe of both the EoS parameter *w* = *p*/ and *cs* = *dp*/*d*. In particular, similarly to the conclusions of Ref. [74], the GW frequency <sup>∼</sup> <sup>10</sup>−<sup>8</sup> Hz corresponding to the Hubble scale during the QCD phase transition appears to be in the range of the planned GW detectors.

Let us now turn to a description based on the fundamental theory. Using the Lagrangian of the SM, one can extract the thermodynamical quantities either directly from the lattice calculations [89], deduce them from lattice simulations using a dimensionallyreduced EFT [90] or use the perturbation theory. With the temperature *T* depending on the lattice spacing *a* and the number of lattice points in the temporal direction *Nt* as *T* = (*aNt*)−1, the reliable Monte Carlo simulations of QFT in the high-temperature regime need very fine lattices. At the same time, as the lattice spacing *a* is reduced, the autocorrelation times for zero temperature simulations rise, and the costs of these simulations explode beyond feasibility [89]. This makes the LQCD simulations in the region of temperatures higher then the few GeVs practically impossible. Alternatively, one can vary the gauge coupling *gs*, which leads to changing *T* as well, although the spacial and temporal dimensions do not [90]. In addition to that, even at the temperatures of the EW phase transition *T*EW *<sup>c</sup>* , the dynamics needs to be treated with lattice methods. This is due to the fact that when the particle momenta in the range *<sup>k</sup>* <sup>∼</sup> *<sup>g</sup>*2*T*/*<sup>π</sup>* are considered, then the dynamics of the system is non-perturbative [194,271,272]. Values much below and far above this value *p*/*T*<sup>4</sup> can be determined by a direct perturbative computation [90].

In the SM framework, the basic thermodynamical observable is the pressure [273], which can be formally defined through the grand canonical partition function Z as

$$Z \equiv \exp\left[\frac{p\_B(T)V}{T}\right], \quad p\_B(T) = p\_E(T) + p\_M(T) + p\_G(T) \,. \tag{61}$$

where *pB* denotes the "bare" result related to the physical (renormalized) pressure as *p*(*T*) = *pB*(*T*) − *pB*(0). The pressure terms *pE*(*T*), *pM*(*T*) and *pG*(*T*) appearing in Equation (61) collect the contributions from the momentum scales *<sup>k</sup>* <sup>∼</sup> *<sup>π</sup>T*, *<sup>k</sup>* <sup>∼</sup> *gT*, and *<sup>k</sup>* <sup>∼</sup> *<sup>g</sup>*2*T*/*π*, respectively. The couplings *g* relevant in the SM are *g* ∈ {*ht*, *g*1, *g*2, *g*3}, where *ht* denotes the Yukawa coupling between the top quark and the Higgs boson, and *g*1, *g*2, *g*<sup>3</sup> are the SM couplings related to U(1)*Y*, SU(2)L and SU*c*(3) gauge groups, respectively. Note that in Equation (61), the thermodynamic limit *V* → ∞ is implied.

Using this technique, the SM calculations of the dimensionless function *p*(*T*)/*T*<sup>4</sup> and of the trace anomaly <sup>Θ</sup>(*T*) (53) up to <sup>O</sup>(*g*5) were performed in Ref. [90]. It was found that similarly to EoS Equation (54), Higgs dynamics induces a peak in heat capacity *c*(*T*) = *d*/*dT* occurring around *T*EW *<sup>c</sup>* ≈ 160 GeV. This leads to a short period of slower temperature change, and correspondingly, a mildly increased abundance of produced particles. However, in general, the largest radiative corrections originate from QCD effects, reducing the energy density by a couple of percent from the free value even at *T* > 160 GeV [90]. It is worth emphasizing that the above-mentioned effects do not exhaust all possible phenomena relevant for dynamics of the early universe. In particular, as already discussed in Section 2.4, at high temperatures, the saturation phenomena become ubiquitous. In this case, the gluon

DoFs are frozen in the classical field (condensate) and do not contribute to the density of the QGP. To shed more light on the microscopic theory in the canonical picture, the dynamical aspects of relativistic gluon and hot meson plasmas, both coupled to the homogeneous condensate, are discussed in detail in Refs. [274–277], respectively.

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−<sup>2</sup> <sup>≤</sup> <sup>≤</sup> <sup>10</sup><sup>16</sup> GeV fm−<sup>3</sup> 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 *<sup>p</sup>*GUT = (0.330 <sup>±</sup> 0.024) valid for 108 <sup>≤</sup> 1016 GeV·fm−<sup>3</sup> (the red triangles in Figure 8) appears only slightly below the ideal gas limit, the hadronic-era EoS *ph* = (0.003 <sup>±</sup> 0.002)+(0.199 <sup>±</sup> 0.002) characterizes the region 1 GeV·fm<sup>−</sup>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* + *c<sup>d</sup>* 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 *<sup>c</sup>* defined implicitly by the equality of the pressures in hadronic and QGP phases *ph*(*c*) = *<sup>p</sup>*SM(*c*) reads: *<sup>c</sup>* (1.2 <sup>±</sup> 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*<sup>1</sup> is always positive, the second pressure *<sup>p</sup>*<sup>2</sup> is negative up to (<sup>7</sup> <sup>−</sup> <sup>13</sup>) GeV·fm<sup>−</sup>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 - (<sup>3</sup> <sup>−</sup> <sup>4</sup>) GeV·fm−3, and for very large energy densities, it falls as <sup>Θ</sup> <sup>∼</sup> *<sup>g</sup>*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 *μ<sup>B</sup>* = 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 *<sup>d</sup>* <sup>=</sup> <sup>−</sup>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 <sup>B</sup>() = <sup>−</sup>(*<sup>a</sup>* <sup>+</sup> *<sup>c</sup>d*), 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 *μ<sup>B</sup>* = 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 <sup>∼</sup><sup>5</sup> <sup>×</sup> <sup>10</sup><sup>7</sup> 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, *<sup>ρ</sup><sup>c</sup>* <sup>10</sup><sup>12</sup> 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<sup>γ</sup>* 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*,*<sup>s</sup>* 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*<sup>0</sup> = 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)−<sup>1</sup> [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) = <sup>−</sup>*p*(*<sup>T</sup>* <sup>=</sup> <sup>0</sup>). 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*σ*˜ <sup>2</sup> <sup>=</sup> *η*2*σ*˜ <sup>2</sup>, 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*<sup>2</sup> *σσ*˜ = 0, etc. Those fluctuations after quantization correspond to physical mesons with masses *m*<sup>2</sup> *<sup>σ</sup>* = 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*<sup>2</sup> *<sup>σ</sup>* and <sup>M</sup><sup>2</sup> <sup>≡</sup> *<sup>v</sup>*<sup>2</sup> <sup>+</sup> *<sup>g</sup>*2*σ*˜ <sup>2</sup> 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 *<sup>d</sup>*2FNE/*dv*<sup>2</sup> <sup>≥</sup> 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*=<sup>0</sup> 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*<sup>0</sup> < *Tc*, when the *σ*-meson fluctuation becomes massless *mσ*(*T*0) = 0, while *v*(*T*0) = 0 is still valid, and *T*<sup>0</sup> 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*<sup>1</sup> 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*<sup>0</sup> = 402 MeV and then grows much faster than the masses of other mesons (such that *m<sup>σ</sup>* > 2*m<sup>π</sup>* 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 *<sup>p</sup>*(*T*), energy density (*T*) and the EoS ((*T*) <sup>−</sup> <sup>3</sup>*p*(*T*) <sup>−</sup> *<sup>A</sup>*)/*T*4, where *<sup>A</sup>* <sup>=</sup> (*<sup>T</sup>* <sup>=</sup> 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 <sup>∝</sup> *<sup>T</sup>*<sup>4</sup> (dashed lines), with the numerical coefficient *<sup>α</sup>* 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*<sup>2</sup> = 0.4, *Tc* = 438 MeV, *T*<sup>0</sup> = 402 MeV, *T*<sup>1</sup> = 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*) <sup>−</sup> <sup>3</sup>*p*(*T*) <sup>−</sup> *<sup>A</sup>*)/*T*<sup>4</sup> as a function of temperature, where *<sup>A</sup>* <sup>≡</sup> (*<sup>T</sup>* <sup>=</sup> <sup>0</sup>) <sup>−</sup> <sup>3</sup>*p*(*<sup>T</sup>* <sup>=</sup> <sup>0</sup>) 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 <sup>0</sup> 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") <sup>0</sup> 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, *<sup>M</sup>*PL <sup>∼</sup> 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 *F<sup>a</sup> μν* 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 <sup>Λ</sup> 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/*a<sup>n</sup>* (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 <sup>Λ</sup> small and positive?) and the "new" Cosmological Constant problem (why is <sup>Λ</sup> 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 Λ−<sup>1</sup> 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, Λ−<sup>1</sup> 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*¯ *<sup>a</sup> <sup>μ</sup>* (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*¯ *<sup>a</sup> <sup>μ</sup>*(*x*) ≡ 0|*A<sup>a</sup> <sup>μ</sup>*(*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 <sup>A</sup>*<sup>a</sup> <sup>μ</sup>* are the rescaled SU(*N*) connections, <sup>F</sup>*<sup>a</sup> μν* are the rescaled field-strength components, and the equality entering the covariant field-strength <sup>F</sup>*<sup>a</sup> μν* in a curved background <sup>∇</sup>*μ*A*<sup>a</sup> <sup>ν</sup>* − ∇*ν*A*<sup>a</sup> <sup>μ</sup>* <sup>=</sup> *∂μ*A*<sup>a</sup> <sup>ν</sup>* <sup>−</sup> *∂ν*A*<sup>a</sup> <sup>μ</sup>* has been accounted for. Furthermore, *g* ≡ det(*gμν*), where *<sup>g</sup>μν* <sup>=</sup> *<sup>a</sup>*(*η*)<sup>2</sup> diag(1, <sup>−</sup>1, <sup>−</sup>1, <sup>−</sup>1) is the FLRW metric as a function of the conformal time *<sup>η</sup>*. 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 *<sup>E</sup>*2 dominate over the averaged CM terms *<sup>B</sup>*2; (ii) vice versa, that is the case of a CM-dominated state with J <sup>&</sup>lt; <sup>0</sup> 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 <sup>J</sup> <sup>&</sup>lt; <sup>0</sup> and, hence, <sup>F</sup> <sup>&</sup>gt; 0. In addition, one refers to a positive *<sup>g</sup>*¯<sup>2</sup> <sup>&</sup>gt; <sup>0</sup> 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*¯<sup>2</sup> 1(*μ*<sup>4</sup> <sup>0</sup>) > 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*¯<sup>2</sup> <sup>1</sup>(<sup>J</sup> ) <sup>&</sup>gt; 0 when −J <sup>&</sup>gt; *<sup>λ</sup>*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 <sup>Z</sup>2-symmetric under simultaneous sign changes of <sup>J</sup> and *<sup>g</sup>*¯2. The invariance of the Lagrangian under this Z<sup>2</sup> 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*¯<sup>2</sup> <sup>∗</sup>) ←→ −*β*(−*g*¯<sup>2</sup> <sup>∗</sup>). 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 <sup>∗</sup> > 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 *μ*<sup>4</sup> <sup>0</sup> <sup>=</sup> |J∗|. Then, for the CM branch, with *<sup>g</sup>*¯<sup>2</sup> <sup>∗</sup> > 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 Z<sup>2</sup> mirror symmetry, the minima of the effective Lagrangian on the two symmetry branches where *g*¯<sup>2</sup> <sup>∗</sup> 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*¯<sup>2</sup> <sup>1</sup>(J∗) < 0.

As is clear from Equation (100) and from the discussion in Ref. [66], the minimum for which J∗ <sup>&</sup>gt; <sup>0</sup> appears in the non-perturbative region defined by <sup>0</sup> <sup>&</sup>lt; J∗ <sup>&</sup>lt; *<sup>λ</sup>*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∗ <sup>&</sup>lt; 0, appears in the perturbative region where |J∗| <sup>&</sup>gt; *<sup>λ</sup>*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 <sup>∼</sup> *<sup>λ</sup>*−<sup>1</sup> <sup>±</sup> <sup>∼</sup> <sup>Λ</sup>−<sup>1</sup> 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 <sup>∼</sup> <sup>Λ</sup>−<sup>1</sup> 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 *<sup>A</sup>*¯*<sup>μ</sup>* and a fluctuating field *<sup>a</sup>μ*: <sup>A</sup>*<sup>μ</sup>* <sup>=</sup> *<sup>A</sup>*¯*<sup>μ</sup>* <sup>+</sup> *<sup>a</sup>μ*. 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 <sup>A</sup>*<sup>a</sup>* <sup>0</sup> = 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*˜ *<sup>a</sup> <sup>k</sup>*(*t*,*x*) ≡ d3*x A*˜ *<sup>a</sup> <sup>k</sup>*(*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 ∼ −<sup>5</sup> <sup>×</sup> <sup>10</sup><sup>9</sup> 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 <sup>∼</sup> <sup>3</sup> <sup>×</sup> <sup>10</sup>−<sup>35</sup> 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 Λ−<sup>1</sup> 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<sup>μ</sup> <sup>μ</sup>* and the total energy density *T*<sup>0</sup> <sup>0</sup>. 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,*<sup>U</sup>* <sup>0</sup> <sup>=</sup> *<sup>T</sup>*0,*<sup>U</sup>* <sup>0</sup> (*U*, *<sup>U</sup>*˙ , *<sup>a</sup>*). After excluding *<sup>T</sup>*0,*<sup>U</sup>* <sup>0</sup> 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 is<sup>6</sup>

$$\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*<sup>∗</sup> ≡ *a*(*t* = *t*0) and the total energy density <sup>0</sup> <sup>≡</sup> *<sup>T</sup>*<sup>0</sup> <sup>0</sup>(*t* = *t*0), respectively. Note that CC 0.

The total energy density, *T*<sup>0</sup> <sup>0</sup>(*t*), and the trace of the EMT, *<sup>T</sup><sup>μ</sup> <sup>μ</sup>*(*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*<sup>0</sup> 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*<sup>∗</sup> <sup>0</sup> = *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,*<sup>U</sup>* <sup>0</sup> (*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,*<sup>U</sup>* <sup>0</sup> , *<sup>T</sup>μ*,*<sup>U</sup> <sup>μ</sup>* . 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*<sup>0</sup> <sup>0</sup>(*t*) and *<sup>T</sup><sup>μ</sup> <sup>μ</sup>*(*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*<sup>0</sup> <sup>0</sup>(*t*) (**left**), the trace of the total QCD EMT *<sup>T</sup><sup>μ</sup> <sup>μ</sup>*(*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*<sup>0</sup> = 0, *U*˙ <sup>0</sup> = (*ξ*ΛQCD)2/ <sup>√</sup>3*e*, *<sup>Q</sup>*<sup>0</sup> <sup>&</sup>gt; 1, *<sup>ξ</sup>* <sup>=</sup> 4.0 <sup>Λ</sup>QCD <sup>=</sup> 332 MeV and *<sup>κ</sup>* <sup>=</sup> <sup>10</sup>−<sup>7</sup> MeV−2. The energy density and the trace are plotted in dimensionless units, rescaled by Λ<sup>4</sup> 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 <sup>Γ</sup>*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) *<sup>k</sup>* , 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) *<sup>k</sup>* <sup>+</sup> *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 Γ*<sup>k</sup>* 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* = <sup>1</sup> 4 *dx F<sup>a</sup> μνF<sup>a</sup> μν*, which here corresponds to the UV limit of the effective theory, the integration constant can be fixed by requiring that Γ*<sup>k</sup>* 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) = *<sup>k</sup>*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 <sup>D</sup>¯ <sup>T</sup> 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, <sup>J</sup> <sup>=</sup> <sup>2</sup>*g*<sup>2</sup> 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 <sup>&</sup>lt; <sup>J</sup> <sup>∗</sup> <sup>&</sup>lt; *<sup>λ</sup>*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 <sup>J</sup> /*λ*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 (<sup>J</sup> <sup>&</sup>lt; *<sup>λ</sup>*4), the two minima in the vicinity of J∗ <sup>&</sup>gt; <sup>0</sup> coincide. Here, <sup>=</sup> 0.01 has been used. (**Middle**) The inverse running coupling *<sup>g</sup>*¯−<sup>2</sup> 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−<sup>5</sup> 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 *<sup>N</sup>* 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<sup>ν</sup>* − ∇*νA<sup>μ</sup>* = *∂μA<sup>ν</sup>* − *∂ν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*<sup>0</sup> = 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*<sup>2</sup> <sup>=</sup> *dt*<sup>2</sup> <sup>−</sup> *<sup>a</sup>*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*) *<sup>i</sup>* 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*) *<sup>i</sup>* , 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/ <sup>√</sup>*N*, hence justifying consistency with the experimental data.

Finally, accounting for an inflationary slow-roll phase, with *B*˙ *<sup>i</sup>* 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, *e<sup>a</sup> <sup>i</sup>* denotes the triads, which recombine into the space metric *hij* = *e<sup>a</sup> i ea <sup>j</sup>* 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* = <sup>1</sup> <sup>3</sup> 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*<sup>∗</sup> 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 *γ<sup>μ</sup>* = *e μ I γI* , with *γ<sup>I</sup>* Dirac matrices, *e μ <sup>I</sup>* 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 *<sup>C</sup><sup>μ</sup>* <sup>=</sup> *<sup>A</sup><sup>μ</sup>* <sup>−</sup> <sup>1</sup> *<sup>m</sup> ∂μθ* and its field strength *<sup>G</sup>μν* <sup>=</sup> *<sup>∂</sup>*[*μCν*] <sup>=</sup> *<sup>∂</sup>*[*μAν*] <sup>=</sup> *<sup>F</sup>μν*, and the fermionic vector current <sup>J</sup> *<sup>μ</sup>* <sup>=</sup> *<sup>q</sup>ψγ*¯ *μψ*. 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*<sup>2</sup> *μναγ* <sup>−</sup> <sup>4</sup>*R*<sup>2</sup> *μν* + *R*<sup>2</sup> and the Weyl tensor *Cμνρσ*.

Considering a spacetime with topology <sup>S</sup><sup>3</sup> <sup>×</sup> <sup>S</sup>1, 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*¯<sup>T</sup> 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 μ <sup>i</sup>* , *i* = *B*, *Q*, *L*, ... denotes the conserved currents corresponding to baryon number *B*, electric charge *Q*, lepton number *L*, etc. Both *Tμν* and *j μ <sup>i</sup>* can be decomposed into time-like

and space-like components using natural projection operators, the local flow four-velocity *<sup>u</sup>μ*, and the second-rank tensor perpendicular to it <sup>Δ</sup>*μν* <sup>=</sup> *<sup>g</sup>μν* <sup>−</sup> *<sup>u</sup>μu<sup>ν</sup>* [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 <sup>=</sup> *<sup>u</sup>μTμνu<sup>ν</sup>* is the energy density, *<sup>p</sup>* <sup>≡</sup> *ps* <sup>+</sup> <sup>Π</sup> <sup>=</sup> <sup>−</sup><sup>1</sup> <sup>3</sup>Δ*μνTμν* is the total (hydrostatic *ps* + bulk <sup>Π</sup>) pressure, *<sup>W</sup><sup>μ</sup>* = <sup>Δ</sup>*<sup>μ</sup> <sup>α</sup>Tαβu<sup>β</sup>* is the energy (or heat) current, *ni* = *uμj μ <sup>i</sup>* is the charge density, *V<sup>μ</sup> <sup>i</sup>* <sup>=</sup> <sup>Δ</sup>*<sup>μ</sup> ν j ν <sup>i</sup>* is the charge current, and *<sup>π</sup>μν* <sup>=</sup> *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 <sup>μ</sup>* <sup>≡</sup> *<sup>j</sup> μ <sup>B</sup>*. 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<sup>μ</sup>* = 0 (*V<sup>μ</sup>* = 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, *<sup>u</sup>μTμν id* <sup>=</sup> *u<sup>ν</sup>* and *<sup>u</sup>μ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<sup>μ</sup>* = *su<sup>μ</sup>* 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 <sup>∇</sup>*<sup>μ</sup>* <sup>=</sup> <sup>Δ</sup>*μν∂ν*. Neglecting the charge current *V<sup>μ</sup>* in Equation (A6), the first-order expansion of *S<sup>μ</sup>* 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* = <sup>1</sup> <sup>3</sup> , 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 *<sup>r</sup>*<sup>⊥</sup> <sup>=</sup> *x*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*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 <sup>∇</sup>*μu<sup>μ</sup>* 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.

Let us provide a few examples of this approach. The 2nd-order theory version of boost-invariant one-dimensional flow, cf. Equation (A13), can be found in Ref. [395]. Due to its length, we do not reproduce it here and refer the interested reader to the original publication. The second example can be found in Ref. [406], where the relaxation time *τπ* proportional to the shear viscosity parameter *η* was used to study the evolution of the universe filled with QGP with nonzero shear viscosity. The authors argue that in general relativity, the following modification of the shear-stress tensor

$$
\pi^{\mu\nu} \to \pi^{\mu\nu} + \pi\_\pi \left[ u^a \pi^{\mu\nu}\_{,\mu} + \frac{4}{3} \pi^{\mu\nu} \nabla\_a u^a \right] = 2\eta \langle \nabla^{\mu} u^{\nu} \rangle, \qquad \pi^{\mu\nu}\_{,a} \equiv \partial\_a \pi^{\mu\nu} + \Gamma^{\mu}\_{a\beta} \pi^{\beta\nu} + \Gamma^{\nu}\_{a\beta} \pi^{\mu\beta}, \tag{A15}
$$

where *πμν* ;*<sup>α</sup>* is a covariant derivative of *<sup>π</sup>μν* and <sup>Γ</sup>*<sup>μ</sup> αβ* are the Christoffel symbols, makes the resulting Navier-Stokes equations causal. Using the FLRW metric and taking into account that the compatibility with the isotropy and homogeneity of the universe demands *πμν* to be diagonal, the solution of Equation (A15) reads [406]

$$
\pi^{00}(t) = \pi^{00}(t\_0) \left[ \frac{a(t\_0)}{a(t)} \right]^4 e^{-\frac{t-t\_0}{\tau\_H}}, \qquad \pi^{ij}(t) = \pi^{ij}(t\_0) \left[ \frac{a(t\_0)}{a(t)} \right]^6 e^{-\frac{t-t\_0}{\tau\_H}} \delta\_{ij}. \tag{A16}$$

In the Friedmann equations, the effect of the traceless viscosity tensor shows up in the modification of the initial energy density (*t*0) and in the behavior of the energy density at times *t t*<sup>0</sup> + *τπ*, which at later times goes over to the standard expression [406]

$$
\varepsilon(t) = \left[\varepsilon(t\_0) + \pi^{00}(t\_0)\right] \left[\frac{a(t\_0)}{a(t)}\right]^4 - \pi^{00}(t\_0) \left[\frac{a(t\_0)}{a(t)}\right]^4 e^{-\frac{t-t\_0}{\overline{\tau}\_T}}.\tag{A17}$$

The third example is provided in Ref. [407] where causal general relativistic viscous fluid theory with the inclusion of all dissipative contributions (shear viscosity *η*, bulk viscosity *ζ*, and heat flow *Wμ*) and the effects from nonzero baryon number are discussed. According to the authors, the applicability of this theory ranges from the modeling of viscous effects in neutron star mergers to low-energy HICs.

#### **Appendix C. YM Equations of Motion in the Effective Action Approach**

This section provides a short summary on the derivation of the YM equations of motion by means of variational methods with respect to the connections <sup>A</sup>*<sup>a</sup> <sup>μ</sup>* when applied to the effective Lagrangian of Equation (84), as was performed in Ref. [66].

When varying the effective action with respect to <sup>A</sup>*<sup>a</sup> <sup>μ</sup>* and *∂ν*A*<sup>a</sup> <sup>μ</sup>*, one arrives at the Euler-Lagrange equations of motion:

$$\frac{\partial \mathcal{L}\_{\mathrm{eff}}}{\partial \mathcal{A}\_{\mu}^{a}} - \nabla\_{\nu} \frac{\partial \mathcal{L}\_{\mathrm{eff}}}{\partial(\partial\_{\nu} \mathcal{A}\_{\mu}^{a})} = 0, \qquad \nabla\_{\nu} \frac{\partial \mathcal{L}\_{\mathrm{eff}}}{\partial(\partial\_{\nu} \mathcal{A}\_{\mu}^{a})} = \frac{1}{\sqrt{-g}} \partial\_{\mathbb{V}} \left[ \sqrt{-g} \frac{\partial \mathcal{L}\_{\mathrm{eff}}}{\partial(\partial\_{\nu} \mathcal{A}\_{\mu}^{a})} \right]. \tag{A18}$$

It is straightforward to compute the derivatives of the effective Lagrangian:

$$\frac{\partial \mathcal{L}\_{\text{eff}}}{\partial \mathcal{A}^{a}\_{\mu}} = \frac{1}{4\xi^{2}} \left[ \frac{\partial \mathcal{J}}{\partial \mathcal{A}^{a}\_{\mu}} - \frac{\mathcal{J}}{\mathcal{J}^{2}} \frac{\partial \mathcal{J}^{2}}{\partial \mathcal{A}^{a}\_{\mu}} \right], \qquad \frac{\partial \mathcal{J}}{\partial \mathcal{A}^{a}\_{\mu}} = \frac{4f^{abc} \mathcal{F}^{b} \mathcal{W}\_{\text{\mu}}}{\sqrt{-g}}, \qquad \frac{\mathcal{J}}{\mathcal{J}^{2}} \frac{\partial \mathcal{g}^{2}}{\partial \mathcal{A}^{a}\_{\mu}} = \frac{\mathcal{J}}{\mathcal{J}^{2}} \frac{\partial \mathcal{g}^{2}}{\partial \mathcal{J}} \frac{\partial \mathcal{J}}{\partial \mathcal{A}^{a}\_{\mu}}, \tag{A19}$$

$$\implies \quad \frac{\partial \mathcal{L}\_{\rm eff}}{\partial \mathcal{A}^{a}\_{\mu}} = \frac{1}{\delta^{2}} \frac{f^{\rm abc} \mathcal{F}^{b} \mu^{\mu} \mathcal{A}^{c}\_{v}}{\sqrt{-g}} \left[1 - \frac{\beta}{2}\right],\tag{A20}$$

$$\frac{\partial \mathcal{L}\_{\text{eff}}}{\partial(\partial\_{\nu} \mathcal{A}\_{\mu}^{a})} = \frac{4 \mathcal{F}^{u \mu \nu}}{\sqrt{-g}}.\tag{A21}$$

$$\Rightarrow \quad \nabla\_{\boldsymbol{V}} \frac{\partial \mathcal{L}\_{\mathrm{eff}}}{\partial(\partial\_{\boldsymbol{V}} \mathcal{A}\_{\mu}^{a})} = \frac{1}{\sqrt{-\mathcal{g}}} \partial\_{\boldsymbol{V}} \left( \sqrt{-\mathcal{g}} \frac{1}{\overline{\mathcal{g}}^{2}} \frac{\mathcal{F}^{\mu \nu \boldsymbol{\nu}}}{\sqrt{-\mathcal{g}}} \left[ 1 - \frac{\mathcal{\boldsymbol{\mathcal{G}}}}{2} \right] \right), \tag{A22}$$

where the RG equation for the exact *β*-function that conveniently recasts Equation (86) as

$$\frac{\mathcal{J}}{\mathcal{J}^2} \frac{\partial \mathcal{g}^2}{\partial \mathcal{J}} \equiv \frac{\mathcal{B}}{2} = \frac{\operatorname{d} \ln |\mathcal{J}^2|}{\operatorname{d} \ln |\mathcal{J}| / \mu\_0^4}, \quad \mathcal{\beta} = \beta(\mathcal{g}^2), \tag{A23}$$

has been inserted. If *β* is known to all-loop order, the running of the coupling and hence the solutions to the equations of motion may be found to all-loop order accuracy. In the expression above, an arbitrary dimensionful renormalisation parameter *μ*<sup>0</sup> has been explicitly introduced as a reference scale. The natural boundary condition is *g*¯(J ) → *g*YM when |J | → *<sup>μ</sup>*<sup>4</sup> 0.

Inserting Equations (A20) and (A22) into Equation (A18), the resulting equations of motion is

$$\frac{1}{\hat{\mathcal{S}}^2} \frac{f^{abc} \mathcal{F}^{b \, \mu \nu} \mathcal{A}\_v^c}{\sqrt{-\mathcal{g}}} \left[1 - \frac{\beta}{2}\right] - \frac{1}{\sqrt{-\mathcal{g}}} \partial\_\nu \left(\sqrt{-\mathcal{g}} \frac{1}{\overline{\mathcal{g}}^2} \frac{\mathcal{F}^{a \, \mu \nu}}{\sqrt{-\mathcal{g}}} \left[1 - \frac{\beta}{2}\right]\right) = 0. \tag{A24}$$

This can be rewritten on the operator form

$$\mathcal{D}\_{\nu}^{ab} \left[ \frac{\mathcal{F}^b \mu \nu}{\bar{\mathcal{S}}^2 \sqrt{-g}} \left( 1 - \frac{\beta}{2} \right) \right] = 0 \,,\tag{A25}$$

where the differential operator <sup>D</sup><sup>ˆ</sup> is given by

$$\mathcal{D}^{ab}\_{\nu} \equiv \delta^{ab} \frac{\partial\_{\nu} \sqrt{-\mathcal{g}}}{\sqrt{-\mathcal{g}}} - f^{abc} \mathcal{A}^c\_{\nu}. \tag{A26}$$

The action of this differential operator on a function *h*(*x*) is defined as follows:

$$\left[\left[\hat{\mathcal{D}}h(\mathbf{x})\right]\_{\nu}^{ab}\equiv\delta^{ab}\frac{\partial\_{\nu}\left[\sqrt{-\mathcal{G}}h(\mathbf{x})\right]}{\sqrt{-\mathcal{G}}}-f^{abc}\mathcal{A}\_{\nu}^{c}h(\mathbf{x})\,. \tag{A27}$$

#### **Appendix D. One-Loop Effective YM Lagrangian**

Let us briefly discuss the effective YM theory at the one-loop order. The usefulness of studying the one-loop case is further motivated by a comparison of the one-loop and all-loop order expansion in Section 4.5. The standard one-loop SU(*N*) *β*-function reads (see, e.g., Ref. [408])

$$\beta\_1 \equiv -B\_1 \mathfrak{g}\_1^2, \quad B\_1 = \frac{bN}{48\pi^2}, \quad b = 11,\tag{A28}$$

and the corresponding solution of the RG equation (Equation (A23)) is given by

$$\bar{\varrho}^2(\mathcal{I}) = \frac{\mathcal{J}\_1^2(\mu\_0^4)}{1 + \frac{B\_1}{2}\mathcal{J}\_1^2(\mu\_0^4)\ln\left(|\mathcal{I}|/\mu\_0^4\right)}.\tag{A29}$$

Substituting this expression into the effective all-order Lagrangian in Equation (84), we obtain

$$\mathcal{L}\_{\rm eff}^{(1)} = \frac{\mathcal{I}}{4\bar{\chi}\_1^2(\mu\_0^4)} \left[ 1 + \frac{B\_1}{2} \bar{\varrho}\_1^2(\mu\_0^4) \ln\left(\frac{|\mathcal{I}|}{\mu\_0^4}\right) \right]. \tag{A30}$$

Making trivial substitutions,

$$
\mathcal{J} \to -4\mathfrak{g}\_1^2(\mu\_0^4)\mathcal{F}, \quad \mu\_0^4 \to 2\varepsilon\mu^4, \quad \mathfrak{g}\_1^2(\mu\_0^4) \to \mathfrak{g}\_{\text{YM}}^2. \tag{A31}
$$

one arrives at another form of the one-loop effective Lagrangian frequently used in the literature (e.g., Ref. [345] and references therein),

$$\mathcal{L}\_{\rm eff}^{(1)} = -\mathcal{F} - \frac{bN}{96\pi^2} g\_{YM}^2 \mathcal{F} \left[ \ln \left( \frac{2|g\_{YM}^2 \mathcal{F}|}{\mu^4} \right) - 1 \right]. \tag{A32}$$

The compact form of the all-order effective Lagrangian used earlier in Equation (84) straightforwardly produces the standard representation of the one-loop effective Lagrangian upon the redefinitions of Equation (A31), which is reassuring. The usual covariant renormalization condition on the effective Lagrangian [345]

$$\frac{\partial \mathcal{L}\_{\rm eff}}{\partial \mathcal{F}}\Big|\_{t=0} = -1 \, , \qquad t \equiv \frac{1}{2} \ln \left( \frac{2|\mathcal{G}\_{\rm YM}^2 \mathcal{F}|}{\mu^4} \right) \, , \tag{A33}$$

is apparently satisfied for Equation (A32). Indeed,

$$\frac{\partial \mathcal{L}^{(1)}\_{\text{eff}}}{\partial \mathcal{F}} = -1 - \frac{bN}{96\pi^2} g\_{\text{YM}}^2 \ln\left(\frac{2|g\_{\text{YM}}^2 \mathcal{F}|}{\mu^4}\right) \to -1 \qquad \text{for} \qquad \ln\left(\frac{2|g\_{\text{YM}}^2 \mathcal{F}|}{\mu^4}\right) \to 0. \tag{A34}$$

This condition has been originally employed in Refs. [326,344] to derive the generic form of the one-loop effective Lagrangian in Equation (A32) (see Ref. [345] and references therein, for a more elaborate review). In a compact notation of Equation (84), the latter condition reads

$$\frac{\partial \mathcal{L}\_{\text{eff}}}{\partial \mathcal{J}}\Big|\_{t=0} = \frac{1}{4\mathcal{g}\_1^2(\mu\_0^4)}, \qquad t \equiv \frac{1}{2}\ln\left(\frac{e|\mathcal{J}|}{\mu\_0^4}\right). \tag{A35}$$

**Notes**

<sup>1</sup> The light cone four-vectors are related to Minkowski four-vectors in a standard way *<sup>k</sup>* = (*k*0, *<sup>k</sup>*⊥, *kz*)=[*k*−, *<sup>k</sup>⊥*, *<sup>k</sup>*+] where *<sup>k</sup>*<sup>±</sup> <sup>=</sup> *<sup>k</sup>*<sup>0</sup> <sup>±</sup> *<sup>k</sup>z*. The Minkowski dot product in light-cone coordinates is *<sup>k</sup>* · *<sup>p</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> (*k*<sup>+</sup> *<sup>p</sup>*<sup>−</sup> <sup>+</sup> *<sup>k</sup>*<sup>−</sup> *<sup>p</sup>*+) <sup>−</sup> *<sup>k</sup><sup>⊥</sup>* · *<sup>p</sup>⊥*.


$$f - h(t)(f)^2 + Ah(t) = 0,$$

with ˜ *h*(*t*) = <sup>1</sup> 2 *g*(*t*) − 3 and *A* = *κ*CC <sup>3</sup> . The introduction of *<sup>m</sup>*(*t*) <sup>≡</sup> ˙ *f* , results in a first-order equation that may be solved. It is explicitly

$$m - \ddot{h}(t)m^2(t) + A\ddot{h}(t) = 0\dots$$

The scale factor is therefore found in terms of the integral of the solution for *m*(*t*) as

$$a(t) = a^\* \exp\left[\int\_{t\_0}^t \mathbf{d}t \, m(t)\right].$$

.

<sup>7</sup> For the case of the simple background considered in Section 4.4, <sup>A</sup>*<sup>μ</sup>* <sup>=</sup> *<sup>A</sup>*¯*<sup>μ</sup>* <sup>+</sup> *<sup>a</sup><sup>μ</sup>* , then <sup>Γ</sup>(*n*,*m*) *k A*¯, *a* = *<sup>δ</sup><sup>n</sup>* (*δA*¯)*<sup>n</sup> δm* (*δa*)*<sup>m</sup>* Γ*<sup>k</sup> A*¯, *a* 

#### **References**

