*Article* **Invisible QCD as Dark Energy**

#### **Andrea Addazi 1,2,\*, Stephon Alexander <sup>3</sup> and Antonino Marcianò 4,5,\***


Received: 23 May 2020; Accepted: 28 May 2020; Published: 31 May 2020

**Abstract:** We account for the late time acceleration of the Universe by extending the Quantum Chromodynamics (QCD) color to a *SU*(3) invisible sector (IQCD). If the Invisible Chiral symmetry is broken in the early universe, a condensate of dark pions (dpions) and dark gluons (dgluons) forms. The condensate naturally forms due to strong dynamics similar to the Nambu–Jona-Lasinio mechanism. As the Universe evolves from early times to present times the interaction energy between the dgluon and dpion condensate dominates with a negative pressure equation of state and causes late time acceleration. We conclude with a stability analysis of the coupled perturbations of the dark pions and dark gluons.

**Keywords:** dark energy; non-Abelian gauge theory; condensate

#### **1. Introduction**

A confluence of cosmological data tell us that the Universe is currently accelerating (see e.g., Refs. [1–3] and references therein). While the data can be explained with a cosmological constant, it is also possible that the Universe is dominated by a form of Dark Energy (DE) with a negative pressure barotropic index *w*∼ −1. Another possibility is that general relativity is modified in the infrared (IR) admitting self accelerating solutions, and such models are still under active research Ref. [4]. In this work we take the position that the Einstein-Hilbert action is the correct IR description of gravity and the cosmological constant is zero, by some yet unknown mechanism. We present a model where DE emerges from an Invisible QCD (IQCD) sector due to the interaction of invisible pions and gluons which were present in the early Universe.

There are two ways to theoretically motivate an Invisible QCD sector: it has been known for a long time that the *E*<sup>8</sup> × *E*<sup>8</sup> Heterotic Superstring Theory, both in the weak and strong coupling limit necessarily gives rise to non-abelian gauge theories that do not directly couple to the standard model *SU*(3)*<sup>c</sup>* × *SU*(2)*<sup>W</sup>* × *U*(1)*<sup>Y</sup>* gauge group. While String Theory has this property, it is also possible to obtain a dark copy in another way, which corresponds to modifying general relativity by enlarging the spin connection *ω*(*e*) *I J <sup>μ</sup>* to a larger group *G* ∼ *SU*(*N*), which breaks to *SU*(2) × *SU*(*N* − 2) where the *SU*(*N* − 2) factor is identified with the dark sector. This philosophy was pursued in Refs. [5–9]. Here we take a phenomenological perspective and simply assume that this dark sector exists and focus on the cosmological consequences.

As a topic contiguous to the one treated in this investigation, we shall provided a short review on the most relevant studies addressing a possible relation of QCD to dark energy. For sure relevant to our discussion is the perspective proposed by E. Witten, who has shown that the small value of the cosmological constant could be explained in terms of a vacuum state with unbroken supersymmetry. Specifically, as summarized in Ref. [10], Witten has argued that, while compactifying Calabi-Yau manifolds, the string tension in higher dimension can be invoked to cancel the vacuum expectation value of QCD and other fields, and hence reproduce the observed value of the cosmological constant. An extensive review on the same argument was offered in Ref. [11] by R. Bousso. While an analysis more focused on the the link between QCD condensate and the observed value of the cosmological constant is the one contained in Ref. [12], where the authors studied quark and gluon condensates in QCD, as associated to the internal dynamics of hadrons, and then succeeded in showing that these condensates actually provide a vanishing contribution to the cosmological constant.

Related to an extended invisible/mirror framework is the investigation carried out in Ref. [13], which appeared after the first submission of this work on the electronic archives [14]. The analysis in Ref. [13] has remarkable consequences at the level of the mirror symmetry, which some of us further investigated in Ref. [15]. Specifically, the authors of Ref. [13] explored the possible compensation of the negative contribution, due to the QCD vacuum, to the ground state energy density of the Universe, thanks to the (opposite in sign) contribution that arises from the chromomagnetic gluon condensate in the invisible/mirror QCD sector. Although not directly related to the fate of the QCD vacuum, it is nonetheless worth mentioning the study achieved by some of us in Ref. [16], on a physical mechanism similar to the one we developed here. This mechanism led to a quasi-de Sitter expansion of the Universe, which inspired the analysis of this paper, by considering a fields backreaction (in the expanding Universe) that sources an effective cosmological constant. This is due to the interaction energy among the gauge hypercharge field and the fermion field.

In this model, late time acceleration arises from extending the color sector of QCD to have an "invisible-copy". IQCD has similar quantum field theoretic properties of QCD, in that it is confining in the IR. It is well known that pions arise as Goldstone modes associated to Chiral Symmetry Breaking (CSB), and in turn the microphysical description of CSB, the Nambu–Jona-Lasinio mechanism Ref. [17,18], is a strongly coupled version of superconductivity induced by hard gluon exchange. A key feature of our DE model uses the same physics of CSB in the IQCD sector. During the matter and radiation era, the dark pions and gluons have negligible effects. However, we will show that at late times the interaction energy between the dark pions and gluons become more significant because they remain nearly constant, mimicking an effective cosmological constant. Through the consistency of the coupled field equations, this interaction energy naturally leads to late time acceleration and we find an interesting connection between the scale of CSB and the scale of DE.

#### **2. The Theory**

We assume that Chiral symmetry in the IQCD sector is broken at some scale *fD* corresponding to the mass of a dark pion<sup>1</sup> (dpion) *π<sup>D</sup>* . At the renormalizable level, such a hidden sector can communicate with the standard model only through gauge interactions: *SU*(3)*c*, *SU*(2)*<sup>W</sup>* or *U*(1)*Y*—see e.g., Refs. [19–21]. In this work, for clarity of presentation, we assume no gauge coupling to the standard model. However, this model can easily be extended such that the gauge-confined quarks are coupled vectorially to *SU*(2)*W*. The authors of Ref. [21] showed that a dark sector with purely vectorial coupling to the standard model has remarkable universal features. A specific parity symmetry (known as G-parity) acting only on the hidden sector fields is left unbroken and stable weakly interacting dark-pions become a compelling candidate for a dark-matter particle. It would be interesting to revisit the constraints on the coupling to the visible sector that is simultaneously consistent with dark-matter and DE. We find it intriguing that our model has the possibility of connecting late time acceleration to Dark Matter and will pursue this possibility in Ref. [22].

<sup>1</sup> From now on, we will remove the subscript *D* when referring to dark pions.

The *SU*(*N*)*<sup>D</sup>* gauge theory we are considering is assumed to have a behavior similar to QCD. The gauge coupling becomes strong in the IR limit, triggering confinement and chiral symmetry breaking at a scale Λ*D*. Below Λ*D*, the effective theory is described by "pions" which are pseudo-Nambu–Goldstone bosons (pNGBs) associated with the spontaneously broken global flavor symmetry of the hidden sector.

For concreteness and without loss of generality, we consider the subgroup, *SU*(2)*<sup>L</sup>* × *SU*(2)*<sup>R</sup>* → *SU*(2)*<sup>V</sup>* with its gauge field *A<sup>a</sup> <sup>μ</sup>*, where *a*, *b* = 1, 2, 3 and *μ*, *ν* = 0, 1, 2, 3 are for the dark color and space-time indices, respectively. The gauge field strength *F* is

$$F^{a}\_{\mu\nu} = \partial\_{\mu}A^{a}\_{\ \nu} - \partial\_{\nu}A^{a}\_{\ \mu} - \g{\epsilon}\epsilon^{a}\_{\ \nu\epsilon}A^{b}\_{\ \mu}A^{c}\_{\ \nu} \tag{1}$$

where *abc* is the totally antisymmetric Levi-Civita symbol, the structure constant of the *SU*(2) algebra. We are led to consider the most general gauge invariant action coupled to "dark quarks", with masses *mi* and Lagrangian

$$\mathcal{L}\_{\text{IQCD}} \equiv \mathcal{L}\_{\text{A}} + \mathcal{L}\_{\text{D}} = -\frac{1}{4g^2} F^{a}\_{\mu\nu} F^{\mu\nu}\_{a} + \overline{\psi}\_i \left( \imath D\_{\mu} \gamma^{\mu} - m\_i \right) \psi\_i \,. \tag{2}$$

Here *D<sup>μ</sup>* stands for the covariant derivative with respect to the dark SU(3) sector and the gravitational connection, *γ<sup>μ</sup>* =*γ<sup>I</sup> e μ <sup>I</sup>* and the metric field is decomposed in tetrad, namely *<sup>g</sup>μν* = *<sup>e</sup><sup>I</sup> <sup>μ</sup> e<sup>I</sup> ν*, the inverse of which is denoted as *e μ <sup>I</sup>* and the internal *SO*(3, 1) indices of which are *I* = 1, 2...4. Given that we are working in a system where the dpion forms as a result of CSB, the decay constant *fD* is defined through the coupling of the axial current to the dpion. In particular, dpions can be created by the axial isospin currents.

Matrix elements of J<sup>5</sup> *<sup>I</sup> <sup>a</sup>* (*x*) between the vacuum and an on-shell dpion state can be parametrized as (see e.g., Ref. [23])

$$<\langle 0 \vert \beta\_a^{5I}(\mathbf{x}) \vert \pi^b(p) \rangle = -\text{tr}\, \delta\_a^b \left. p^I f\_D \, e^{-i p \cdot \mathbf{x}} \right. \tag{3}$$

where the chiral current is J<sup>5</sup> *<sup>I</sup> <sup>a</sup>* (*x*) = *ψ*(*x*) *γ<sup>I</sup> γ*<sup>5</sup> *τ<sup>a</sup> ψ*(*x*) and |*π*(*p*) is the pion ground-state with normalization *π*(*p* )|*π*(*p*) <sup>=</sup> <sup>2</sup>(2*π*)<sup>3</sup> *<sup>p</sup>*0*δ*3(*<sup>p</sup>* − *p*). Relation (3) can be recast in terms of the dpion field as

$$
\langle 0 | \mathcal{J}\_a^{5I}(\mathbf{x}) | \pi^b \rangle = f\_D \, \delta\_a^b \, \partial^I \pi\_b(\mathbf{x}) \,. \tag{4}
$$

We can rotate the expectation value of 0|J<sup>5</sup> *<sup>I</sup> <sup>a</sup>* (*x*)|*πb* within the internal Lorentz indices' space, so to accomplish an explicitly space-like axial vector with vanishing temporal component. The symmetry of the vacuum state on the FLRW background allows to further reduce Equation (4) to a homogenous axial vector:

0|J<sup>5</sup> *<sup>I</sup> <sup>a</sup>* (*x*)|*πb* <sup>=</sup> *<sup>f</sup>* <sup>2</sup> *D δi <sup>a</sup> π*(*t*, 0), (5)

where *<sup>π</sup>*(*t*)≡||*πa*(*t*)||, with respect to the internal indices. The interaction of the axial current with the gauge field L*int*. *<sup>π</sup>* = *g A<sup>a</sup> μ* J 5 *μ <sup>a</sup>* is therefore consistent with homogenous and isotropic metrics.

The low energy dpion effective Lagrangian reads

$$\mathcal{L}^{0}\_{\pi} = -\frac{1}{2} \partial\_{\mu} \pi\_{a} \partial^{\mu} \pi^{a} + \frac{\lambda}{4} \left( \pi^{a} \,\pi\_{a} - f\_{D}^{2} \right)^{2}. \tag{6}$$

Consequently, the total effective Lagrangian reads

$$\begin{split} \mathcal{L}\_{\text{tot.}} &= \, \mathcal{L}\_{\text{GR}} + \mathcal{L}\_{\text{A}} + \mathcal{L}\_{\pi}^{int.} + \mathcal{L}\_{\pi}^{0} \\ &= \, \, \mathcal{M}\_{p}^{2} \, \mathcal{R} - \frac{1}{4} F\_{\mu\nu}^{a} F\_{a}^{\mu\nu} + \mathcal{g} \, A\_{\mu}^{a} \mathcal{J}\_{a}^{5\mu} + \mathcal{L}\_{\text{B}}^{0} \end{split} \tag{7}$$

in which we have introduced the reduced Planck mass as *M*<sup>2</sup> *<sup>p</sup>* = (8*πG*)−1. Quark fields have been integrated out in the path integral in order to get the effective action. The interaction term L*int*. *<sup>π</sup>* = *g A<sup>a</sup> μ* J 5 *μ <sup>a</sup>* , which entails parity violations of the *SU*(2) subgroup of the dark sector, preserves renormalizability. The total action is S*tot*. = *d*4*x e* L*tot*., with *e* volume density denoting the determinant of the tetrad *e<sup>I</sup> μ*.

#### **3. Field Equations**

Solutions to the field equations that are consistent with a FLRW background can be recovered once a rotationally invariant configuration for the gauge field has been implemented:

$$A^a\_{\;\;\mu} = \begin{cases} \begin{array}{ll} a(t)\,\phi(t)\,\delta^a\_{\;i\;\prime} & \mu = i, \\ 0, & \mu = 0. \end{array} \end{cases} \tag{8}$$

Thanks to Equations (4) and (8), the energy-momentum tensor of the theory is isotropic and homogenous. In particular, the energy-momentum tensor associated to the interaction Lagrangian yields the remarkable feature of having a barotropic index *w* = −1, since energy and pressure densities respectively read

$$\begin{array}{rcl} \rho\_{\text{AJ}} & = & \mathfrak{Z} \, \mathfrak{g} \, f\_{\text{D}}^{2} \, \phi(t) \, \pi(t), \\\ -P\_{\text{AJ}} & = & \mathfrak{Z} \, \mathfrak{g} \, f\_{\text{D}}^{2} \, \phi(t) \, \pi(t), \end{array}$$

*g* denoting above the absolute vale of the coupling constant. This naturally leads towards a de Sitter accelerating phase of the Universe, as soon as the interaction term becomes dominant.

In the rotationally symmetric configuration Equation (8) the gauge field-strength's components simplify and read *F<sup>a</sup>* <sup>0</sup>*<sup>i</sup>* = *<sup>∂</sup>t*(*φ*(*t*)*a*(*t*)*δ<sup>a</sup> <sup>i</sup>* ) and *<sup>F</sup><sup>a</sup> ij* <sup>=</sup> <sup>−</sup>*<sup>g</sup> <sup>a</sup> ij*(*φ*(*t*)*a*(*t*))2, having specified our system in co-moving coordinates *ds*<sup>2</sup> <sup>=</sup>*dt*2−*a*2(*t*) *<sup>d</sup>x*2. Using (8), the total gauge Lagrangian becomes

$$\mathcal{L}^{\mathbf{A}} + \mathcal{L}^{\text{int.}} = \frac{1}{2} \frac{1}{g^2} \left( 3 \left( \dot{\phi} + H\phi \right)^2 - 3 \, g^2 \phi^4 \right) + 3g\phi \, \mathbb{J}(a), \tag{9}$$

where ¯*J*(*a*)<sup>≡</sup> *<sup>f</sup>* <sup>2</sup> *<sup>D</sup> <sup>π</sup>*(*t*). From <sup>1</sup> *a*3 *∂ ∂t*(*a*<sup>3</sup> *<sup>∂</sup>*<sup>L</sup> *∂φ*˙ ) = *<sup>∂</sup>*<sup>L</sup> *∂φ* , we obtain the equation of motion for *φ*, which captures the dynamics of *A<sup>a</sup> <sup>μ</sup>* through Equation (8), namely

$$
\ddot{\phi} + 3H\dot{\phi} + \left(2H^2 + \dot{H}\right)\phi + 2\,\text{g}^2\phi^3 - \text{g}\,\text{J}(a) = 0.\tag{10}
$$

The equation of motion for the dpion field is recovered varying Equation (6), within the assumption of spatial homogeneity. This is plausible, since a previous inflationary epoch of the universe can smooth out the dpion field. In the next section, we show that the dpion field remains homogenous against perturbations.

Using the decomposition in a homogeneous absolute value (in the internal space) times a space-dependent unit vector, *i.e. π<sup>a</sup>* = ||*πa*|| *na* = *π*(*t*) *na*(*x*), we recover for the pion field

$$
\ddot{\pi} + 3H\dot{\pi} + \lambda \,\pi(\pi^2 - f\_D^2) - 3\,\text{g}f\_D^2 \phi = 0 \,\text{.}\tag{11}
$$

To gain some insight as to why we might expect to see late time acceleration, consider the slow roll regime of the dpion field, which is obtained by neglecting the acceleration term. In this approximation, when the dpion field exhibits an inverse scaling with time, *π* = *π*<sup>0</sup> *a*−1(*t*), the equation of motion reduces to

$$3H\frac{d}{a^2} = \frac{\lambda}{a} \left(\frac{\pi\_0^2}{a^2} - f\_D^2\right) \,. \tag{12}$$

*Universe* **2020**, *6*, 75

Solving this latter results in a power law acceleration of the Universe, namely *H*(*t*) *t* <sup>−</sup>1, provided that2 *π*(*t*0) = *π*<sup>0</sup> >> *fD* and, as customary when taking into account cosmological scalar fields, the slow roll condition holds: 3*Hπ*˙ *V* >>*π*¨. When the interaction term Lint =−g Œ(t)f 2 <sup>D</sup> ß(t) between the dpion and the gauge field is considered, we will see that this term persists to have a nearly constant energy density yielding a negative pressure equation of state. Finally, it is straightforward to show that a slightly different behavior in the time dependence of the dpion, *i.e. π* = *π*<sup>0</sup> *a*−*n*(*t*) with *n* > 0, would yield the same late time-behavior *H*(*t*) *t* <sup>−</sup>1.

Late time acceleration is recovered when the gauge field asymptotically evolves in time as the scale factor, and the pion field approaches the constant value *π fD*. Below we show how these self consistent solutions to the equations of motion for the pion and gauge field can be recovered, by working in comoving coordinates and assuming the expansion of the universe to be given by de Sitter phase. Finally, given the non-linearities in the coupled differential equations, we pursue a numerical analysis of the full system of equations.

First we consider the energy density *ρ* and the pressure *P* for our low energy effective system that emerges from the dark sector, namely

$$
\rho = \rho\_{\rm YM} + \rho\_{\rm Af} \, \, \, \qquad P = \frac{1}{3} \rho\_{\rm YM} - \rho\_{\rm Af} \, \, \, \, \tag{13}
$$

where

$$\rho\_{\rm YM} = \frac{3}{2} (\dot{\phi} + H\phi)^2 + \frac{3}{2} \mathbf{g}^2 \phi^4, \qquad \rho\_{\rm Al} = 3 \,\mathrm{g} \,\phi \,\bar{f}(a) \,, \tag{14}$$

and recall that the Friedmann equations are given by

$$\begin{aligned} H^2 \, \mathcal{M}\_p^2 &= \frac{1}{2} (\dot{\phi} + H\phi)^2 + \frac{1}{2} g^2 \phi^4 + g \, \phi f(a) \\ &+ \frac{1}{6} \dot{\pi}^2 + \frac{\lambda}{12} \left(\pi^2 - f\_D^2\right)^2, \\\ (\dot{H} + H^2) \, \mathcal{M}\_p^2 &= -\frac{1}{2} (\dot{\phi} + H\phi)^2 - \frac{1}{2} g^2 \, \phi^4 + g \, \phi \, \overline{f}(a) \\ &- \frac{1}{3} \dot{\pi}^2 + \frac{\lambda}{12} (\pi^2 - f\_D^2)^2. \end{aligned} \tag{15}$$

Having derived the differential equations that govern the dynamical system, we can now proceed to solve it.

#### **4. Field Dynamics**

Unlike usual gauge field theories, where the gauge fields dilute during cosmic expansion, the coupling of the gauge field to the dquark current leads to a growth of its homogenous component. This can be understood from the inspection of the equations of motion on the FLRW background. The growth of the gauge field will generically occur as it scales with *a*(*t*), while the interaction energy *ρA J* =*g φ* ¯*J*(*a*) remains nearly constant at late times. We can then get ahead with our purpose of solving the full dynamical system for the fields involved, so as to plot the evolution of barotropic index

<sup>2</sup> Notice that assuming a very large initial value of the pion field may actually induce to consider the tower of all the higher dimensional operators. Nevertheless, taking into account these terms would just strengthen our argument, since the crucial ingredient that ensures the quasi-de Sitter solution at the cosmological level is the mass of the pion. Thus any other higher order term would imply a power law time evolution at earlier cosmological times. Even more, at very high energy the very same concept of pion ground-state looses meaning, and the asymptotic behaviour for the pion field that is sustaining the accelerated phase of the Universe shall no longer be considered. Bearing in mind that our intent in this section is only to show the power law behavior of the background at earlier cosmological times, we avoid further comments.

*Universe* **2020**, *6*, 75

$$\begin{split} w &= \frac{P}{\rho} = \frac{P\_{\mathbf{A}\mathbf{j}} + P\_{\mathbf{Y}\mathbf{M}} + P\_{\pi}}{\rho\_{\mathbf{A}\mathbf{j}} + \rho\_{\mathbf{Y}\mathbf{M}} + \rho\_{\pi}} = \\ &= \frac{\frac{1}{2}(\not{p} + H\phi)^{2} + \frac{1}{2}\mathbf{g}^{2}\phi^{4} - 3\mathbf{g}\phi\overline{\mathbb{J}}(a) + \frac{1}{2}\dot{\pi}^{2} - \frac{\lambda}{4}(\pi^{2} - f\_{\text{D}}^{2})^{2}}{\frac{3}{2}(\not{p} + H\phi)^{2} + \frac{3}{2}\chi^{2}\phi^{4} + 3\chi\phi\overline{\mathbb{J}}(a) + \frac{1}{2}\dot{\pi}^{2} + \frac{\lambda}{4}(\pi^{2} - f\_{\text{D}}^{2})^{2}}. \end{split} \tag{16}$$

Under customary assumption, we are able to solve for the coupled system of differential equations in the configuration space {*φ*(*t*), *π*(*t*)}, and to find solutions consistent with a de Sitter expanding phase. We assume in Equation (11) the slow roll condition for *π*. Furthermore, we assume that at any time energy densities are dominated by the coupling of the gauge field to the axial current, *ρπ* <<*ρ***AJ** and *ρ***YM** << *ρ***AJ**, and show later that these assumptions are consistent with the solutions obtained. In this heuristic analysis, the dynamical system is analytically solved imposing initial conditions at recombination3. From Equations (11) and (15) we find, imposing slow roll conditions on *π*, and then considering the late-time dominant contribution to the solution,

$$
\pi(t) \simeq f\_D \left[ 1 - c\_0 \exp\left( -\frac{2\lambda f\_D^2 t}{3H} \right) \right]^{-\frac{1}{2}},\tag{17}
$$

which asymptotically reaches the value *π*<sup>∞</sup> = *fD*, and in which *c*<sup>0</sup> is an initial constant. We can then solve for the gauge field, and within a similar assumption on its time derivatives than the condition imposed for *π*(*t*) we find

$$\Phi(t) \simeq \frac{f\_D}{(2g)^{\frac{1}{\beta}}} \left[ 1 - c\_0 \exp\left( -\frac{2\lambda f\_D^2 t}{3H} \right) \right]^{-\frac{1}{\theta}}.\tag{18}$$

Both solutions Equayion (17) and (18) are monotonically decreasing and converge asymptotically towards values that are proportional to *fD*; thus their product conspire to provide an accelerating solution well approximated by a de Sitter phase, the effective cosmological constant of which assumes the asymptotic value

$$
\Lambda \simeq \frac{f\_D^4}{M\_p^2} \,. \tag{19}
$$

Supernovae data, which entail at current times *<sup>H</sup>* <sup>10</sup>−<sup>42</sup> GeV, are consistent with the asymptotic value for the gauge field *<sup>A</sup> fD* <sup>10</sup>−<sup>3</sup> eV, the coupling constant *<sup>g</sup>* having been assumed to be order unity. This suggests a fascinating conclusion: cosmic acceleration is the consequence of the CSB in the dark sector, since it occurs at the same scale of energy, *i.e. MDE fD*. We remark that our heuristic analysis has ben focusing only on the radiation energy density of the dark SU(3) sector, with the aim of showing how the evolution of the dark energy component — namely, the interaction energy density among the pion condensate field and the dark radiation — would differ from that one of radiation, rather than dealing with the observed amount of (visible and dark) matter of the Universe.

Conservation of the energy-momentum tensor is easily checked, and the attractor behavior of the solutions is also recovered. Indeed, specifying the numerical values *g* = *λ* = 0.1, the system of differential equations provides the unique fixed point:

$$(\Phi\_f, \pi\_f, H\_f) = (2.17 \cdot 10^{-3} \text{eV}, 2.04 \cdot 10^{-3} \text{eV}, 1.08 \cdot 10^{-35} \text{eV}).$$

<sup>3</sup> Integration over the dark quark degrees of freedom may provide higher-order derivative terms in the dark gauge sector, exactly as it happens for the Euler-Heisenberg effective action in QED. Higher order derivative terms then provide corrections to the radiation energy-density redshift dependence, which are subdominant in the asymptotic time limit. In our case, these corrections will be further suppressed by powers of the masses of the heaviest dark quarks present in our picture and powers of the coupling constants.

Upon linearization, the first order dynamical system recast in term of derivatives in the variable *N* = ln *a*(*t*), provides a matrix the eigenvalues of which have negative real components. This ensures that the fixed point is an attractor, and that the system asymptotically converges towards a de Sitter phase. A detailed analysis of the linearization of the first order dynamical system, and of the related eigenvalues of the matrix that is hence recovered, is provided in Ref. [22].

A more cogent numerical analysis corroborates the analytical investigations reported above. Related plots are shown in Figure 1, which makes evident the transition from a radiation dominated epoch (at the recombination *z*1100) to our present time (*z* = 0) dominated by DE. Coupling with Standard Model matter would affect the time-evolution of *w* at *z*>>10, without changing its qualitative behaviour and the asymptotic value at recombination epoch.

**Figure 1.** Plot of *w* against the redshift *z*, from current time up to recombination. For the blue line, the coupled system of non-linear differential Equations (10), (11) and (15) has been solved numerically for *<sup>g</sup>* <sup>=</sup>10−<sup>1</sup> and *<sup>λ</sup>* <sup>=</sup>1, and the initial conditions on the functions *<sup>H</sup>*(0) =10−<sup>42</sup> GeV and *<sup>φ</sup>*(0) *<sup>π</sup>*(0) *fD* =10−<sup>3</sup> eV , and on their derivatives *φ* (0)*π* (0) =10−<sup>5</sup> (eV)<sup>2</sup> . Transition from DE to (dark sector) radiation happens for *z* 2 in the blue line. The plotted red line entails a coupling constant *g*, the value of which is one order of magnitude smaller than in the case of the blue line.

Notice that in our toy-model, which provides a mechanism for an effective dynamical cosmological constant, the strong energy condition is clearly violated Ref. [24]. Indeed, the (predominant, with respect to the other contributions) interaction energy density *ρ***AJ** equals the opposite sign of the (predominant, with respect to the other contributions) interaction pressure density *p***AJ**. Nonetheless, the null energy condition is not violated, the interaction energy density providing a positive contribution to the total energy density.

Finally, we wish to emphasize the relevance of the value of the coupling constant *g*, as already outlined in Figure 1. At this purpose, considering the renormalization group flow of the coupling constant *g* will provide improvement of the toy-model we present here, and shall be considered in forthcoming more refined and realistic analyses. On a different foot, the relevance of the running of the coupling constant was also studied for QCD in Ref. [15]. Nonetheless, differently that in this latter investigation, we do not expect here that novelties in the attempt to reproduce the cosmological dynamics of the accelerating Universe will be provided by the emergence of a discrete symmetry, but rather that the absolute value of the coupling constant will affect the precise determination of the transition epoch to the dark energy cosmological behavior.

#### **5. Perturbation Analysis**

The analysis of perturbations can be developed on the FLRW background by implementing the choice of conformal coordinates *gμν* = *a*(*η*)<sup>2</sup> *ημν*. Varying action Equation (7), we recover the gauge field equation of motion:

$$
\eta^{\mu\rho} (\partial\_{\mu} F^{a}\_{\rho\nu} + \g. \epsilon^{a}\_{bc} \, A^{b}\_{\mu} F^{c}\_{\rho\nu}) = \g. a^{4} \partial^{5}\_{\nu} a
$$

From the above equation, the gauge field spatial perturbations immediately follow:

$$\begin{split} \left(\Box \delta A^{a}\_{i} + \operatorname{g}\left(\epsilon^{a}\_{\operatorname{bc}} A^{b}\_{\slash} \partial^{j} \delta A^{c}\_{i} + \epsilon^{a}\_{\operatorname{bc}} \delta A^{b}\_{\slash} \partial^{j} A^{c}\_{i}\right) + \\ + \operatorname{g}^{2}\Big(\delta A^{k}\_{b} A^{b}\_{[k} A^{a}\_{i]} + A^{k}\_{b} \delta A^{b}\_{[k} A^{a}\_{i]} + A^{k}\_{b} A^{b}\_{[k} \delta A^{a}\_{i]}\Big) &= \\ = \operatorname{g} f\_{D} a(\eta)^{4} \,\,\partial\_{i} \delta \pi^{a}(\vec{x}, \eta) \,. \end{split} \tag{20}$$

The time component perturbations (denoted by = *d*/*dη*) of the gauge field are found to be:

$$
\Box \delta A\_0^a + \mathcal{g} \epsilon^a{}\_{bc} A\_k^b \, \partial^k \delta A\_0^c - \mathcal{g}^2 A\_b^i \delta A\_{\left[0\right.}^b A\_{\left[\right.]}^a = \mathcal{g} \, a^4 \, \delta \pi^{a\prime}(\vec{x}, \eta) \, .
$$

As stated above while using co-moving coordinates, the background solutions for *A<sup>a</sup> <sup>j</sup>*(*η*) are subjected to the gauge *A<sup>a</sup> <sup>j</sup>*(*η*) = *<sup>A</sup>*(*η*) *<sup>δ</sup><sup>a</sup> <sup>j</sup>* , which reduces the spatial perturbation equation to

$$\begin{split} &\Box \delta A^{a}\_{i}(\vec{\text{x}},\eta) + \text{g}\, A(\eta) \left[\nabla \wedge \delta \vec{A}(\vec{\text{x}},\eta)\right]^{a}\_{i} + \\ &+ \text{g}^{2} \left[2A^{2}(\eta)\,\text{Tr}[\delta A^{a}\_{i}(\vec{\text{x}},\eta)]\,\delta^{a}\_{i}\right] = \text{g}\, a(\eta)^{4}\,\partial\_{i}\delta\pi^{a}(\vec{\text{x}},\eta)\,, \\ &\Box \delta A^{a}\_{0}(\vec{\text{x}},\eta) + \text{g}\, A(\eta)\,\epsilon^{a}\_{\text{bc}}\,\partial^{b}\delta A^{c}\_{0}(\vec{\text{x}},\eta) + \\ &+ 2\text{g}^{2}A(\eta)^{2}\delta A^{a}\_{0}(\vec{\text{x}},\eta) = \text{g}\, a(\eta)^{4}\,\delta\pi^{a'}(\vec{\text{x}},\eta)\,, \end{split} \tag{22}$$

Taking the trace of the gauge-field perturbations and writing *δA<sup>a</sup> <sup>i</sup>* = Tr[*δA<sup>a</sup> <sup>i</sup>* (*x*, *<sup>η</sup>*)] *<sup>δ</sup><sup>a</sup> <sup>i</sup>* /3 yield

$$\begin{split} \left[ \square \operatorname{Tr} [\delta A\_i^a(\vec{x}, \eta)] \right] &+ 6g^2 A^2(\eta) \operatorname{Tr} [\delta A\_i^a(\vec{x}, \eta)] \\ &= g \, a(\eta)^4 \, \operatorname{Tr} [\partial\_i \delta \pi^a(\vec{x}, \eta)] \,. \end{split} \tag{23}$$

Implementing the same gauge, we recover the equation for the perturbations to the dpion field,

$$\begin{split} \frac{\Box}{a} \delta \pi^{a}(\vec{\text{x}}, \eta) + 2 \frac{g'}{a^{2}} \delta \pi^{a'}(\vec{\text{x}}, \eta) + \lambda a \, \delta \pi^{a}(\vec{\text{x}}, \eta) \left( 3 \, \pi(\eta)^{2} - f\_{\text{D}}^{2} \right) \\ = -gf\_{\text{D}} \, \left( \partial^{i} \delta A\_{i}^{a}(\vec{\text{x}}, \eta) + \delta A\_{0}^{a'}(\vec{\text{x}}, \eta) \right) \, . \end{split} \tag{24}$$

We are interested in seeing if sub-horizon modes of the dpion field, (24), develop instabilities and cluster. Therefore we focus on wavelengths which are either sub-horizon or above the binding energy of the dpions involved, *H*>>| *k*|>>*fD* , by studying the evolution of plane-waves *δπ*(*x*, *<sup>η</sup>*)*α*(*η*) exp (*<sup>ı</sup> k*·*x*) solutions to (24). The equation of motion further reduces to *α* + 2*Ha α* + (*k*<sup>2</sup> + *λ a*<sup>2</sup> *f* <sup>2</sup> *<sup>D</sup>*) 0, solutions of which are superpositions of spherical Bessel functions of the first and second type,

$$\begin{aligned} \alpha(\eta) &= \alpha\_1 \, \dot{y}\_\nu \left( \frac{k}{Ha(\eta)} \right) + \alpha\_2 \, \, \dot{y}\_\nu \left( \frac{k}{Ha(\eta)} \right), \\\nu &= \, \frac{-H \pm \sqrt{H^2 - 4\lambda f\_D^2}}{2H} \end{aligned}$$

in which *k*=| *k*|, and are convergent to zero at late times for a proper choice of the initial conditions, *<sup>α</sup>*<sup>1</sup> <sup>∈</sup> <sup>R</sup> and *<sup>α</sup>*<sup>2</sup> <sup>=</sup> 0. Indeed, assuming *<sup>λ</sup>* << 1 we see that sub-horizon modes are oscillatory and bounded over all the time axis. This behavior mimics the behavior of super-horizon modes of scalar fields during inflation and is a consequence of the acceleration of space-time. Zero modes are constant and not evolving in time.

Finally, perturbations of the gauge field decrease exponentially (in comoving time), *i.e.* show the conformal time behavior

$$\operatorname{Tr}[\overleftarrow{\delta}\overleftarrow{A}\_i^{\sf a}(k,\eta)] \simeq 1/a(\eta)\,.$$

We therefore conclude that all sub-horizon perturbations are suppressed. In future work we will perform a full perturbation analysis taking into account metric perturbations Ref. [25].

#### **6. Conclusions and Discussion**

We have provided a model of late time acceleration from minimal assumptions, in that aside from instantiating a dark non-abelian copy, we have not introduced any new physics. In fact, we have employed the well known physics of the Nambu–Jona-Lasinio mechanism of chiral symmetry breaking in ordinary QCD applied to IQCD, which is well motivated from string theory. Late time acceleration emerges from the interaction between gravity, a chiral condensate and an invisible gluon that fills the universe today. A preliminary perturbation analysis shows that DE does not cluster on sub-horizon modes.

Finally, we leave to detailed investigations (see e.g., Ref. [22]) the analysis of constraints on the coupling to the visible sector, and the eventual behavior as dark-matter of dquark and dpions in this scenario. Indeed we find intriguing that our model has the possibility of connecting late time acceleration to Dark Matter (see e.g., Refs. [26,27]).

**Author Contributions:** A.A., S.A., A.M., contributed equally to this article. All authors have read and agreed to the published version of the manuscript.

**Funding:** S.A. is supported by the US Department of Energy under grant DE-SC0010386. A.M. wishes to acknowledge support by the Shanghai Municipality, through Grant No. KBH1512299, by Fudan University, through Grant No. JJH1512105, and by NSFC, through Grant No. 11875113.

**Acknowledgments:** We are indebted to David Spergel for enlightening discussions during the early stage of this work. We also thank Robert Caldwell, John Collins, Gia Dvali, Daniel Grinstein and Diego Guadagnoli for their useful feedback.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).
