Andrea Addazi ^1,
, Antonino Marcianò ^1,
, Roman Pasechnik ^2,
, 1.Department of Physics & Center for Field Theory and Particle Physics, Fudan University, Shanghai 200433, China
2.Department of Astronomy and Theoretical Physics, Lund University, Lund SE-223 62, Sweden
Received Date:2019-01-27
Available Online:2019-06-01
Abstract:We propose a novel mechanism for the production of gravitational waves in the early Universe that originates from the relaxation processes induced by the QCD phase transition. While the energy density of the quark-gluon mean-field is monotonously decaying in real time, its pressure undergoes a series of violent oscillations at the characteristic QCD time scales that generate a primordial multi-peaked gravitational waves signal in the radio frequencies’ domain. The signal is an echo of the QCD phase transition that is accessible by planned measurements at the FAST and SKA telescopes.

HTML

--> --> -->

The intriguing possibility that prompt phase transitions in the early Universe might have imprinted signatures in the background of gravitational radiation will be testable through the next generation of gravitational interferometers. The idea was firstly suggested in Refs. [1-5]. Important developments on understanding the primordial gravitational waves (GW) production in the early Universe were made, in particular, in Refs. [6,7].
Dynamical effects of quarks confinement in baryons and mesons are related to the dimensional-scale transmutation as much as first-order phase transition (FOPT) phenomena, which are characterized by a dynamically generated energy scale $ \Lambda_{\rm QCD}\simeq 200\, {\rm MeV} $ [8]. This suggested the possibility that the Quantum Chromodynamics (QCD) phase transition may generate a GW signal in the hot Early Universe, at a temperature of $ T\simeq \Lambda_{\rm QCD}\simeq 200\, {\rm MeV} $. That a FOPT related to strong interactions may emit GWs was initially, although qualitatively, proposed by Witten [1], and then quantitatively re-elaborated in Refs. [9,10].
The GW signal associated to the QCD phase transition (QCDPT) cannot be detected in GW terrestrial interferometers, such as LIGO/VIRGO [11] and KAGRA [12], and cannot be either measured in future space experiments, such as LISA [13], U-DECIGO [14], BBO [15], TAIJI [16] and TianQin [17] projects. The GWs frequency range of a QCDPT is around $ 10^{-8}\div 10^{-9}\, {\rm Hz} $, which is 5-6 digits lower than the one provided by space experiments, and 9-10 digits away from LIGO/VIRGO/KAGRA limits [10]. Furthermore, QCDPT does not leave any smoking-gun imprinting in the Cosmic Microwave Background as being sensitive to very low frequency modes (5-6 digits less [10]).
Nonetheless, a nHz phase transition such as a QCDPT can be detected, with high precision, from radio astronomical observation of pulsar timing: the GW backgrounds propagating through pulsar systems alter the radio signal, leaving an imprint that is in principle observable^①. This opens a pathway towards the exciting possibility of testing fundamental particle physics with current and future radio astronomy experiments, including FAST [20] and SKA [21]. Within previous QCDPT analyses, the role of possible relaxation phenomena and gravitational back-reactions were completely neglected. Note, after the QCDPT, the prompt and violent relaxation effects around the QCD vacuum energy state are expected, which retains a broad analogy with the reheating mechanism in inflationary models.
In this letter, we study in detail the possible effects of the gluon condensate relaxation phenomena. We analyze the non-linear field equations for the gluonic condensate, coupled to the Einstein equations, in a Friedmann-Lema?tre-Robertson-Walker (FLRW) cosmological background. Our formalism is based upon the effective Einstein–Yang-Mills equations of motion for real-time evolution of the homogeneous gluon condensate in expanding Universe initially formulated in Ref. [22]. During the relaxation phase, a surprising non-equilibrium phenomenon arises: the gluonic condensate field violently oscillates during the relaxation phase, inducing fast oscillations of the energy-momentum tensor trace for a transient time of $ \tau \simeq 10\div 20 \Lambda_{\rm QCD}^{-1} $. The oscillating solution is a classical non-perturbative solution of the Yang-Mills field equations coupled to the Einstein field equations. The emergence of spikes, localized in a characteristic QCD time lapse $ \Delta t\simeq \Lambda_{\rm QCD}^{-1} $, and extended in the space dimensions, reveals the presence of a ordered pattern of space-like soliton/domain-walls solutions. We dub these new solutions chronons. At large cosmological times $ t > 20 \Lambda_{\rm QCD}^{-1} $, the spikes' periodicity disappears, and the energy density approaches the QCD vacuum energy minimum. The time-ordered classical solution that we found is a time crystal, i.e. a periodic classical solution spontaneously breaking time invariance down to a discrete time shift symmetry $ T_{n}:t \rightarrow t+n \Lambda_{\rm QCD}^{-1} $, n denoting a natural number. The concept of time crystal has been first proposed by Wilczek in Refs. [23,24] within the context of superconductors and superfluids physics^①. For a review of time crystals, see e.g. Ref. [27]. The experimental discovery of time crystals was achieved in Ref. [28]. The spontaneous symmetry breaking of T-invariance from the localization of chronons is associated to the appearance of Nambu-Goldstone bosons, as time-like moduli excitations over the classical background.
During the relaxation stage, a new characteristic feature is the produced GW signal. While the energy-density part of the energy momentum tensor does not exhibit so violent transitions, the condensate pressure provides the main contribution to the energy-momentum tensor trace variation. These pressure kinks inject kinetic energy into the primordial plasma, inducing turbulence and sound/shock waves in the plasma very efficiently. In analogy with the case of a bubble propagating in the plasma, the gravitational radiation is emitted from magnetohydrodynamical (MHD) turbulence and sound waves. From our numerical simulations, which we compare with semi-analytical estimates, we show that such gravitational background signal can be tested in future radio observatories form pulsar timing effects. The spectrum that is predicted not only lies within the SKA sensitivity, but it further displays very peculiar features of the shape form that cannot be reproduced in any other known mechanism. In other words, time crystallization of the QCD medium during the relaxation phase can be tested in close future, which implies a radical reconsideration of our picture of QCD confinement itself from the prospective of dynamical cosmological evolution.

A standard static domain-wall can be easily obtained from a scalar field theory that is $ Z_{2} $ invariant. With a simple sombrero-like Higgs potential, $ Z_{2} $ can be spontaneously broken when the scalar field rolls down to one of the two possible minima $ \phi_{vac} = \pm v $. These internal field configurations can be localized in the space direction z as kink profiles. The kink profile interpolates the two minima, namely, $ \phi(z = -\infty) = -v $ and $ \phi(z = \infty) = v $. A domain-wall configuration, as a xy-plane orthogonal to the z-direction, is achieved through the kink profile transition region, and its characteristic thickness in z-direction is directly related to the kink shape. For a $ \lambda \phi^{4} $ theory with a sombrero potential, one can find a simple analytic kink solution, specified by $ \phi(z) = v\tanh\Big[\displaystyle\frac{\lambda v}{\sqrt{2}}(z-z_{0})\Big] $, with $ z_{0} $ being the kink center.
As is well-known, for the standard domain-walls the translational invariance is spontaneously broken, being the barrier localized in a $ z_{0} $ point. This corresponds to the appearance of a Nambu-Goldstone modulus boson $ z_{0}(t,x,y) $, localized on the surface of the domain wall, as a low energy-excitation of its surface in the z-direction.
Intriguingly and exotically, one may consider a kink profile that, despite of being localized in a space direction, is also localized in time. A new domain wall extended in three spatial dimension but localized in a time lapse, which we dub as a chronon, may correspond to this solution. By just replacing the z-coordinate with the time variable, one can consider a kink solution, such as $ \phi(t) = v\tanh\Big[\displaystyle\frac{\lambda v}{\sqrt{2}}(t-t_{0})\Big] $, centered at a time instant $ t = t_{0} $ and interpolating the two vacuum states in the asymptotic time limits $ \phi(t = -\infty) = -v $ and $ \phi(t = \infty) = v $. This solution is associated to a spontaneous symmetry breaking of the time invariance and to the appearance of a Nambu-Goldstone boson localized on the xyz-surface $ t_{0}(x,y,z) $.
In the case of the gluon condensate field equation coupled to gravity, in a FLRW cosmological background one can decompose the gluonic field in a classical background field $ U(t) $ plus a non-homogeneous part—see the Appendix for more technical details. Let us consider the limit of a static FLRW space time ($ a = {\rm const} $). A branch of solutions for the U-field satisfies the equation

${U'^2} - \frac{1}{4}{U^4} = {\rm{const}}{\mkern 1mu} ,$

(1)

where U′ is the field derivative with respect to the Cartesian coordinate time, which we denote here with $ x_0\equiv \eta $. A branch of solutions of these equations, obtained by $ U^{2}\rightarrow U^{2}-U_{0}^{2} $ energy density vacuum shift, corresponds to kink (anti-kink) profiles

$U(\eta ) \simeq \frac{v}{{\sqrt 2 }}\tanh \left[\frac{v}{{\sqrt 2 }}(\eta - {\eta _0})\right]{\mkern 1mu} ,$

(2)

where $ v\simeq \Lambda_{\rm QCD} $. A space-like domain wall corresponds to a kink profile of this type. Time-translational invariance is spontaneously broken, and a $ \eta_{0}(x,y,z) $ moduli field arises, with U acquiring the form, $ U(\eta-\eta_{0}(x,y,z)) $. The coordinates x,y,z are the domain-wall's worldsheet coordinates. The effective corresponding action reads as

$\begin{split}S =& \int {{\rm d}^4}x{\mkern 1mu} \frac{1}{2}\left[{\left(\frac{{\partial \phi }}{{\partial {\eta _0}}}\frac{{\partial {\eta _0}}}{{\partial {x^a}}}\right)^2} - V(U)\right]\\ =& {\rm{const}} + \frac{{{T_W}}}{2}\int {{\rm d}^3}x{\mkern 1mu} {\left(\frac{{\partial {\eta _0}({x^a})}}{{\partial {x^a}}}\right)^2}{\mkern 1mu} .\end{split}$

(3)

This shows that the moduli field is massless, according to the Nambu-Goldstone theorem.
When the gravitational dynamics is taken into account, and the scale factor time-dependence is considered, a more complicated time pattern for the space-like domain walls arises— see Appendix for technical details. In this latter case, time-translation is not only broken down to a $ Z_{2} $ symmetry involving the $ T_{+} $ and $ T_{-} $ branches, but more interestingly is spontaneously broken down to a discrete time-translation symmetry $ T_{n}:\eta\rightarrow \eta+n\Lambda_{\rm QCD}^{-1} $. In other words, the system behaves as a time crystal.

The general coupled field equations of the gluon field with gravity reads

$\begin{split}\frac{1}{\varkappa}\left( {R_\mu ^\nu - \frac{1}{2}\delta _\mu ^\nu R} \right) =& \frac{b}{{32{\pi ^2}}}\frac{1}{{\sqrt { - g} }}\left[\left( - {\cal F}_{\mu \lambda }^a{\cal F}_a^{\nu \lambda } +{\mkern 1mu} \frac{1}{4}\delta _\mu ^\nu {\cal F}_{\sigma \lambda }^a{\cal F}_a^{\sigma \lambda }\right)\right. \\&\times\left.\ln \frac{{e|{\cal F}_{\alpha \beta }^a{\cal F}_a^{\alpha \beta }|}}{{\sqrt { - g} {\mkern 1mu} {\lambda ^4}}} - \frac{1}{4}\delta _\mu ^\nu {\mkern 1mu} {\cal F}_{\sigma \lambda }^a{\cal F}_a^{\sigma \lambda }\right]{\mkern 1mu} ,\end{split}$

(4)

$\begin{split}\left( {\frac{{{\delta ^{ab}}}}{{\sqrt { - g} }}{\partial _\nu }\sqrt { - g} - {f^{abc}}{\cal A}_\nu ^c} \right)\left( {\frac{{{\cal F}_b^{\mu \nu }}}{{\sqrt { - g} }}{\mkern 1mu} \ln \frac{{e|{\cal F}_{\alpha \beta }^a{\cal F}_a^{\alpha \beta }|}}{{\sqrt { - g} {\mkern 1mu} {\lambda ^4}}}} \right) = 0{\mkern 1mu} ,\end{split}$

where $ \lambda\equiv \xi \Lambda_{\rm QCD} $ is related to the QCD scale parameter $ \Lambda_{\rm QCD} $ by an arbitrary scaling constant $ \xi $, and e is the base of the natural logarithm. In the FLRW background, here cast in conformal coordinates defined by $ {\rm d}s^2 = a^2(\eta)({\rm d}\eta^2 + {\rm d}\vec{x}^{\,2}) $, the dynamical system is simplified to

$\begin{split}&\frac{6}{\varkappa} \frac{{a''}}{{{a^3}}} = T_\mu ^{\mu ,{\rm{U}}}{\mkern 1mu} ,\\&T_\mu ^{\mu ,{\rm{U}}} = \frac{{3b}}{{16{\pi ^2}{a^4}}}\left[{(U')^2} - \frac{1}{4}{U^4}\right]{\mkern 1mu} ,\\&\frac{\partial }{{\partial \eta }}\left(U'{\mkern 1mu} \ln \frac{{6e|{{(U')}^2} - \frac{1}{4}{U^4}|}}{{{a^4}{\lambda ^4}}}\right)\\&\quad + {\mkern 1mu} \frac{1}{2}{U^3}{\mkern 1mu} \ln \frac{{6e|{{(U')}^2} - \frac{1}{4}{U^4}|}}{{{a^4}{\lambda ^4}}} = 0{\mkern 1mu} .\end{split}$

(5)

In order to account for thermal bath effects, we have to consider the thermal loop correction to the classical equations. The leading order corrections are proportional to $ T^{2}U^{2} $, where T is the temperature of the early Universe plasma.
Solving the dynamical system of equations specified above, at T=0, we find a $ U(\eta) $ profile characterized by a time-ordered pattern of spikes—see the Appendix for more details. Turning on thermal corrections, the only relevant $ O(1) $ corrections arise in correspondence of the first spike, close to the QCD phase transition scale $ \Lambda_{\rm QCD} $. After $ 10\,\Lambda_{\rm QCD}^{-1} $, thermal correction will reduce to a small effect. The trace of energy-momentum tensor follows the spike series solution, with violent oscillations while having relaxation. Relaxation is induced by the linear dissipative terms of the gluonic condensate in the field equations of motion, while spikes can be understood as a back-reaction effect of the gluonic condensate with the gravitational background. The energy-density part $ \rho_{U} $ has only a suppressed modulation over the relaxation decay profile in time. Since the gravitational waves emission is related to the time variation of the energy-density, it turns out that the GW spectrum is actually suppressed. Most of the trace tensor variation is provided by the pressure component $ p_{U} $. The pressure kinks can be efficiently transmitted to the primordial plasma, since the gluonic condensate is strongly interacting with it. Consequently, the pressure kinks pattern chaotized (thermalised) the early Universe plasma, inducing turbulence and shock sound waves in it. As it happens in the standard picture of FOPTs, turbulence wiggles and sound waves efficiently produce the gravitational radiation. The pressure distribution of one kink in the plasma thermalizes very efficiently, within a standard deviation of the order of the QCD scale. Numerical results are shown in Fig. 1 and Fig. 2, and put in comparison with the current FAST experiment sensitivity curves and future SKA predicted bounds. The model efficiency factors enter crucially the pressure transfer mechanism from the classical condensate to the Universe plasma. Although clearly affected by several uncertainties, one can show several GW profiles in the reasonable efficiency range of $ 10^{-2}\div 10^{-3} $. We remark that an efficiency lower than $ 10^{-3} $ seems to be nearly impossible, since the gluonic condensate cannot be so weakly coupled to the plasma. This observation carries important phenomenological consequences. Indeed, we can show that, for an efficiency higher than 3%, FAST will definitely probe this model. On the other hand, future FAST data, one can reach the 0.1% efficiency scale.
Figure1. (color online) The gravitational waves spectrum is displayed for different efficiency factors, in comparison with the FAST sensitivity curve [20]. The efficiency factor is considered to be $\kappa = 0.03\div 0.1$.

Figure2. (color online) The gravitational waves spectrum is displayed for different efficiency factors, in comparison with FAST sensitivity curve [21]. The efficiency factor is considered to vary within $\kappa = 10^{-3}\div 3\times 10^{-3}$.

We provide below simple semi-analytic estimates, which are nonetheless in agreement with our numerical analysis.
The red-shift due to the expansion of the gravitational background must be taken into account while comparing the GW signals ab origine with measurements at present time. At this purpose, we recall that the ratio between the scale factor of the Universe today $ R_{0} $ and the scale factor of the Universe during the GW production is expressed by

$\frac{{{R_*}}}{{{R_0}}} = 8.0 \times {10^{ - 14}}{\left( {\frac{{100}}{{{g_*}}}} \right)^{\frac{1}{3}}}\left( {\frac{{1{\mkern 1mu} {\rm{GeV}}}}{{{T_*}}}} \right){\mkern 1mu} ,$

(6)

$ g_* $ denoting the effective number of degrees of freedom. Assuming that the Universe expanded adiabatically implies that the entropy $ S\sim R^{3}\, T^{3} $ remained constant. The characteristic frequency of the GW signal today, denoted as $ f_{0} $, is related to the one on the GW emission time $ f_{*} $ by

${f_0} = {f_*}\left( {\frac{{{R_*}}}{{{R_0}}}} \right) = 1.65 \times {10^{ - 7}}{\rm{Hz}}\left( {\frac{{{f_*}}}{{{H_*}}}} \right)\left( {\frac{{{T_*}}}{{1{\mkern 1mu} {\rm{GeV}}}}} \right){\left( {\frac{{{g_*}}}{{100}}} \right)^{\frac{1}{6}}}.$

(7)

The order of magnitude of the GW energy density today, denoted as $ \Omega_{\rm GW} $, is related to the one during the emission time, namely $ \Omega_{\rm GW *} $, through

$\begin{split}{\Omega _{{\rm{GW}}}} =& {\Omega _{{\rm{GW*}}}}{\left( {\frac{{{R_*}}}{{{R_0}}}} \right)^4}{\left( {\frac{{{H_*}}}{{{H_0}}}} \right)^2}\\ =& 1.67 \times {10^{ - 5}}{h^{ - 2}}{\left( {\frac{{100}}{{{g_*}}}} \right)^{\frac{1}{3}}}{\Omega _{{\rm{GW*}}}}{\mkern 1mu} ,\end{split}$

(8)

where h is the current value of the Hubble parameter in units of $ 100\, {\rm km}\,/ ({\rm sec}\, {\rm Mpc}) $, and

${H_*} = \frac{{8\pi G{\rho _{{\rm{rad}}}}}}{3} = \frac{{8{\pi ^3}{g_*}T_*^4}}{{90M_{Pl}^2}}{\mkern 1mu} ,$

(9)

which is the Hubble contribution in the radiation dominated epoch.
The sound and turbulence spectrum induced by the spiky pressure kinks is in general very complicated. However, it will display a characteristic series of peaks, related to the pressure peaks. The magnitude of these GW peaks can be estimated very easily, thanks to semi-analytical estimates. The turbulence GW peaks are described by

${h^2}{\Omega _{{\rm{turb}}}} = 3.35 \times {10^{ - 4}}\sum\limits_i^{{N_{{\rm{eff}}}}} \left( {\frac{{{H_{*,i}}}}{{{\beta _i}}}} \right){\left( {\frac{{{\kappa _{\rm turb}}\alpha }}{{1 + \alpha }}} \right)^{\frac{3}{2}}}{\left( {\frac{{100}}{{{g_{*,i}}}}} \right)^{\frac{1}{3}}}v,$

(10)

${f_{{\rm{turb}}}} = 2.7 \times {10^{ - 2}}{\mkern 1mu} {\rm{mHz}}\frac{1}{v}\left( {\frac{{{\beta _i}}}{{{H_{*,i}}}}} \right)\left( {\frac{{{E_{*,i}}}}{{{\rm{GeV}}}}} \right){\left( {\frac{{{g_{*,i}}}}{{100}}} \right)^{\frac{1}{6}}},$

(11)

where the sum is over the number of peaks that contribute significantly to the GW spectrum, and

$\alpha = \frac{{{\rho _U}}}{{{\rho _{{\rm{rad}}}}}}$

(12)

is the ratio between the energy density U and the radiation energy density.
The sound waves spectrum is characterized by the expressions

${h^2}{\Omega _{{\rm{sound}}}} = 2.65 \times {10^{ - 6}}\sum\limits_{i = 1}^{{N_{{\rm{eff}}}}} \left( {\frac{{{H_{*,i}}}}{{{\beta _i}}}} \right){\left( {\frac{{{\kappa _v}\alpha }}{{1 + \alpha }}} \right)^2}{\left( {\frac{{100}}{{{g_{*,i}}}}} \right)^{\frac{1}{3}}}{v_i},$

(13)

${f_{{\rm{sw}}}} = 1.9 \times {10^{ - 2}}{\rm{mHz}}\sum\limits_{i = 1}^{{N_{{\rm{eff}}}}} {\frac{1}{{{v_i}}}} \left( {\frac{{{\beta _i}}}{{{H_{*,i}}}}} \right)\left( {\frac{{{T_{*,i}}}}{{100{\mkern 1mu} {\rm{GeV}}}}} \right){\left( {\frac{{{g_{*,i}}}}{{100}}} \right)^{\frac{1}{6}}}.$

(14)

We estimated the U peak rapidity as $ v_{i} = \Delta (T_{\mu}^{\mu}) \beta_{i} $, which is close to 1 (fast), and the inverse time scale of the peaks as $ \beta_{i} \simeq \Lambda $ ($ \Lambda $ being the confinement scale), while $ \kappa_{\rm v, turb} $ are the efficiency factors of the energy transfer from the condensate peaks to the early Universe plasma. We find that: i) the energy-scale of the first peaks is around $ 200\, {\rm MeV} $; ii) $ \beta_{i} \simeq H_{*,i} $, since $ H_{*} $ is related to the energy density and pressure of the Universe and follows the oscillations of the condensate; iii) $ \alpha \simeq 1 $, if the condensate dominates in the radiation epoch.
Assuming the efficiency factors $ \kappa_{\rm turb,v}\simeq $ 0.1%÷1%—in analogy to typical response factors in the plasma during FOPTs—from the estimates specified above, we obtain an energy frequency within the frequency range $10^{-9}\, \div 10^{-8}\, {\rm Hz}$ for both turbulence and sound waves. This implies that the GW signal here predicted lies in the radio-astronomy pulsar timing scale, while the energy-density of the GW signal for both the contributions is around $10^{-8}\div 10^{-12}$. Thus, SKA must be definitely able probe to this GW signal.

We have shown that the relaxation dynamics of the gluon condensate close to the QCD phase transition behaves as a time-crystal within a time range of $ 1\div 20\Lambda_{QCD^{-1}} $. This is an effect obtained as a non-equilibrium solution from the classical gluon equations of the gauge fields coupled with gravity. Time translation is spontaneously broken in a time translation discrete symmetry. We have shown that this proposal to model gluonic condensates in the early Universe cosmology is already testable. The model we developed predicts a gravitational radiation background that can alter the pulsar timing system, and thus can be tested in radio-astronomy experiments. Specifically, predictions of the model lay in sensitivity curves of FAST and SKA. A clear understanding of QCD confinement dynamics remains the most challenging problem of the Standard Model of particles physics. Possible new information about the confinement dynamics in the early Universe from radio-astronomy opens a new pathway towards the frontiers of strong interaction physics.

We discuss in this section how a gauge-invariant description of spatially homogeneous isotropic Yang-Mills (YM) condensates that depend only on time can be obtained. For this purpose, it is most useful to work in the ghost-free temporal (Hamilton or Weyl) gauge, fixed ${\rm{by}} $ the condition

$\tag{A1}A_0^a = 0{\mkern 1mu} ,$

which is the basis of the Hamiltonian formulation. In this gauge, the asymptotic states of the S-matrix automatically contain transverse modes only, which enables to formulate the YM theory in the Heisenberg representation, consistently beyond the phase transitions—for more details, see e.g. Refs. [29,30]. Under the condition (A1), the Gauss-laws, which enforce the symmetry under the residual time-independent gauge transformations, have to be implemented additionally to the Hamiltonian equations of motion. These equations can be brought into an unconstrained form by imposing a further time-independent gauge. While in the perturbative phase the non-Abelian Coulomb gauge is a useful choice, in the case of spatially homogeneous condensates the symmetric gauge for SU(2) [22,31-34] and $ SU(3) $ [35-40] gauge theories is a suitable choice as will be discussed below.
2

-->

A.1.Gauge invariant description of a coherent SU(2) condensate

In the SU(2) gauge theory, due to the local isomorphism of the isotopic SU(2) gauge group and the SO(3) group of spatial 3-rotations, the unique (up to a rescaling) SU(2) YM configuration can be parameterized in terms of a scalar time-dependent spatially-homogeneous field (see e.g. Refs. [33,34,41-43]). Indeed, in the case of the SU(2) Yang-Mills fields

$\tag{A2}{A_k} = {A_{ak}}\frac{{{\sigma _a}}}{2}{\mkern 1mu} ,\qquad \left[\frac{{{\sigma _a}}}{2},{\mkern 1mu} \frac{{{\sigma _b}}}{2}\right] = i{\epsilon _{abc}}\frac{{{\sigma _c}}}{2}{\mkern 1mu} ,$

one employs the polar decomposition (symmetric gauge) [33,34] of the gauge field

$\tag{A3}{A_{ak}}(q,S) = {O_{ai}}(q){\mkern 1mu} {S_{ik}} - \frac{1}{{2\bar g}}{\epsilon _{abc}}{\mkern 1mu} {(O(q){\partial _k}{O^T}(q))_{bc}}{\mkern 1mu} ,$

into an orthogonal matrix $ O $, depending on three gauge angles, and a physical symmetric (positively definite) tensor field $ S_{ik} $ ($ \epsilon_{ijk}\,S_{jk} = 0 $) with two spatial indices consisting of one spin-0 component $ S_{ik}^{(0)} = (1/3)\delta_{ik}\,\mathrm{Tr}\,S $ and five spin-2 components. Using the symmetric gauge, one therefore obtains a unique and gauge-invariant decomposition of the gauge field into a spatially homogeneous isotropic part (the YM condensate) and non-isotropic/non-homogeneous parts (the YM waves), namely

$\tag{A4}\begin{split}{A_{ak}}(t,\vec x) =& {\delta _{ak}}U(t) + {{\tilde A}_{ak}}(t,\vec x){\mkern 1mu} ,\\\langle {{\tilde A}_{ik}}(t,\vec x)\rangle =& \int {{{\rm d}^3}} x{{\tilde A}_{ik}}(t,\vec x) = 0{\mkern 1mu} ,\end{split}$

where the YM condensate is positively definite $ U(t)>0 $. In the QFT formulation, the inhomogeneous YM wave modes $ \widetilde{A}_{ik} $ are interpreted as YM quanta (e.g. gluons), while $ U(t) $ contributes to the ground state of the theory. This is consistent with taking the expectation value of the fields on a coherent state $ |\alpha\rangle $ that encodes the symmetry of the background and represents a long-wavelength cosmological condensate state $ \langle \alpha| A_{ak} | \alpha \rangle = \delta_{ak} U(t) $, having picked at $ \vec k = 0 $ the macroscopic semiclassical-state.
In the absence of gravity, the spatially homogeneous isotropic part $ \delta_{ik}U(t) $ of Eq. (A4) satisfies the classical YM equations

$\tag{A5}{(\dot U)^2} + {\bar g^2}{\mkern 1mu} {U^4} = {\rm{const}}{\mkern 1mu} ,$

which can be integrated analytically [22]. In the FLRW Universe, the classical YM condensate behaves as radiation medium with $ a(t)\propto t^{1/2} $ and $ p_{\rm YM} = \epsilon_{\rm vac}/3 $, characteristic features of the classical YM field behaviour. The semi-classical dynamics of the homogeneous SU(2) condensate and the small YM waves has been thoroughly studied in the Minkowski spacetime in Ref. [31].
From the group-theoretical viewpoint, the separation into spatially homogeneous and inhomogeneous components (A4) in the Minkowski spacetime has certain similarities with an analogical procedure in the conventional QCD instanton theory in Euclidean spacetime [44-46]. Namely, in both cases one performs a mapping of 3-space onto SU(2) subgroup elements of the color SU(3)$ _c $. Moreover, as will be discussed below, both cases exhibit remarkable dynamical similarities (or duality).
Notice furthermore that the homogeneous YM condensate can be introduced for every gauge group, which contains at least one SU(2) subgroup (e.g. SU(N)), whereas the condensate embedding or the extraction procedure can be different in every case. These condensates obey similar equations of motion which may differ by a rescaling of the coupling constant. Hence, the frequency of the classical YM condensate oscillations for a SU(N) gauge theory depends on integer N, while the overall dynamical properties are the same.
2 -->

A.2.Inclusion of matter fields: effective QCD energy-momentum tensor

The effective QCD energy-momentum tensor that includes also quarks is expressed by relations that are similar to ones derived for the case of pure gluodynamics. The only difference between the two expressions arises in the effective beta-function coefficient b. When only gluons are taken into account, the one-loop $ \beta $-function coefficient of the pure gluodynamics provides the value $ b = b(0) = 11 $. When quarks are also included, the effective value is found to be [47]

$\tag{A6}{b_{{\rm{eff}}}} = b + 8{L_g}({m_u} + {m_d} + {m_s}) \simeq 9.6{\mkern 1mu} ,$

in which $ L_{g} = (1500\pm 300\, {\rm MeV})^{-1} $ is the correlation length of the fluctuations, recovered by the experimental data on quark and gluon condensates (and supported by the lattice QCD calculations), and with a value close to the minimal scale of quantum-topological fluctuations that contribute to the QCD vacuum. It is indeed renown that for quantum-topological quark-gluon fluctuations [46] the equality holds

$\tag{A7}\begin{split}&\langle 0|:\bar ss:|0\rangle \simeq \langle 0|:\bar uu:|0\rangle \langle 0|:\bar uu:|0\rangle \\ &\quad= - \langle 0|:\frac{{{\alpha _s}}}{\pi }F_{\mu \nu }^aF_a^{\mu \nu }:|0\rangle {L_g} = - {(225 \pm 25{\mkern 1mu} {\rm{MeV}})^3}{\mkern 1mu} .\end{split}$

From a phenomenological point of view, we may deploy reasonable assumptions, and imagine non-perturbative quantum-wave (hadron) fluctuations to occur at the same space-time scales as quantum topological fluctuations. This implies that they should satisfy a functional relation in analogy to (A7). The operator relation between quark and gluon fluctuations can be then established in terms of the trace of the quark energy-momentum tensor, once the vacuum average $ \frac{1}{4} \langle 0 | T^\mu_{\ \mu \,{\rm (QCD)}} | 0\rangle $ is performed, and thus the topological contribution to the energy density recovered. The trace is hence derived from the trace anomaly [48-50], namely

$\tag{A8}T_{\;\mu {\kern 1pt} ({\rm{QCD}})}^\mu = \frac{{\beta (g_s^2)}}{2}{\mkern 1mu} F_{\mu \nu }^aF_a^{\mu \nu } + \sum\limits_{q = u,d,s} {{m_q}} \bar qq{\mkern 1mu} .$

The characteristic topological instanton-type contribution to the energy density of the QCD vacuum, here denoted with $ \epsilon^{\rm QCD}_{\rm top.} $, is then recovered to be

$\tag{A9}\begin{split}\varepsilon _{{\rm{top}}.}^{{\rm{QCD}}} =& - \frac{9}{{32}}\langle 0|:\frac{{{\alpha _s}}}{\pi }{\mkern 1mu} F_{\mu \nu }^aF_a^{\mu \nu }:|0\rangle + \frac{1}{4}[\langle 0|:{m_u}\bar uu:|0\rangle \\ &+ \langle 0|:{m_d}\bar dd:|0\rangle + \langle 0|:{m_s}\bar ss:|0\rangle ]\\& \simeq - (5 \pm 1){\mkern 1mu} \times {10^9}{\mkern 1mu} {\rm{MeV}}{\mkern 1mu} .\end{split}$

This expression is due to the gluons and the light sea u, d, s quark contributions, providing the maximal value of the topological contribution to the QCD vacuum energy density. Nonetheless, thanks to our assumptions both the relations (A7) and (A9) are also valid for the quantum-wave contributions, hence providing the effective quark contribution to the QCD energy momentum tensor operators, i.e.

$\tag{A10}T_{\;\nu {\kern 1pt} ({\rm{q}}),{\kern 1pt} {\rm{eff}}}^\mu = \frac{{8{L_g}}}{{b(3)}}({m_u} + {m_d} + {m_s}){\mkern 1mu} T_{\;\nu {\kern 1pt} ({\rm{g}}),{\kern 1pt} {\rm{eff}}}^\mu {\mkern 1mu} ,$

with obvious meaning of the labels. From this latter, and considering the expression of the operator energy-momentum tensor of the gluon field

$\tag{A11}T_{\;\nu {\kern 1pt} ({\rm{q}}),{\kern 1pt} {\rm{eff}}}^\mu = \frac{{8{L_g}}}{{b(3)}}({m_u} + {m_d} + {m_s}){\mkern 1mu} T_{\;\nu {\kern 1pt} ({\rm{g}}),{\kern 1pt} {\rm{eff}}}^\mu {\mkern 1mu} ,$

one can recover in the one-loop approximation the phenomenologically motivated complete QCD energy-momentum tensor

$\tag{A12}\begin{split}T_{\;\nu {\kern 1pt} ({\rm{QCD}})}^\mu \simeq & \frac{{{b_{{\rm{eff}}}}}}{{32{\pi ^2}}}\left( { - F_{\nu \rho }^a{\mkern 1mu} F_a^{\mu \rho } + \frac{1}{4}\delta _\nu ^\mu {\mkern 1mu} F_{\mu \nu }^a{\mkern 1mu} F_a^{\mu \nu }} \right)\ln \frac{{e{\cal J}}}{{{\lambda ^4}}}\\ &- \delta _\nu ^\mu {\mkern 1mu} \frac{{{b_{{\rm{eff}}}}}}{{128{\pi ^2}}}{\mkern 1mu} F_{\mu \nu }^a{\mkern 1mu} F_a^{\mu \nu }{\mkern 1mu} ,\end{split}$

with $ \mathcal{J} = -F^a_{\mu\nu}\, F_a^{\mu\nu}/\sqrt{-g} $, the constant $ \lambda $ related to the QCD scale parameter by a parameter $ \xi $, namely $ \lambda\equiv \xi \Lambda_{\rm QCD} $, e denoting the base of the natural logarithm and $ b_{\rm eff} $ specified in (A6).

Now we come to an analysis of the equations of motion for physical time evolution of the homogeneous YM condensate in the cosmological environment. For this purpose, we first consider the effective QCD theory in the one-loop approximation, as in Appendix A.2. As was shown in Ref. [51], an extrapolation of the one-loop approximated Lagrangian of the SU(2) gauge theory into deeply infrared (strongly-coupled) regime for the YM vacuum is numerically justified by a comparison to the all-loop result. This argument enables us to pursue the same path for QCD revealing the most important features of the generic non-perturbative QCD. Besides, the results of such a model can be useful to describe the properties of the homogeneous YM condensates at large $ \mathcal{J} $, away from the ground state relevant for cosmological QCD phase transition.
2

-->

B.1.Einstein-YM equations

We take a simplistic approach assuming that before the QCD transition epoch the gluon-field energy density is dominated by positively-valued chromoelectric components, while negatively-valued chromomagnetic components are negligibly small. For convenience, in what follows we re-label the corresponding energy-momentum tensor of the homogeneous YM condensate as $ T^{\nu}_{\mu} \to T^{\nu,{\rm U}}_{\mu} $. As discussed in Ref. [51], both components approach asymptotic ground-state attractors, with exactly opposite densities, thus yielding the exact cancellation in the IR limit of the theory. In this section, we look at the dominant chromoelectric solution in the quark-gluon plasma phase, and study its real time evolution relevant for the QCD transition epoch, away from the asymptotic ground-state attractor.
By the variational principle, one obtains the EYM system of operator equations of motion in a non-trivial spacetime

$\tag{B1}\begin{split}\frac{1}{\varkappa}\left( {R_\mu ^\nu - \frac{1}{2}\delta _\mu ^\nu R} \right) =& \frac{b}{{32{\pi ^2}}}\frac{1}{{\sqrt { - g} }}\Biggl[\left( - {\cal F}_{\mu \lambda }^a{\cal F}_a^{\nu \lambda }+ {\mkern 1mu} \frac{1}{4}\delta _\mu ^\nu {\cal F}_{\sigma \lambda }^a{\cal F}_a^{\sigma \lambda }\right)\\& \ln \frac{{e|{\cal F}_{\alpha \beta }^a{\cal F}_a^{\alpha \beta }|}}{{\sqrt { - g} {\mkern 1mu} {\lambda ^4}}} - \frac{1}{4}\delta _\mu ^\nu {\mkern 1mu} {\cal F}_{\sigma \lambda }^a{\cal F}_a^{\sigma \lambda }\Biggl] ,\end{split}\!\!\!\!\!\!\!\!\!\!\!\!$

$\begin{split}\left( {\frac{{{\delta ^{ab}}}}{{\sqrt { - g} }}{\partial _\nu }\sqrt { - g} - {f^{abc}}{\cal A}_\nu ^c} \right)\left( {\frac{{{\cal F}_b^{\mu \nu }}}{{\sqrt { - g} }}{\mkern 1mu} \ln \frac{{e|{\cal F}_{\alpha \beta }^a{\cal F}_a^{\alpha \beta }|}}{{\sqrt { - g} {\mkern 1mu} {\lambda ^4}}}} \right) = 0{\mkern 1mu} .\end{split}$

Notice that from now on, in all the derivations we will perform a rescaling of the gluon condensate, namely $ \bar{g}\, U(t) \to U(t) $. The energy density and the pressure of other forms of matter are irrelevant for the discussion of the dynamical properties of the gluon condensate in the early Universe (in particular, for cosmological inflation), and will thus be omitted in practical calculations.
In what follows, we work in the flat FLRW conformal metric, characterized by the relations

$\sqrt { - g} = {a^4}(\eta ){\mkern 1mu} ,\qquad t = \int a (\eta ){\rm d}\eta {\mkern 1mu} ,$

the comoving time being defined in $ {\rm d}s^2 = {\rm d}t^2-a^2(t)\vec{{\rm d}x}^2 $. Besides this, we neglect quantum-wave fluctuations $ \widetilde{S}_{ak} $, assuming that the homogeneous gluon condensate $ U(t) $ dominates at considered spacetime scales. Under these conditions, the system of equations of motion describing conformal time evolution of the gluon condensate $ U = U(\eta) $ and the scale factor $ a = a(\eta) $ read

$\tag{B2}\begin{split}&\frac{6}{\varkappa}\frac{{a''}}{{{a^3}}} = T_\mu ^{\mu ,{\rm{U}}},\quadT_\mu ^{\mu ,{\rm{U}}} = \frac{{3b}}{{16{\pi ^2}{a^4}}}\left[{(U')^2} - \frac{1}{4}{U^4}\right] ,\\&\frac{\partial }{{\partial \eta }}\left(U'{\mkern 1mu} \ln \frac{{6e|{{(U')}^2} - \frac{1}{4}{U^4}|}}{{{a^4}{\lambda ^4}}}\right) + {\mkern 1mu} \frac{1}{2}{U^3}{\mkern 1mu} \ln \frac{{6e|{{(U')}^2} - \frac{1}{4}{U^4}|}}{{{a^4}{\lambda ^4}}} = 0 .\end{split}$

An additional coefficient 1/2 appears in front of the QCD coupling constant—which has been absorbed into the definition of the gluon field—as compared to the SU(2) condensate case considered in Ref. [22]. The first integral of Eq. (B2) is the Einstein (0,0)-equation, and reads

$\tag{B3}\begin{split}\frac{3}{\varkappa}\frac{{{{(a')}^2}}}{{{a^4}}} =& T_0^{0,{\rm{U}}}{\mkern 1mu} ,\\T_0^{0,{\rm{U}}} =& \frac{{3b}}{{64{\pi ^2}{a^4}}}{\mkern 1mu} \left(\left[{(U')^2} + \frac{1}{4}{U^4}\right]{\mkern 1mu} \ln \frac{{6e|{{(U')}^2} - \frac{1}{4}{U^4}|}}{{{a^4}{\lambda ^4}}}\right.\\ &\left.+ {\mkern 1mu} {(U')^2} - \frac{1}{4}{U^4}\right) .\end{split}$

2 -->

B.2.Asymptotic attractor solution and Z2 symmetry restoration

In full analogy to the SU(2) condensate case considered in Ref. [22], in the QCD case the system of equations (B2) and (B3) has the exact solution corresponding to the vanishing logarithm or, equivalently, satisfies the transcendent equation

$\tag{B4}\begin{split}|Q| = 1{\mkern 1mu} ,\qquad Q \equiv & \frac{{32{\pi ^2}e}}{{11{{(\xi {\Lambda _{{\rm{QCD}}}})}^4}}}T_\mu ^{\mu ,{\rm{U}}*}\\ = &\frac{{6e[{{(U')}^2} - \frac{1}{4}{U^4}]}}{{{a^4}{{(\xi {\Lambda _{{\rm{QCD}}}})}^4}}}{\mkern 1mu} ,\end{split}$

which yields two distinct cases $ Q = \pm 1 $.
One of these cases corresponding to $ Q = +1 $ has been discussed in the broad literature before (see e.g. Refs. [22,32,52-55]). In particular, it provides a positive constant energy density of the gluon condensate

$\tag{B5}{(U')^2} - \frac{1}{4}{U^4} > 0{\mkern 1mu} ,\;\;T_0^{0,{\rm{U}}*} \equiv \frac{{3b}}{{64{\pi ^2}}}\frac{{{{(\xi {\Lambda _{{\rm{QCD}}}})}^4}}}{{6e}} > 0,$

and thus has been considered as a potential driver of early-time acceleration epochs in cosmology.
The second solution, characterized by $ Q = -1 $, has not been sufficiently discussed in the literature. It is worth noticing that it appears due to the symmetry of the corresponding RG equation with respect to the $ \mathcal{J} \leftrightarrow -\mathcal{J} $ interchange—see e.g. Ref. [51]. Thus, this must be considered on the same footing as the solution with $ Q = +1 $. Such a contribution to the energy density of the universe is negative and has the same absolute value as in Eq. (B5), i.e.

$\tag{B6}{(U')^2} - \frac{1}{4}{U^4} < 0,\;\;T_0^{0,{\rm{U}}*} \equiv - \frac{{3b}}{{64{\pi ^2}}}{\mkern 1mu} \frac{{{{(\xi {\Lambda _{{\rm{QCD}}}})}^4}}}{{6e}} < 0.$

Provided such a negative-energy vacuum solution is stable, a positive CC-term should be present to compensate this negative contribution in order to comply with cosmological observations. One may notice, however, that the corresponding negative energy density for $ \xi\simeq 4 $ coincides with the quantum-topological term $ \epsilon^{\rm QCD}_{\rm top.} $ provided in Eq. (A9). As was advocated in Ref. [51], such a cancellation is due to the emergence of a $ \mathbb{Z}_2 $ symmetry in the asymptotic ground state. The metastable positive-energy vacuum configurations are trapped in pockets of 3-space, which are separated from the stable chromomagnetic pockets by means of the infinite barriers, the so called domain walls.
Such a compensation mechanism grossly reduces or eliminates the QCD vacuum effect on the macroscopic late-time Universe expansion. Indeed, under the conditions specified in Ref. [51], the macroscopic evolution of the Universe reduces to the standard Friedmann equation, driven only by matter fields and a small uncompensated observable term $ \epsilon_{\rm CC}\lll T^{0,{\rm U}*}_{0} $, while the evolution of the gluon condensate happens at characteristic microscopic scales corresponding to the QCD confinement scale $ \Lambda_{\rm QCD} $, i.e.

$\tag{B7}\begin{split}\frac{3}{\varkappa}\frac{{{{(a')}^2}}}{{{a^4}}} =& \epsilon + {\epsilon _{{\rm{CC}}}}{\mkern 1mu} ,\\{(U')^2} - \frac{1}{4}{U^4} =& {a^4}\frac{{{{(\xi {\Lambda _{{\rm{QCD}}}})}^4}}}{{6e}}{\mkern 1mu} ,\qquad \xi \simeq 4{\mkern 1mu} .\end{split}$

Consequently, such a relatively slow macroscopic evolution of the Universe for $ a = a(\eta) $, and the rapid fluctuations of the gluon condensate $ U = U(\eta) $ at the characteristic QCD timescale, get practically separated and are independent from each other.
In the present Universe, for which we pick $ a\equiv a_0 = 1 $, the exact (implicit) partial solutions for the homogeneous gluon condensate read

$\tag{B8}\begin{split}Q =& \pm 1{\mkern 1mu} ,\qquad \int_{{{\tilde U}_0}}^{\tilde U} {\frac{{{\rm d}u}}{{\sqrt {\frac{1}{4}{u^4} \pm 1} }}} = \tilde \eta {\mkern 1mu} ,\\\tilde U =& U\frac{{{{(6e)}^{1/4}}}}{{4{\Lambda _{{\rm{QCD}}}}}}{\mkern 1mu} ,\tilde \eta {\mkern 1mu} = \eta \frac{{4{\Lambda _{{\rm{QCD}}}}}}{{{{(6e)}^{1/4}}}}{\mkern 1mu} ,\end{split}$

corresponding to (B5) and (B6) solutions. Thus, the cosmological evolution of the gluon field in its ground state can be interpreted as a regular sequence of quantum tunneling transitions through the “time-barriers” represented by the regular singularities in the quantum vacuum solution of the effective YM theory. In this sense, the homogeneous gluon condensate in Minkowski spacetime is analogous to the topological condensate in the instanton theory of the QCD vacuum in Euclidean spacetime. This latter can be interpreted in terms of spatially-inhomogeneous gluon field fluctuations, which are induced by the quantum tunneling of the field through topological (spatial) barriers between different classical vacua.
It is worth noticing that the well-known't Hooft-Polyakov monopole [56,57] is analogous to the classical YM condensate (A4). While the monopole is introduced by means of an antisymmetric matrix with mixed Lorentz-isotopic indices in the Weyl gauge, the considered YM condensate provides a symmetric analogue of this solution. The singularities emerge at the level of the 4-potential of the YM condensate $ U(t) $, while all the observable quantities (such as pressure and energy density), are finite. This situation reminds the Dirac monopole solution, where a singularity in the 4-potential emerges along the Dirac string, while the magnetic flux is finite.
As was noticed above, the compensation may not be exact at the early stages of the Universe evolution, and may be fulfilled only at asymptotically large times and on the average. Notice furthermore that the macroscopic evolution of the present Universe, in this case, is only affected by a small remnant of the above cancellation, the observable cosmological constant term. The latter can be generated e.g. by an uncompensated quantum-gravity correction to the ground state energy in QCD, which may explain both its absolute value and sign, as was argued in Ref. [47]. Besides that, a positive uncompensated component of the YM vacuum, corresponding to the positive-valued homogeneous YM (non-QCD) condensate contribution in the early Universe, can be a natural cause for the cosmic inflation. Below, we will briefly discuss a possible scenario of the early-time accelerated expansion naturally provided by the homogeneous YM condensate $ U(t) $.
2 -->

B.3.Numerical analysis and tracker solutions of the EYM system

We now consider a deviation from the exact partial solution given by Eq. (B8), and study the general solution of the equations of motion (B2) and (B3), numerically. Let us first choose the subset of initial conditions satisfying $ Q_0\equiv Q(t = t_0)>1 $, and discuss the results of the numerical analysis qualitatively.
For this choice of initial conditions, Fig. B1 (upper panel) illustrates the physical time evolution of the total energy density (in dimensionless units) of the homogeneous gluon condensate $ U = U(t) $
FigureB1. (color online) We illustrate, as functions of the physical time $t = \int a{\rm d}\eta$, and using units of the characteristic time scale $\Lambda_{\rm QCD}^{-1}$, respectively the total energy density $T^0_0(t)$ of the homogeneous gluon condensate (upper panel), the trace of the total QCD energy-momentum tensor $T^\mu_\mu(t)$ (middle panel) and the logarithm of the scale factor $a(t)$ (lower panel). We indicate the total energy density and the trace values for $Q_0\equiv Q(t_0) = 1$, respectively, with horizontal lines in the up and middle panels, where the initial conditions are chosen as $U_0 = 0$, $\dot{U}_0 = (\xi \Lambda_{\rm QCD})^2/\sqrt{3e}$, $Q_0>1$, $\xi\simeq 4$, and the gravitational constant is set to be $\varkappa = 10^{-7} {\rm MeV}^{-2}$, for simplicity of the numerical analysis. We plot both $T^0_0(t)$ and $T^\mu_\mu(t)$ in dimensionless units, and rescale them by $\Lambda_{\rm QCD}^4$. It is evident that the amplitude of the quasi-periodic oscillations of $Q = Q(t)$ happens to decrease at large $t\gg \Lambda_{\rm QCD}^{-1}$, and to approach asymptotically unity.

$\tag{B9}T_0^0(t) \equiv \bar \epsilon + T_0^{0,{\rm{U}}}(t){\mkern 1mu} ,$

where $ T^{0,{\rm U}}_0 $ and $ \bar{\epsilon} $ are given respectively by Eq. (B3) and by $ \bar{\epsilon} = \epsilon_{\rm top.}^{\rm (QCD)}+\epsilon_{\rm CC} $, the $ \epsilon_{\rm CC} $ term denoting the cosmological constant contribution. In Fig. B1 (middle panel), we show the corresponding result for the trace of the total gluon energy-momentum tensor

$\tag{B10}T_\mu ^\mu (t) \equiv 4\bar \epsilon + T_\mu ^{\mu ,{\rm{U}}}(t){\mkern 1mu} ,$

in dimensionless units, and the corresponding solution for the logarithm of the scale factor is given in Fig. B1 (lower panel). The period of the $ T^\mu_\mu(t) $ oscillations is practically time independent, which can also be proven analytically (see below), while a small residual time-dependence appears due to a possibly large deviation from $ Q = 1 $. In this analysis we used $ \xi\simeq 4 $, providing the exact compensation of the QCD vacuum energy density at $ t\gg t_0 $, such that a change of $ \xi $ will only affect the asymptotic values of $ T^0_0(t) $ and $ T^\mu_\mu(t) $ at large t.
Although the amplitude of the condensate $ U(t) $ possesses quasi-periodic singularities, as is seen in the upper panel of Fig. B2, the evolution of its energy density $ T^{0,{\rm U}}_0(t) $ as well as pressure (or $ T^{\mu,{\rm U}}_\mu(t) $) remains continuous in time. One immediately notices that the general solution asymptotically reaches the analytic solution corresponding to Eq. (B5). This happens after a number of oscillations of function $ Q(t) $ whose amplitude approaches unity at large physical times t, $ Q(t\to \infty)\Rightarrow 1 $, for any initial conditions satisfying $ Q_0>0 $. During such a relaxation regime, the total energy density of the QCD vacuum continuously decreases and eventually vanishes in the asymptotic limit $ t\gg t_0 $. Note, this regime is accompanied by a decelerating expansion of the Universe.
FigureB2. (color online) We show the homogeneous QCD condensate amplitude oscillations $U = U(t)$, manifesting quasi-periodic singularities in the physical time $t = \int a{\rm d}\eta$, for the solutions labelled by $Q(U) = 1$ and the $Q(U) = -1$. These are shown respectively in the upper and lower panels, using units of the characteristic time scale $\Lambda_{\rm QCD}^{-1}$. The spikes that are displayed are localized in time-lapse, along the space-like directions, and must be interpreted as new solitonic solutions. In a previous work, we dubbed these solutions chronons, or $\chi$-solutions.

The same quantities have also been studied in the opposite case of initial conditions, i.e. for $ 0<Q_0<1 $, and the results (without accounting for a $ \bar{\epsilon} $ term) are illustrated in the lower panel of Fig. B2. In this case the general solution asymptotically approaches the de-Sitter regime as well, in full analogy with the $ Q_0>1 $ case. A qualitatively similar situation is realized for $ Q_0<0 $ as well. The de-Sitter solution, therefore, appears as an attractor (or tracker) solution of the EYM system. This provides a dynamical mechanism for the elimination of the gluon vacuum component of the ground state energy of the Universe, asymptotically at macroscopic space-time scales and for arbitrary initial conditions and parameters of the model. This fact provides a generalization of the basic result of Ref. [22] to the case of arbitrary initial conditions and to a gauge group possessing a $ SU(2) $ subgroup.

In Section B.3, the system of equations (B2) and (B3) was investigated numerically in the general case. We assumed arbitrary initial conditions and found the universal asymptotics corresponding to the (partial) exact analytic solutions with minimal energy $ |Q| = 1 $. In this Section, we will construct the general analytic approximate solution of this system, in the case of $ Q(t) $ being not too far from unity which reproduces the results of numerical simulations and provides us with an extra important information about the cosmological evolution of the QCD (or YM, in general) vacuum. In particular, such a solution will provide an additional proof for the observed asymptotic behaviour of $ a(t) $, $ T_0^0(t) $ and $ T_{\mu}^{\mu}(t) $. By comparison with the numerical analysis performed above, we will notice that all the characteristic features of the approximate analytic solution, corresponding to the one-loop effective action, turn out to be valid also in the case of a large deviation of $ Q(t) $ from unity, as well as beyond the perturbation theory. In this sense, the aforementioned numerical analysis is of important guidance, and supplements the approximated analytic results obtained below.
2

-->

C.1.Period of oscillations and relaxation time

One could expect that the oscillatory behavior of the $ U(\eta) $ function found analytically for $ |Q| = 1 $, will also be present for $ Q_0\neq 1 $, as long as the deviation of $ Q_0 $ from unity is not very large. Indeed, if the initial energy density is not too large than the period of oscillations of the positive-energy solution for the condensate $ U = U(\eta) $ in conformal time, $ T_{\eta}^{>} $ can be roughly estimated using Eq. (B8) as follows

$\tag{C1}\begin{split}T_\eta ^ > \simeq \frac{{4k{\mkern 1mu} {{(6e)}^{1/4}}}}{{a\xi {\Lambda _{{\rm{QCD}}}}}}{\mkern 1mu} ,\qquad k \equiv & \int_0^\infty {\frac{{{\rm d}u}}{{\sqrt {1 + \frac{1}{4}{u^4}} }}} \\ =& \frac{{\Gamma {{(1/4)}^2}}}{{2\sqrt {2\pi } }} \approx 2.622{\mkern 1mu} .\end{split}$

The corresponding period in physical time reads approximately

$\tag{C2}{T_t} \equiv T_t^ > \simeq \frac{{4k{\mkern 1mu} {{(6e)}^{1/4}}}}{{\xi {\Lambda _{{\rm{QCD}}}}}} \simeq \frac{{5.3}}{{{\Lambda _{{\rm{QCD}}}}}}{\mkern 1mu} ,$

which is close to the numerical result discussed above. Notice that the period of YM condensate oscillations corresponding to the negative-energy solution (B6) $ T_t^{<} $ differs from Eq. (C2), being

$\tag{C3}T_t^ < \simeq \frac{{\sqrt 2 k{\mkern 1mu} {{(6e)}^{1/4}}}}{{\xi {\Lambda _{{\rm{QCD}}}}}} \simeq \frac{{1.86}}{{{\Lambda _{{\rm{QCD}}}}}}{\mkern 1mu} .$

The conformal time derivative of the gluon condensate energy density for $ Q_0 > 1 $ in QCD, given by

$\tag{C4}\frac{{\partial T_0^{0,{\rm{U}}}}}{{\partial \eta }} = - \frac{{33\varkappa }}{{16{\pi ^2}}}\frac{{a'}}{{{a^5}}}[{U'^2} + \frac{1}{4}{U^4}]\ln \frac{{6e|{{(U')}^2} - \frac{1}{4}{U^4}|}}{{{a^4}{{(\xi {\Lambda _{{\rm{QCD}}}})}^4}}},$

is negative in the initial moment of time $ t = t_0 $ and its absolute value decreases such that $ T_0^{0,{\rm U}}(t) $ indeed approaches the value $ T_0^{0,{\rm U}^*} $, corresponding to the exact analytical (de-Sitter) solution given by Eq. (B5). On dimensional grounds, a rough estimate for the relaxation time of the YM condensate energy density is approximately given by

$\tag{C5}{t_r} \simeq \frac{1}{{\sqrt {\varkappa {\epsilon _0}} }}{\mkern 1mu} ,\qquad {\epsilon _0} \equiv T_0^0(t = {t_0}){\mkern 1mu} .$

As will be shown analytically below, the scale factor approaches the de-Sitter solution at late times

$\tag{C6}t > \frac{1}{{\sqrt {\varkappa {\epsilon _{{\rm{CC}}}}} }} \gg {t_r}{\mkern 1mu} .$

Under the viable hypothesis of co-existence of chromoelectric and chromomagnetic condensates around the QCD phase transition epoch, their initial values at $ t = t_0 $ before the hadronization of the quark-gluon plasma may be different (due to different formation energy scales, for example), and thus they do not identically compensate each other. Due to the attractor nature of the corresponding solutions the compensation may happen only at late times (in particular, after the Universe gets hadronized), at typical time scales $ t\gg t_r $. While the compensation can be (nearly) exact when the ground-state energy density is averaged over macroscopically large spacetime separations, such a compensation may not be exact locally, yielding a possibly non-vanishing effect (within a typical hadron scale one is sensitive to) in e.g. hadronic reactions.
Now let us consider the oscillations of the trace $ T_{\mu}^{\mu,{\rm U}}(t) $, providing an information about the pressure of the gluon condensate. The unobservable function $ U(t) $, and hence the combination $ (U')^2+\displaystyle\frac{1}{4}U^4 $, oscillate with non-physical quasi-periodic singularities. Nonetheless, the energy density $ T_0^{0,{\rm U}}(t) $, which is a physical observable given by Eq. (B3), has to be continuous, as confirmed by the numerical analysis. The latter condition can only be realized if $ Q = Q(t) $ oscillates as well with the same period as the condensate, and periodically reaches unity (such that its logarithm vanishes). In such a way $ Q = Q(t) $ compensates the corresponding singularities in the coefficient $ (U')^2+\displaystyle\frac{1}{4}U^4 $ term. Besides that, the $ Q = Q(t) $ function can not cross unity, since otherwise it would lead to a non-continuity in the behavior of $ T_0^{0,{\rm U}}(t) $, due to a sign change in the first term in Ref. [32] when $ Q = 1 $. From these qualitative arguments it is clear that the time evolution of $ Q = Q(t) $ function satisfies the constraints $ 0<Q(t)\leqslant1 $ or $ Q(t)\geqslant1 $ for the initial conditions $ 0<Q_0<1 $ or $ Q_0>1 $, respectively. The period of its oscillations can then be estimated as

$\tag{C7}T \equiv \frac{{{T_t}}}{2} \simeq \frac{{2k{\mkern 1mu} {{(6e)}^{1/4}}}}{{\xi {\Lambda _{{\rm{QCD}}}}}}{\mkern 1mu} ,$

and its relaxation time is the same as for $ T_0^{0,{\rm U}}(t) $ given in Eq. (C5). Consequently, these basic features of the general solution of the EYM equations for the YM condensate (as the periods of oscillations and the relaxation times of the condensate), its energy density and the pressure can be described qualitatively, without a reference either to the numerical calculations or to a particular choice of model parameters.