Emergence and Late-Time Evolution of SU(N) Symmetric Multiplet of Pseudoscalar Fields as an Origin of Multi-Component Dark Matter

1. Introduction

The discovery of the accelerated expansion of the Universe has radically reconstructed a set of trends in the study of models of the Universe evolution [1]. One of the most obvious new trend: the introduction of some hypothetical cosmic substratum, the dark energy with negative pressure, (see, e.g., [2,3,4,5]) was supplemented by a more understandable alternative, which can be formulated as a Modified Theory of Gravity [6,7,8,9,10]. In the context of this theory one can find a lot of versions of the so-called Scalar-Tensor and Vector-Tensor theories of gravity, which are now elaborated in detail (see, e.g., [11,12,13,14,15,16]). However, there is a certain lack of research in the area, that is associated with the SU(N) symmetric models in Modified Gravity (see, e.g., [17] for references and formulation of the problems). In particular, the hypothesis that fundamental field configurations, known today and experimentally studied during last century, could have had their own SU(N) symmetric prototypes, still remains underdeveloped. This hypothesis was voiced by V.A. Rubakov in his lectures [18,19], and it means, in particular, that in addition to the canonic Vector-Tensor Gravity theory one needs to formulate, e.g., the SU(N) symmetric gauge invariant Vector-Tensor Gravity theory, as well as, the SU(N) symmetric gauge invariant (Pseudo)Scalar-Tensor Gravity theory. Why is this problem interesting in the cosmological context?

First of all, one should emphasize, that many cosmologists support the idea that the dark matter, the important element of all cosmological models, is multi-component, and one of its fractions is formed by axions, massive pseudo-Goldstone bosons [20,21,22,23,24,25,26,27,28,29,30,31]. One can assume that in the early Universe the SU(N) symmetric multiplet of pseudoscalar fields was born, and now it forms the relic multi-component dark matter. In particular, the relic light axions form the cold dark matter, while the pseudo-particles with a greater mass, known as WIMPs and ALPs, form the warm and hot fractions of the dark matter (see, e.g., [18,19,31] for explanation of terminology).

Similarly, now there exists the model indicated as dynamic aether [32,33,34,35,36,37,38,39,40,41], which was introduced as a special version of the Vector-Tensor branch of the Modified Gravity, and the aether, the new hypothetical cosmic substratum, was indicated as a candidate for the role of a dark energy. This theory introduces into consideration a timelike unit four-vector field $U^j$, which is interpreted as the velocity four-vector of the aether. Following the Rubakov hypothesis one can assume that in the early Universe there existed a SU(N) symmetric multiplet of vector fields $U^{\left(a\right)j}$, which degenerated into a singlet $U^j$ in a later era [42,43,44].

No extended SU(N) symmetric model can be operational without the introduction of the gauge field, which supports the invariance with respect to the SU(N) symmetry. This means that one needs to include into the Lagrangian the extended covariant derivatives of the vector and pseudoscalar field multiplets [45,46], thus providing the multiplet of potential co-vectors $A_m^{\left(a\right)}$, which describes the Yang-Mills field, to appear inevitably as the structural element of the theory. This gauge fields are associated with the color aether and multi-axion configuration, exclusively.

So, we have the SU(N) symmetric system, which contains three interacting multiplets: the pseudoscalar one $\left\{{\varphi }^{\left(a\right)}\right\}$, the vectorial multiplet $\left\{U^{\left(a\right)j}\right\}$ and the multiplet of the Yang-Mills potentials, $\left\{A_m^{\left(a\right)}\right\}$ (the index $\left(a\right)$ runs over $1,...,N^2-1$, indicating the fields in the adjoint representation of the SU(N) group). In other words, we consider the extension of the U(1) symmetric Einstein - Maxwell - aether - axion theory [47,48,49,50] and establish the SU(N) symmetric model, which can be indicated as Einstein - Yang - Mills - color aether - multi-axion theory.

One of the goals of the presented work is to describe the mechanism of generation of the SU(N) symmetric multiplet of pseudoscalar fields in the process of evolution of the color aether. In other words we intend to work with a multi-axionic system, which is associated now with a multi-component cosmic dark matter.

To be more precise from the mathematical point of view, we take the well-known Peccei-Quinn term ${\displaystyle \frac{1}{4}}\varphi F_{mn}^{\ast \left(a\right)}{F}^{\left(b\right)mn}{G}_{\left(a\right)\left(b\right)}$ with the pseudoscalar (axion) singlet $\varphi$ in front (see [20]), and try to extend this interaction term using the pseudoscalar multiplet $\left\{{\varphi }^{\left(c\right)}\right\}$. Clearly, if we use the product of the Yang-Mills field strength $F^{\left(b\right)mn}$ with its dual pseudotensor $F_{mn}^{\ast \left(a\right)}$, and the metric in the group space $G_{\left(a\right)\left(b\right)}$ only, we cannot provide the appropriate inclusion of the multiplier ${\varphi }^{\left(c\right)}$, which adds the third (extra) group index. We need a new supplementary color object, which has to be a vector in the group space and a scalar in the physical spacetime. For this purpose in this work we suggest to use the extended covariant divergence of the vector field ${\Omega }^{\left(a\right)}={\widehat{D}}_m{U}^{\left(a\right)m}$, which gives us the required SU(N) symmetric multiplet of scalars ${\omega }^{\left(a\right)}={\displaystyle \frac{{\Omega }^{\left(a\right)}}{\sqrt{{\Omega }^{\left(d\right)}{\Omega }^{\left(b\right)}{G}_{\left(d\right)\left(b\right)}}}}$. This new internal instrument stored in the toolbox of the color aether and previously unused, expands our modeling capabilities, since now we can use in the Lagrangian modeling the following three new SU(N) symmetric invariants: ${\varphi }^{\left(d\right)}{\omega }_{\left(d\right)}{F}_{mn}^{\ast \left(a\right)}{F}_{\left(a\right)}^{mn}$, ${\varphi }_{\left(d\right)}{\omega }^{\left(b\right)}{F}_{mn}^{\ast \left(d\right)}{F}_{\left(b\right)}^{mn}$, and ${\varphi }^{\left(b\right)}{\omega }_{\left(d\right)}{F}_{mn}^{\ast \left(d\right)}{F}_{\left(b\right)}^{mn}$. Respectively, variation of the listed terms with respect to ${\varphi }^{\left(h\right)}$ generates three new sources in the master equations of the pseudoscalar fields, thus supporting the idea that we have to deal with the axionic multiplet instead of the axionic singlet.

From the physical point of view, the following scheme of pseudoscalar multiplet emergence can be suggested. If we admit that in the early Universe the multiplet of vector fields $\left\{U^{\left(a\right)j}\right\}$ existed, the Yang-Mills field potentials $A_m^{\left(a\right)}$ were inevitably generated in order to keep the SU(N) symmetry of the color aether. Why does it turn out this way? To ensure the SU(N) invariance, it is not enough for the Lagrangian of the vector fields to contain only the covariant spacetime derivative ${\nabla }_m{U}^{\left(a\right)j}$; it has to include the extended covariant derivative ${\widehat{D}}_m{U}^{\left(a\right)j}={\nabla }_m{U}^{\left(a\right)j}+g{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_m^{\left(b\right)}{U}^{\left(c\right)j}$ based on the group constants $f_{\left(b\right)\left(c\right)}^{\left(a\right)}$ and gauge potentials $A_m^{\left(b\right)}$[45,46]. Thus, the corresponding Yang-Mills field strength $F_{mn}^{\left(a\right)}$ and its dual object $F_{mn}^{\ast \left(d\right)}$ appear in the total Lagrangian to describe properly the gauge field. If the vector fields have non-vanishing covariant divergence, i.e., ${\widehat{D}}_m{U}^{\left(a\right)m}\ne 0$, then conditions for the emergence of the pseudoscalar multiplet inevitably arise. In other words, the multi-axionic configuration can appear as a "by-product" of the color aether evolution.

In order to explain the late-time state of the dynamic aether, we consider a special ansatz, that in one of the early epochs of the Universe expansion a Phase Transition took place, as a result of which the multiplet $\left\{U^{\left(a\right)j}\right\}$ turned into the singlet $U^j$. In other words, we assume that at the moment $t_{\ast }$ of the cosmological time the SU(N) symmetric multiplet of vector fields was converted into the one-dimensional bundle of vectors parallel in the group space, i.e., $U^{\left(a\right)j}=q^{\left(a\right)}{U}^j$, where $q^{\left(a\right)}$ is a constant vector in the group space. In the works [43,44] we indicated the transition from the SU(N) symmetric aether to the canonic dynamic aether as a spontaneous polarization of the color aether. This idea was based on the analogy with Phase Transitions in ferroelectric materials under the influence of the temperature drop [51], since, passing the Curie point, the ferroelectric material leaves the symmetrical phase and enters the dissymmetrical one.

A question arises: what the Yang-Mills fields do, when the color aether converts into the canonic dynamic aether? When we were dealing with an axionic singlet and the gauge field was generated solely by a multiplet of vector fields [43,44], it was quite reasonable to assume that the multiplet of gauge fields also became one-dimensional in the group space (or quasi-Abelian), i.e., $A_m^{\left(a\right)}=Q^{\left(a\right)}{A}_m$. Such mechanism of the gauge field restructuring was discussed in many works, we would like to quote only two of them [52,53]. In the simplest case $Q^{\left(a\right)}$ and $q^{\left(a\right)}$ are parallel, i.e., $Q^{\left(a\right)}=q^{\left(a\right)}$, however, more interesting models exist, e.g., when these vectors in the group space are orthogonal, $q_{\left(a\right)}{Q}^{\left(a\right)}=0$.

The situation is becoming noticeably more complicated, when we work with the SU(N) symmetric pseudoscalar multiplet ${\varphi }^{\left(a\right)}$, which also contributes the source of the gauge field and requires that the gauge field remains non-Abelian. One can imagine at least two outcomes. The first one means that the decay of color aether is the dominating process, and it causes that the pseudoscalar multiplet also become parallel, ${\varphi }^{\left(a\right)}={\tilde{Q}}^{\left(a\right)}\varphi$, providing the parallelization of the gauge field. The second outcome is that the gauge field becomes quasi-Abelian, however the restructuring of the pseudoscalar multiplet is so masterful, that it remains non-parallel in the group space. This quite unexpected outcome is analyzed in this work.

The paper is organized as follows. In Section II we discuss the extended mathematical formalism of the presented theory, i.e., we reconstruct the total action functional, derive the master equations for the SU(N) symmetric vector, gauge and pseudoscalar fields, as well as, the corresponding equations of the gravity field. In Section III we reduce the developed formalism in the framework of Bianchi-I anisotropic homogeneous cosmological model. In Section IV we analyze the late-time behavior of the pseudoscalar multiplet. Section V contains conclusions.

2. The Formalism

2.1. Action Functional

2.1.1. SU(N) Symmetric Generalization of the Theory of the Dynamic Aether

We consider the model, which can be described by the following action functional

$$-S=\int d^{4}x\sqrt{-g}\left\{{\displaystyle \frac{1}{2\kappa }}\left[R+2\mathrm{\Lambda }+\lambda \left(U_{\left(a\right)}^m{U}_m^{\left(a\right)}-1\right)+{\mathcal{K}}_{\left(a\right)\left(b\right)}^{ijmn}{\widehat{D}}_i{U}_m^{\left(a\right)}{\widehat{D}}_j{U}_n^{\left(b\right)}\right]+\right.$$

$$\left.+{\displaystyle \frac{1}{4}}{\mathcal{C}}_{\left(a\right)\left(b\right)}^{mnls}{F}_{mn}^{\left(a\right)}{F}_{ls}^{\left(b\right)}+{\displaystyle \frac{1}{2}}{\mathrm{\Psi }}_0^2\left[V\left(\mathrm{\Phi }\right)-{\mathcal{K}}_{\left(a\right)\left(b\right)}^{sl}{\widehat{D}}_s{\varphi }^{\left(a\right)}{\widehat{D}}_l{\varphi }^{\left(b\right)}\right]\right\}.$$

(1)

Here, R is the Ricci scalar, $\mathrm{\Lambda }$ is the cosmological constant, $\kappa =8\pi G$ is the Einstein constant (we use the system of units with $c=1$). The SU(N) symmetric multiplet of vector fields $\left\{U^{\left(a\right)j}\right\}$ describes the object, which is logical to indicate as the color aether, for short; the Lagrange multiplier $\lambda$ introduces the term, which guarantees that the extended normalization condition exists

$$g_{mn}{G}_{\left(a\right)\left(b\right)}{U}^{\left(a\right)m}{U}^{\left(b\right)n}=1,$$

(2)

as the generalization of the relationship $g_{mn}{U}^m{U}^n=1$ in the canonic theory of the dynamic aether [32]. The object ${\mathcal{K}}_{\left(a\right)\left(b\right)}^{ijmn}$ is a simplest SU(N) symmetric generalization of the Jacobson constitutive tensor

$${\mathcal{K}}_{\left(a\right)\left(b\right)}^{ijmn}=G_{\left(a\right)\left(b\right)}\left[C_1{g}^{ij}{g}^{mn}+C_2{g}^{im}{g}^{jn}+C_3{g}^{in}{g}^{jm}\right]+C_4{U}_{\left(a\right)}^i{U}_{\left(b\right)}^j{g}^{mn}.$$

(3)

$G_{\left(a\right)\left(b\right)}$ is the metric in the group space, the indices $\left(a\right)$ take values $1,2,...,N^2-1$. Here we use the generalization, which does not add new coupling constants to the constitutive tensor; in fact the group indices $\left(a\right)$ and $\left(b\right)$ appear either due to the metric coefficients $G_{\left(a\right)\left(b\right)}$, or due to replacement of the canonic term $U^i{U}^j$ by the color one, $U_{\left(a\right)}^i{U}_{\left(b\right)}^j$. For discussion of the full-format extension of the constitutive tensor ${\mathcal{K}}_{\left(a\right)\left(b\right)}^{ijmn}$ see, e.g., the work [42].

2.1.2. The Yang-Mills Field Associated with the Color Aether

The term ${\widehat{D}}_m{U}_n^{\left(a\right)}$ denotes the covariant derivative of the vector fields

$${\widehat{D}}_m{U}_n^{\left(a\right)}={\nabla }_m{U}_n^{\left(a\right)}+g{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_m^{\left(b\right)}{U}_n^{\left(c\right)},$$

(4)

where ${\nabla }_m$ is the spacetime covariant derivative, g is the coupling constant associated with the potentials $A_j^{\left(c\right)}$ of the Yang-Mills field. $f_{\left(b\right)\left(c\right)}^{\left(a\right)}$ is the set of the SU(N) group constants, which are skew-symmetric with respect to the indices $\left(b\right)$ and $\left(c\right)$; as well, the set of group constants $f_{\left(d\right)\left(b\right)\left(c\right)}=G_{\left(d\right)\left(a\right)}{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}$ is skew-symmetric with respect to the indices $\left(b\right)$ and $\left(d\right)$. We keep in mind that the metric in the group space $G_{\left(a\right)\left(b\right)}$ and the structure constants $f_{\left(a\right)\left(c\right)}^{\left(d\right)}$ are the gauge covariant constant tensors, i.e.,

$${\widehat{D}}_m{G}_{\left(a\right)\left(b\right)}=0,{\widehat{D}}_m{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}=0.$$

(5)

The tensor $F_{mn}^{\left(a\right)}$ is connected with the Yang-Mills potentials $A_i^{\left(a\right)}$ as follows:

$$F_{mn}^{\left(a\right)}={\nabla }_m{A}_n^{\left(a\right)}-{\nabla }_n{A}_m^{\left(a\right)}+g{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_m^{\left(b\right)}{A}_n^{\left(c\right)}.$$

(6)

The dual tensor $F_{\left(a\right)}^{\ast ik}={\displaystyle \frac{1}{2}}{ϵ}^{ikpq}{F}_{\left(a\right)pq}$ satisfies the relations

$${\widehat{D}}_k{F}_{\left(a\right)}^{\ast ik}=0.$$

(7)

Here ${ϵ}^{ijmn}={\displaystyle \frac{E^{ijmn}}{\sqrt{-g}}}$ is the Levi-Civita pseudotensor; $E^{ijmn}$ is the completely anti-symmetric symbol with $E^{0123}=1$. The tensor ${\mathcal{C}}_{\left(a\right)\left(b\right)}^{mnls}$ is the SU(N) symmetric extension of the Tamm tensor, which appeared in the U(1) symmetric electrodynamics [54]. In this work we use the following generalization of the Tamm constitutive tensor:

$${\mathcal{C}}_{\left(a\right)\left(b\right)}^{ijmn}={\displaystyle \frac{1}{2}}{G}_{\left(a\right)\left(b\right)}\left(g^{im}{g}^{jn}-g^{in}{g}^{jm}\right)+{\displaystyle \frac{1}{2}}{ϵ}^{ijmn}\left\{G_{\left(a\right)\left(b\right)}{\omega }_{\left(d\right)}{\varphi }^{\left(d\right)}+\beta \left[{\omega }_{\left(a\right)}{\varphi }_{\left(b\right)}+{\omega }_{\left(b\right)}{\varphi }_{\left(a\right)}\right]\right\}.$$

(8)

$\beta$ is the phenomenologically introduced coupling parameter.

The new element of modeling, the multiplet of scalars ${\omega }^{\left(a\right)}$, associated with the color aether, is introduced as follows. First, we define the multiplet of scalar fields ${\Omega }^{\left(a\right)}$, which forms the color vector in the group space

$${\Omega }^{\left(a\right)}\equiv {\widehat{D}}_m{U}^{\left(a\right)m}={\nabla }_m{U}^{\left(a\right)m}+g{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_m^{\left(b\right)}{U}^{\left(c\right)m}.$$

(9)

Then we introduce the scalar $\Omega$ and the normalized vector in the group space ${\omega }^{\left(a\right)}$

$$\Omega =\sqrt{{\Omega }^{\left(a\right)}{\Omega }_{\left(a\right)}},{\omega }^{\left(a\right)}={\displaystyle \frac{{\Omega }^{\left(a\right)}}{\Omega }}.$$

(10)

We have to emphasize that the new constructive element of the extension of the Yang-Mills equations, ${\omega }^{\left(a\right)}$, appears exclusively due to the SU(N) symmetry of the color aether, and just this element plays the creative role in the emergence of the multiplet of the pseudoscalar fields, which forms the multi-component dark matter in the late-time Universe.

2.1.3. The Multiplet of Pseudoscalar Fields

The symbol ${\varphi }^{\left(a\right)}$ denotes the multiplet of pseudoscalar fields. The covariant derivative of the pseudoscalar fields is standardly defined as

$${\widehat{D}}_k{\varphi }^{\left(a\right)}={\nabla }_k{\varphi }^{\left(a\right)}+g{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_k^{\left(b\right)}{\varphi }^{\left(c\right)}.$$

(11)

The constitutive tensor ${\mathcal{K}}_{\left(a\right)\left(b\right)}^{sl}$, which form the kinetic part of the Lagrangian of the pseudoscalar field multiplet, is chosen to be of the form

$${\mathcal{K}}_{\left(a\right)\left(b\right)}^{sl}=g^{sl}{G}_{\left(a\right)\left(b\right)}-G_{\left(a\right)\left(b\right)}\left[G^{\left(c\right)\left(d\right)}-{\omega }^{\left(c\right)}{\omega }^{\left(d\right)}\right]U_{\left(c\right)}^s{U}_{\left(d\right)}^l+\gamma U_{\left(a\right)}^m{U}_{\left(b\right)m}\left[g^{sl}-U_{\left(h\right)}^s{U}^{\left(h\right)l}\right].$$

(12)

This constitutive tensor contains one new coupling constant $\gamma$.

In this work the potential $V\left(\mathrm{\Phi }\right)$ is considered to be the function of the argument $\mathrm{\Phi }=\sqrt{{\varphi }^{\left(a\right)}{\varphi }_{\left(a\right)}}$, which depends on the pseudoscalar fields only. For modeling of physical results, below we consider three versions of this potential. The first one is the standard quadratic potential $V\left(\mathrm{\Phi }\right)=m_A^2{\mathrm{\Phi }}^2$. In the second version we deal with the periodic potential

$$V\left(\mathrm{\Phi }\right)=2m_A^2\left(1-cos\mathrm{\Phi }\right),$$

(13)

which converts into the standard quadratic one $V\left(\mathrm{\Phi }\right)\to m_A^2{\mathrm{\Phi }}^2$ for small values $\mathrm{\Phi }$. The parameter $m_A$ is associated with a mass of the axion, but we do not discuss its value in the context of this work. The third version of the potential is associated with the Higgs potential

$$V\left(\mathrm{\Phi }\right)={\displaystyle \frac{1}{2}}\rho \left[{\left({\mathrm{\Phi }}^2-{\mathrm{\Phi }}_{\ast }^2\right)}^2-{\mathrm{\Phi }}_{\ast }^4\right],$$

(14)

which is known to have one unstable maximum at $\mathrm{\Phi }=0$ and two stable minima at $\mathrm{\Phi }=\pm {\mathrm{\Phi }}_{\ast }$. The constant ${\mathrm{\Psi }}_0$ in (1) is reciprocal to the coupling constant of the axion-gluon interactions $g_{AG}$, i.e., $g_{AG}={\displaystyle \frac{1}{{\mathrm{\Psi }}_0}}$.

2.2. Master Equations

2.2.1. Master Equations of Evolution of the Vector Fields

Variation of the total action functional with respect to the Lagrange multiplier $\lambda$ gives the normalization condition (2). Variation with respect to the vector fields $U^{\left(a\right)j}$ gives the balance equations

$${\nabla }_k{\mathcal{J}}_{\left(a\right)}^{kj}=g{f}_{\left(b\right)\left(a\right)}^{\left(d\right)}{A}_k^{\left(b\right)}{\mathcal{J}}_{\left(d\right)}^{kj}+\lambda U_{\left(a\right)}^j+{\mathcal{I}}_{\left(a\right)}^j,$$

(15)

where the color tensor ${\mathcal{J}}_{\left(a\right)}^{kj}$ and the Lagrange multiplier $\lambda$ are introduced in analogy with the canonic aether theory [32]:

$${\mathcal{J}}_{\left(a\right)}^{kj}={\mathcal{K}}_{\left(a\right)\left(b\right)}^{kmjn}\left[{\nabla }_m{U}_n^{\left(b\right)}+g{f}_{\left(c\right)\left(d\right)}^{\left(b\right)}{A}_m^{\left(c\right)}{U}_n^{\left(d\right)}\right],$$

(16)

$$\lambda =U_j^{\left(a\right)}\left[{\nabla }_k{\mathcal{J}}_{\left(a\right)}^{kj}-g{f}_{\left(b\right)\left(a\right)}^{\left(d\right)}{A}_k^{\left(b\right)}{\mathcal{J}}_{\left(d\right)}^{kj}\right]-U_j^{\left(a\right)}{\mathcal{I}}_{\left(a\right)}^j.$$

(17)

As for the source four-vectors ${\mathcal{I}}_{\left(a\right)}^j$, they can be written as a sum of three terms

$${\mathcal{I}}_{\left(a\right)}^j={\mathcal{I}}_{\left(a\right)\left(\mathrm{A}\right)}^j+{\mathcal{I}}_{\left(a\right)\left(YM\right)}^j+{\mathcal{I}}_{\left(a\right)\left(\mathrm{PS}\right)}^j.$$

(18)

The first term

$${\mathcal{I}}_{\left(a\right)\left(\mathrm{A}\right)}^j=C_4\left[{\widehat{D}}^j{U}_{\left(a\right)}^m{U}_{\left(b\right)}^n{\widehat{D}}_n{U}_m^{\left(b\right)}+{\widehat{D}}^k{U}_{\left(b\right)}^j{U}_n^{\left(b\right)}{\widehat{D}}_k{U}_{\left(a\right)}^n\right]$$

(19)

is the direct SU(N) generalization of the canonic aetheric current in [32]. The second term

$${\mathcal{I}}_{\left(a\right)\left(\mathrm{YM}\right)}^j=-{\displaystyle \frac{1}{4}}\kappa {\widehat{D}}^j\left\{{\displaystyle \frac{1}{\Omega }}\left[{\delta }_{\left(a\right)}^{\left(c\right)}-{\omega }^{\left(c\right)}{\omega }_{\left(a\right)}\right]\times \right.$$

(20)

$$\left.\times \left[\beta {\varphi }^{\left(h\right)}\left(F_{mn\left(c\right)}{F}_{\left(h\right)}^{mn\ast }+F_{mn\left(h\right)}{F}_{\left(c\right)}^{mn\ast }\right)+{\varphi }_{\left(c\right)}{F}_{mn}^{\left(d\right)}{F}_{\left(d\right)}^{mn\ast }\right]\right\}$$

is quadratic in the Yang-Mills field strength. The third term

$${\mathcal{I}}_{\left(a\right)\left(\mathrm{PS}\right)}^j=\kappa {\mathrm{\Psi }}_0^2{U}_{\left(a\right)}^s\left[G_{\left(d\right)\left(b\right)}+\gamma U_{\left(d\right)}^m{U}_{\left(b\right)m}\right]\left({\widehat{D}}_s{\varphi }^{\left(d\right)}\right)\left({\widehat{D}}^j{\varphi }^{\left(b\right)}\right)-$$

(21)

$$-\kappa {\mathrm{\Psi }}_0^2\gamma \left[U_{\left(b\right)}^j\left(g^{sl}-U_{\left(h\right)}^s{U}^{l\left(h\right)}\right){\widehat{D}}_l{\varphi }^{\left(b\right)}{\widehat{D}}_s{\varphi }_{\left(a\right)}\right]+$$

$$-\kappa {\mathrm{\Psi }}_0^2{U}_{\left(c\right)}^s{\omega }^{\left(c\right)}{\omega }_{\left(a\right)}{\widehat{D}}^j{\varphi }^{\left(d\right)}{\widehat{D}}_s{\varphi }_{\left(d\right)}-$$

$$+\kappa {\mathrm{\Psi }}_0^2{\widehat{D}}^j\left\{{\displaystyle \frac{1}{\Omega }}\left[{\delta }_{\left(a\right)}^{\left(c\right)}-{\omega }^{\left(c\right)}{\omega }_{\left(a\right)}\right]U_{\left(c\right)}^s{U}_{\left(d\right)}^l{\omega }^{\left(d\right)}{\widehat{D}}_s{\varphi }_{\left(b\right)}{\widehat{D}}_l{\varphi }^{\left(b\right)}\right\}$$

contains quadratic combinations of the covariant derivatives of the pseudoscalar fields.

2.2.2. Master Equations for the Yang-Mills Fields

Variation of the total action functional (1) with respect to the potentials $A_i^{\left(a\right)}$ gives the extended Yang-Mills equations

$${\nabla }_k\left[{\mathcal{C}}_{\left(a\right)\left(b\right)}^{ikls}{F}_{ls}^{\left(b\right)}\right]-f_{\left(c\right)\left(a\right)}^{\left(d\right)}{A}_k^{\left(c\right)}\left[{\mathcal{C}}_{\left(d\right)\left(b\right)}^{ikls}{F}_{ls}^{\left(b\right)}\right]={\mathrm{\Gamma }}_{\left(a\right)\left(\mathrm{U}\right)}^i+{\mathrm{\Gamma }}_{\left(a\right)\left(\mathrm{PS}\right)}^i.$$

(22)

The first part of the color current

$${\mathrm{\Gamma }}_{\left(a\right)\left(\mathrm{U}\right)}^i={\displaystyle \frac{g}{\kappa }}{f}_{\left(c\right)\left(a\right)}^{\left(d\right)}{U}_j^{\left(c\right)}{\mathcal{K}}_{\left(d\right)\left(b\right)}^{imjn}{\widehat{D}}_m{U}_n^{\left(b\right)}+$$

(23)

$$+{\displaystyle \frac{g}{4\Omega }}{f}_{\left(c\right)\left(a\right)}^{\left(d\right)}{U}^{i\left(c\right)}\left({\delta }_{\left(d\right)}^{\left(q\right)}-{\omega }^{\left(q\right)}{\omega }_{\left(d\right)}\right)\times$$

$$\times \left\{F_{mn}^{\ast \left(p\right)}{F}^{mn\left(b\right)}{\varphi }^{\left(f\right)}\left[G_{\left(p\right)\left(b\right)}{G}_{\left(q\right)\left(f\right)}+\beta \left(G_{\left(p\right)\left(q\right)}{G}_{\left(b\right)\left(f\right)}+G_{\left(p\right)\left(f\right)}{G}_{\left(b\right)\left(q\right)}\right)\right]\right\}$$

is the contribution of the multiplet of vector fields; it appears since the covariant derivative ${\widehat{D}}_k{U}^{\left(a\right)j}$ contains the potential $A_j^{\left(b\right)}$. The second part

$${\mathrm{\Gamma }}_{\left(a\right)\left(\mathrm{PS}\right)}^i=-g{f}_{\left(c\right)\left(a\right)}^{\left(d\right)}{\mathrm{\Psi }}_0^2{\varphi }^{\left(c\right)}{\mathcal{K}}_{\left(d\right)\left(b\right)}^{im}{\widehat{D}}_m{\varphi }^{\left(b\right)}+$$

(24)

$$-{\mathrm{\Psi }}_0^2{\displaystyle \frac{g}{\Omega }}{f}_{\left(c\right)\left(a\right)}^{\left(d\right)}{U}^{i\left(c\right)}\left({\delta }_{\left(d\right)}^{\left(q\right)}-{\omega }^{\left(q\right)}{\omega }_{\left(d\right)}\right)U_{\left(q\right)}^s{\widehat{D}}_s{\varphi }_{\left(b\right)}{\widehat{D}}_l{\varphi }^{\left(b\right)}{U}_{\left(h\right)}^l{\omega }^{\left(h\right)}$$

is, respectively, the contribution of the multiplet of pseudoscalar fields. The structure of the source terms in the right-hand sides of the equations (22) shows that both the multiplets of vector and pseudoscal fields form the Yang-Mills field of the SU(N) symmetric system under consideration.

2.2.3. Master Equations for the Pseudoscalar Fields ${\varphi }^{\left(a\right)}$

Variation of the total action functional with respect to the pseudoscalar fields yields $N^2-1$ equations, which can be written in the following form:

$${\nabla }_j\left\{{\mathcal{K}}_{\left(h\right)\left(b\right)}^{jl}\left[{\nabla }_l{\varphi }^{\left(b\right)}+g{f}_{\left(e\right)\left(f\right)}^{\left(b\right)}{A}_l^{\left(e\right)}{\varphi }^{\left(f\right)}\right]\right\}-g{f}_{\left(c\right)\left(h\right)}^{\left(d\right)}{A}_j^{\left(c\right)}\left\{{\mathcal{K}}_{\left(d\right)\left(b\right)}^{jl}\left[{\nabla }_l{\varphi }^{\left(b\right)}+g{f}_{\left(e\right)\left(f\right)}^{\left(b\right)}{A}_l^{\left(e\right)}{\varphi }^{\left(f\right)}\right]\right\}+$$

$$+{\displaystyle \frac{{\varphi }_{\left(h\right)}}{2\mathrm{\Phi }}}{V}^{\prime }\left(\mathrm{\Phi }\right)=-{\displaystyle \frac{1}{4{\mathrm{\Psi }}_0^2}}{F}_{mn}^{\left(a\right)\ast }\left[F_{\left(a\right)}^{mn}{\omega }_{\left(h\right)}+2\beta F_{\left(h\right)}^{mn}{\omega }_{\left(a\right)}\right].$$

(25)

The source term in the right-hand side of these equations is of zero order in the pseudoscalar field; it is provided by the Yang-Mills field strength $F_{\left(a\right)}^{mn}$ in combination with the multiplier ${\omega }_{\left(h\right)}$, originated from the covariant divergence of the vector field multiplet.

2.2.4. Master Equations for the Gravitational Field

Variation procedure with respect to metric gives the set of equations for the gravitational field

$$R_{pq}-{\displaystyle \frac{1}{2}}R{g}_{pq}=\mathrm{\Lambda }{g}_{pq}+\lambda U_p^{\left(a\right)}{U}_{\left(a\right)q}+T_{pq}^{\left(\mathrm{CA}\right)}+\kappa T_{pq}^{\left(\mathrm{YM}\right)}+\kappa T_{pq}^{\left(\mathrm{PS}\right)}.$$

(26)

The stress-energy tensor of the color vector fields is of the following form:

$$T_{pq}^{\left(\mathrm{CA}\right)}={\displaystyle \frac{1}{2}}{g}_{pq}{\mathcal{K}}_{\left(a\right)\left(b\right)}^{ijmn}{\widehat{D}}_i{U}_m^{\left(a\right)}{\widehat{D}}_j{U}_n^{\left(b\right)}+$$

(27)

$$+G_{\left(a\right)\left(b\right)}{\widehat{D}}^m\left[U_{(p}^{\left(b\right)}{\mathcal{J}}_{q)m}^{\left(a\right)}-{\mathcal{J}}_{m(p}^{\left(a\right)}{U}_{q)}^{\left(b\right)}-{\mathcal{J}}_{\left(pq\right)}^{\left(a\right)}{U}_m^{\left(b\right)}\right]+$$

$$+G_{\left(a\right)\left(b\right)}{C}_1[{\widehat{D}}_m{U}_p^{\left(a\right)}{\widehat{D}}^m{U}_q^{\left(b\right)}-{\widehat{D}}_p{U}^{m\left(a\right)}{\widehat{D}}_q{U}_m^{\left(b\right)}]+$$

$$+C_4\left\{U_{\left(a\right)}^m{\widehat{D}}_m{U}_p^{\left(a\right)}{U}_{\left(b\right)}^n{\widehat{D}}_n{U}_q^{\left(b\right)}-U_{m\left(b\right)}{U}_{n\left(a\right)}{\widehat{D}}_p{U}^{m\left(a\right)}{\widehat{D}}_q{U}^{n\left(b\right)}\right\}.$$

The parentheses $\left(pq\right)$ denote the symmetrization with respect to the coordinate indices p and q.

The symbol $T_{pq}^{\left(\mathrm{YM}\right)}$ refers to the stress-energy tensor of the Yang-Mills field

$$T_{pq}^{\left(\mathrm{YM}\right)}=G_{\left(a\right)\left(b\right)}\left[{\displaystyle \frac{1}{4}}{g}_{pq}{F}_{mn}^{\left(a\right)}{F}^{mn\left(b\right)}-F_{pn}^{\left(a\right)}{F}_q^{\left(b\right)n}\right]-$$

(28)

$$-{\displaystyle \frac{1}{4}}{g}_{pq}\left[G_{\left(a\right)\left(b\right)}{G}_{\left(c\right)\left(d\right)}+\beta \left(G_{\left(a\right)\left(c\right)}{G}_{\left(b\right)\left(d\right)}+G_{\left(a\right)\left(d\right)}{G}_{\left(b\right)\left(c\right)}\right)\right]\times$$

$$\times {\nabla }_j\left\{U^{\left(h\right)j}{F}_{mn}^{\left(a\right)}{F}^{\ast \left(b\right)mn}{\displaystyle \frac{{\varphi }^{\left(d\right)}}{\Omega }}\left[{\delta }_{\left(h\right)}^{\left(c\right)}-{\omega }^{\left(c\right)}{\omega }_{\left(h\right)}\right]\right\}.$$

The last term describes the effective stress-energy tensor of the multiplet of pseudoscalar fields

$$T_{pq}^{\left(\mathrm{PS}\right)}={\displaystyle \frac{1}{2}}{g}_{pq}{\mathrm{\Psi }}_0^2\left[V\left(\varphi \right)-{\mathcal{K}}_{\left(a\right)\left(b\right)}^{ls}{\widehat{D}}_l{\varphi }^{\left(a\right)}{\widehat{D}}_s{\varphi }^{\left(b\right)}\right]+$$

(29)

$$+{\mathrm{\Psi }}_0^2{\widehat{D}}_p{\varphi }^{\left(a\right)}{\widehat{D}}_q{\varphi }^{\left(b\right)}\left(G_{\left(a\right)\left(b\right)}+\gamma U_{\left(a\right)}^j{U}_{j\left(b\right)}\right)+$$

$$+{\mathrm{\Psi }}_0^2\gamma U_{p\left(a\right)}{U}_{q\left(b\right)}\left(g^{sl}-U_{\left(h\right)}^s{U}^{l\left(h\right)}\right){\widehat{D}}_l{\varphi }^{\left(a\right)}{\widehat{D}}_s{\varphi }^{\left(b\right)}+$$

$$+{\mathrm{\Psi }}_0^2{g}_{pq}{\nabla }_n\left\{{\displaystyle \frac{{\omega }^{\left(b\right)}}{\Omega }}{U}^{\left(d\right)n}\left[{\delta }_{\left(d\right)}^{\left(c\right)}-{\omega }^{\left(c\right)}{\omega }_{\left(d\right)}\right]U_{\left(c\right)}^s{U}_{\left(b\right)}^l{\widehat{D}}_s{\varphi }_{\left(h\right)}{\widehat{D}}_l{\varphi }^{\left(h\right)}\right\}.$$

The consequence of the Bianchi identity

$${\nabla }^q\left[R_{pq}-{\displaystyle \frac{1}{2}}{g}_{pq}R\right]\equiv 0\Rightarrow {\nabla }^q\left[\lambda U_p^{\left(a\right)}{U}_{\left(a\right)q}+T_{pq}^{\left(\mathrm{CA}\right)}+\kappa T_{pq}^{\left(\mathrm{YM}\right)}+\kappa T_{pq}^{\left(\mathrm{PS}\right)}\right]=0$$

(30)

is performed automatically on the solutions to the listed master equations.

3. Decay of the SU(N) Symmetric Configuration

The analysis of the set of coupled master equations (15), (22), (25) allows us to make the following series of remarks. First, if the multiplet of vector fields $\left\{U^{\left(a\right)j}\right\}$ was part of the field configuration in the early Universe, then the gauge field, which is characterized by the potentials $\left\{A_j^{\left(a\right)}\right\}$ and by the Yang-Mills field strength $F_{mn}^{\left(a\right)}$, was inevitably formed to support the SU(N) symmetry of the aether. Second, the multiplet of scalars ${\omega }^{\left(a\right)}$, constructed using the divergence of the vector fields ${\widehat{D}}_k{U}^{\left(a\right)k}$ (see (9) and (10)), provides the multiplet of pseudoscalar sources in the right-hand side of the equations (25) to appear. Third, if these sources are nonvanishing, the multiplet of pseudoscalar fields ${\varphi }^{\left(a\right)}$ inevitably appears as the solution to the set of equations (25). These scheme can serve as a justification for the assumption about the emergence of a multi-axion structure in the field configuration in the early Universe. If this is true, then it is worth considering the history of further evolution of the resulting pseudoscalar multiplet.

3.1. Spontaneous Polarization and Decay of Color Aether

We follow the idea that the color aether that existed at the early stage of the Universe expansion spontaneously disintegrates, and the canonical aether appears. Such an event can occur, for example, as a phase transition of the second kind, which is characterised by the alignment the vector field multiplet in the group space. Mathematically this means that all the vectors $U_{\left(a\right)}^j$ become parallel to the timelike four-vector $U^j$, i.e., $U_{\left(a\right)}^j=q_{\left(a\right)}{U}^j$, where $U^j$ is associated with the velocity four-vector of the dynamic aether, and $U^j{U}_j=1$. The set of parameters $q^{\left(a\right)}$ form a normalized constant vector in the group space, this means that $q_{\left(a\right)}{q}^{\left(a\right)}=1$. In general case, the direction of this color vector in the group space is not fixed. However, below we consider a special model with fixed direction of the color vector $q^{\left(a\right)}$, but this step will be motivated physically.

In the case of the color aether polarization the covariant derivative ${\widehat{D}}_k{U}^{\left(a\right)j}$ and the corresponding divergence take the form

$${\widehat{D}}_k{U}^{\left(a\right)j}=q^{\left(a\right)}{\nabla }_k{U}^j+g{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_k^{\left(b\right)}{U}^j{q}^{\left(c\right)},$$

(31)

$${\widehat{D}}_k{U}^{\left(a\right)k}=q^{\left(a\right)}{\nabla }_k{U}^k+g{f}_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_k^{\left(b\right)}{U}^k{q}^{\left(c\right)}.$$

(32)

We assume that the Yang-Mills potentials satisfy the so-called Landau gauge conditions $U_j{A}^{\left(a\right)j}=0$ with respect to the aether velocity four-vector. Then we obtain that ${\Omega }^{\left(a\right)}=q^{\left(a\right)}\mathrm{\Theta }$, where $\mathrm{\Theta }={\nabla }_k{U}^k$, so that ${\omega }_{\left(a\right)}=q_{\left(a\right)}$. As well, in this situation we obtain that

$$q_{\left(a\right)}{\widehat{D}}_k{U}^{\left(a\right)j}={\nabla }_k{U}^j+g{q}^{\left(a\right)}{f}_{\left(a\right)\left(b\right)\left(c\right)}{q}^{\left(c\right)}{A}_k^{\left(b\right)}{U}^j={\nabla }_k{U}^j.$$

(33)

Now we have to check whether the equations (15) remain compatible, if the spacetime platform is chosen to be of the Bianchi-I form. In other words, we consider the homogeneous anisotropic model with the spacetime metric

$$d{s}^2=d{t}^2-a^2\left(t\right){d{x}^1}^2-b^2\left(t\right){d{x}^2}^2-c^2\left(t\right){d{x}^3}^2.$$

(34)

The symmetry of the model hints that the vector, gauge and pseudoscalar fields can also be considered as spatially homogeneous thus depending on the cosmological time only. Our supplementary ansatz is that the global timelike vector field $U^j$ has to be chosen in the form $U^j={\delta }_0^j$. For this choice the covariant derivative of the aether velocity four-vector

$${\nabla }_k{U}^m={\delta }_1^m{\delta }_k^1\left({\displaystyle \frac{\dot{a}}{a}}\right)+{\delta }_2^m{\delta }_k^2\left({\displaystyle \frac{\dot{b}}{b}}\right)+{\delta }_3^m{\delta }_k^3\left({\displaystyle \frac{\dot{c}}{c}}\right)$$

(35)

is symmetric, i.e., ${\nabla }_k{U}_m={\nabla }_m{U}_k$, and can be decomposed as

$${\nabla }_k{U}^m={\sigma }_k^m+{\displaystyle \frac{1}{3}}\mathrm{\Theta }{\Delta }_k^m.$$

(36)

Here and below the dot denotes the derivative with respect to cosmological time. We use the definitions of the scalar of expansion of the aether flow

$$\mathrm{\Theta }={\nabla }_k{U}^k={\displaystyle \frac{\dot{a}}{a}}+{\displaystyle \frac{\dot{b}}{b}}+{\displaystyle \frac{\dot{c}}{c}},$$

(37)

the standard definition of the projector with respect to the aether velocity four-vector

$${\Delta }_m^k={\delta }_m^k-U^k{U}_m,$$

(38)

and of the shear tensor

$${\sigma }_{jl}={\displaystyle \frac{1}{2}}{\Delta }_j^m{\Delta }_l^n\left({\nabla }_m{U}_n+{\nabla }_n{U}_m\right)-{\displaystyle \frac{1}{3}}{\Delta }_{jl}\mathrm{\Theta }.$$

(39)

The non-vanishing components of the shear tensor are

$${\sigma }_1^1={\displaystyle \frac{\dot{a}}{a}}-{\displaystyle \frac{1}{3}}\mathrm{\Theta },{\sigma }_2^2={\displaystyle \frac{\dot{b}}{b}}-{\displaystyle \frac{1}{3}}\mathrm{\Theta },{\sigma }_3^3={\displaystyle \frac{\dot{c}}{c}}-{\displaystyle \frac{1}{3}}\mathrm{\Theta }.$$

(40)

The acceleration four-vector and the vorticity tensor are vanishing for the model with Bianchi-I spacetime platform

$$\mathcal{D}{U}^j=U^k{\nabla }_k{U}^j=0,{\omega }_{jl}={\Delta }_j^m{\Delta }_l^n\left({\nabla }_m{U}_n-{\nabla }_n{U}_m\right)=0.$$

(41)

Finally, we assume that due to the decay of the color aether the Yang-Mills field potentials also become parallel in the group space, i.e., $A_j^{\left(a\right)}=q^{\left(a\right)}{A}_j$. This means that the Yang-Mills field strength become quasi-Abelian

$$F_{mn}^{\left(a\right)}=q^{\left(a\right)}({\nabla }_m{A}_n-{\nabla }_n{A}_m),$$

(42)

and the potential four-(co)vector $A_j$ has only spatial components depending on time $A_j=(0,A_1\left(t\right),A_2\left(t\right)),A_3\left(t\right))$, since $U^j{A}_j=0$. In this case the color vector $q^{\left(a\right)}$ becomes the covariant constant one, since

$${\widehat{D}}_k{q}^{\left(a\right)}={\nabla }_k{q}^{\left(a\right)}+f_{\left(b\right)\left(c\right)}^{\left(a\right)}{A}_k{q}^{\left(b\right)}{q}^{\left(c\right)}={𝜕}_k{q}^{\left(a\right)}=0.$$

(43)

Mention should be made that the Lorentz gauge condition

$${\nabla }_k{A}^k=0\Rightarrow {\displaystyle \frac{1}{abc}}{\displaystyle \frac{d}{dt}}\left(abc{A}^0\right)=0$$

(44)

is satisfied identically, so that one can speak about the Coulomb gauge, which is the combination of the Landau and Lorentz gauge conditions. For the presented Yang-Mills field potentials the pseudoinvariant $F_{mn}^{\ast \left(a\right)}{F}_{\left(a\right)}^{mn}$ vanishes. This means, in particular, that after the decay of the color aether the sources for the multiplet of the pseudoscalar fields disappear, and we deal with the simple evolution of the multi-axion configuration born earlier.

3.2. Exact Solutions to the Master Equations for the Polarized Multiplet of Vector Fields

The ansatz that $U^{\left(a\right)j}={\delta }_0^j{q}^{\left(a\right)}$ has to be checked from the point of view of compatibility of the equations (15)-(21). One can start with the analysis of the source terms in these equations. First, we see that the reduced term ${\mathcal{I}}_{\left(a\right)\left(\mathrm{A}\right)}^j$ (19)

$${\mathcal{I}}_{\left(a\right)\left(\mathrm{A}\right)}^j=C_4\left[{\widehat{D}}^j{U}_{\left(a\right)}^m\mathcal{D}{U}_m-{\nabla }^k{U}^{j}g{f}_{\left(d\right)\left(b\right)\left(a\right)}{q}^{\left(b\right)}{A}_k{q}^{\left(d\right)}\right]$$

(45)

vanishes, since the acceleration four-vector $\mathcal{D}{U}_m$ is equal to zero, and the group constants are skew-symmetric with respect to indices $\left(d\right)$ and $\left(b\right)$. The reduced source term (20)

$${\mathcal{I}}_{\left(a\right)\left(\mathrm{YM}\right)}^j=-{\displaystyle \frac{1}{4}}\kappa {\widehat{D}}^j\left\{{\displaystyle \frac{1}{\mathrm{\Theta }}}\left[{\delta }_{\left(a\right)}^{\left(c\right)}-q^{\left(c\right)}{q}_{\left(a\right)}\right]F_{mn}{F}^{mn\ast }{\varphi }_{\left(c\right)}\right\}$$

(46)

vanishes, since now $F_{mn}{F}^{mn\ast }=0$. The source term (21) converts into

$${\mathcal{I}}_{\left(a\right)\left(\mathrm{PS}\right)}^j=\kappa {\mathrm{\Psi }}_0^2{q}_{\left(a\right)}{U}^j{\left[q_{\left(d\right)}{\dot{\varphi }}^{\left(d\right)}\right]}^2.$$

(47)

Then, keeping in mind that $C_1+C_3=0$ (see, e.g., [55] for motivation), we obtain the Jacobson tensor in the form

$${\mathcal{J}}_{\left(a\right)}^{kj}=q_{\left(a\right)}{C}_2{g}^{kj}\mathrm{\Theta },$$

(48)

and $4(N^2-1)$ master equations for the vector field happen to be transformed into

$$q_{\left(a\right)}{U}^j\left\{C_2\dot{\mathrm{\Theta }}-\lambda -\kappa {\mathrm{\Psi }}_0^2{\left[q_{\left(d\right)}{\dot{\varphi }}^{\left(d\right)}\right]}^2\right\}=0.$$

(49)

They give us only one meaningful equation, which is in fact the equation for the Lagrange multiplier

$$\lambda =C_2\dot{\mathrm{\Theta }}-\kappa {\mathrm{\Psi }}_0^2{\left[q_{\left(d\right)}{\dot{\varphi }}^{\left(d\right)}\right]}^2.$$

(50)

In other words, the extended Jacobson equations are solved and give the exact solution $U^{\left(a\right)j}=q^{\left(a\right)}{\delta }_0^j$ with $\lambda$ presented in the form (50).

3.3. Solutions to the Reduced Master Equations for the Gauge Field

The analysis of the reduced Yang-Mills equations shows that source term (23) is equal to zero for the polarized color aether and gauge fields, ${\mathrm{\Gamma }}_{\left(a\right)\left(\mathrm{U}\right)}^i=0$. The source term (24) can be written as follows

$${\mathrm{\Gamma }}_{\left(a\right)\left(\mathrm{PS}\right)}^i=-g{f}_{\left(c\right)\left(a\right)}^{\left(d\right)}{\mathrm{\Psi }}_0^2{\varphi }^{\left(c\right)}\left[U^i{\dot{\varphi }}_{\left(d\right)}+g{f}_{\left(d\right)\left(h\right)\left(b\right)}{A}^i{q}^{\left(h\right)}{\varphi }^{\left(b\right)}\right],$$

(51)

and Yang-Mills equations (22) take now the form

$$q_{\left(a\right)}{\nabla }_k{F}^{ik}=f_{\left(c\right)\left(a\right)}^{\left(d\right)}{A}_k{q}^{\left(c\right)}{F}^{\ast ik}\beta {\varphi }_{\left(d\right)}-g{f}_{\left(c\right)\left(a\right)}^{\left(d\right)}{\mathrm{\Psi }}_0^2{\varphi }^{\left(c\right)}\left[U^i{\dot{\varphi }}_{\left(d\right)}+g{f}_{\left(d\right)\left(h\right)\left(b\right)}{A}^i{q}^{\left(h\right)}{\varphi }^{\left(b\right)}\right].$$

(52)

The question arises: when this set of $4(N^2-1)$ equations is self-consistent? There are two evident examples.

1) When the pseudoscalar fields are also self-parallel in the group space, i.e., ${\varphi }^{\left(a\right)}=q^{\left(a\right)}\varphi$, we obtain immediately that

$$q_{\left(a\right)}{\nabla }_k{F}^{ik}=0\Rightarrow {\displaystyle \frac{d}{dt}}\left(abc{F}^{i0}\right)=0\Rightarrow F^{\alpha 0}\left(t\right)=F^{\alpha 0}\left(t_0\right){\displaystyle \frac{a\left(t_0\right)b\left(t_0\right)c\left(t_0\right)}{a\left(t\right)b\left(t\right)c\left(t\right)}},$$

(53)

for the color vector $q^{\left(a\right)}$ any way oriented in the group space.

2) When the non-Abelian multiplet of the pseudoscalar fields is designed so that it is not capable of forming a source for the gauge field. It is possible, first, when the second term in the right-hand side of equations (52) vanishes, i.e.,

$$f_{\left(d\right)\left(c\right)\left(a\right)}{\varphi }^{\left(c\right)}{\dot{\varphi }}^{\left(d\right)}=0.$$

(54)

This equation can be easily integrated yielding

$${\dot{\varphi }}^{\left(d\right)}=\nu \left(t\right){\varphi }^{\left(d\right)}\Rightarrow {\varphi }^{\left(d\right)}\left(t\right)={\varphi }^{\left(d\right)}\left(t_0\right)exp\left\{{\int }_{t_0}^{t}d\tau \nu \left(\tau \right)\right\}.$$

(55)

The function $\nu \left(t\right)$ has to be found from the master equations for the pseudoscalar fields.

Second, the first term in the right-hand side of equations (52) can be chosen as vanishing, since one has to take into account the symmetry of the Bianchi-I model. Indeed, in order to guarantee that the stress-energy tensor of the Yang-Mills field (28) has no non-diagonal components, we have to choose the potential vector $A_{\alpha }$ along one of the spatial axis, say, along $O{x}^3$, providing that only $F_{03}$ is non-vanishing. Then the convolution $A_k{F}^{\ast ik}={\displaystyle \frac{1}{2}}{A}_3{ϵ}^{i3mn}{F}_{mn}$ is equal to zero.

Third, the term

$$f_{\left(c\right)\left(a\right)}^{\left(d\right)}{\varphi }^{\left(c\right)}{f}_{\left(d\right)\left(h\right)\left(b\right)}{q}^{\left(h\right)}{\varphi }^{\left(b\right)}$$

(56)

can be made proportional to the color vector $q_{\left(a\right)}$ by the special choice of this vector in the group space. We can illustrate this statement for the group SU(2), for which the group constants $f_{\left(a\right)\left(b\right)\left(c\right)}$ coincide with the three-dimensional Levi-Civita symbols ${\epsilon }_{\left(a\right)\left(b\right)\left(c\right)}$. Indeed, we see that in this case

$$f^{\left(d\right)\left(c\right)\left(a\right)}{f}_{\left(d\right)\left(h\right)\left(b\right)}{q}^{\left(h\right)}{\varphi }^{\left(b\right)}{\varphi }_{\left(c\right)}=\left({\delta }_{\left(h\right)}^{\left(c\right)}{\delta }_{\left(b\right)}^{\left(a\right)}-{\delta }_{\left(b\right)}^{\left(c\right)}{\delta }_{\left(h\right)}^{\left(a\right)}\right)q^{\left(h\right)}{\varphi }^{\left(b\right)}{\varphi }_{\left(c\right)}={\varphi }^{\left(a\right)}{q}^{\left(h\right)}{\varphi }_{\left(h\right)}-q^{\left(a\right)}{\varphi }^{\left(b\right)}{\varphi }_{\left(b\right)}.$$

(57)

Thus, if the color vector $q^{\left(a\right)}$ is orthogonal to the color vector ${\varphi }^{\left(a\right)}$, i.e., $q^{\left(h\right)}{\varphi }_{\left(h\right)}=0$, we obtain the object proportional to $q^{\left(a\right)}$. In this situation the set of equations (52) can be transformed into

$$q_{\left(a\right)}\left\{{\nabla }_k{F}^{ik}-g^2{\mathrm{\Psi }}_0^2{A}^i{\varphi }^{\left(b\right)}{\varphi }_{\left(b\right)}\right\}=0,$$

(58)

and we deal in fact with only one equation

$${\displaystyle \frac{d}{dt}}\left[\left({\displaystyle \frac{ab}{c}}\right){\displaystyle \frac{d}{dt}}{A}_3\right]=-g^2{\mathrm{\Psi }}_0^2\left({\displaystyle \frac{ab}{c}}\right)A_3\left[{\varphi }^{\left(b\right)}{\varphi }_{\left(b\right)}\right]$$

(59)

for one function $A_3\left(t\right)$.

3.4. Evolutionary Equations for the Set of the Non-Abelian Pseudoscalar Fields

The equations for the pseudoscalar fields (25) loss now the source terms in their right-hand sides, and convert into the following set of equations:

$${\ddot{\varphi }}_{\left(h\right)}+\mathrm{\Theta }{\dot{\varphi }}_{\left(h\right)}+{\displaystyle \frac{g^2}{c^2}}{A}_3^2{\varphi }_{\left(h\right)}+{\displaystyle \frac{{\varphi }_{\left(h\right)}}{2\mathrm{\Phi }}}{V}^{\prime }\left(\mathrm{\Phi }\right)=0.$$

(60)

Taking into account the requirement (55) we obtain the equation for the function $\nu \left(t\right)$.

$$\dot{\nu }+{\nu }^2+\mathrm{\Theta }\nu +{\displaystyle \frac{g^2}{c^2}}{A}_3^2+{\displaystyle \frac{1}{2\mathrm{\Phi }}}{V}^{\prime }\left(\mathrm{\Phi }\right)=0.$$

(61)

The last term in this equation depends on the structure of the potential $V\left(\mathrm{\Phi }\right)$. We could mention three the most known versions of the potential modeling.

a) The simplest potential is quadratic and has the form $V\left(\mathrm{\Phi }\right)=m_A^2{\mathrm{\Phi }}^2$. In this case the last term in (61) is constant and the equation for $\nu \left(t\right)$ is the inhomogeneous Riccati differential equation

$$\dot{\nu }+{\nu }^2+\mathrm{\Theta }\nu =-{\displaystyle \frac{g^2}{c^2\left(t\right)}}{A}_3^2\left(t\right)-m_A^2.$$

(62)

b) If the potential is of the periodic form, $V\left(\varphi \right)=2m_A^2(1-cos\mathrm{\Phi })$, the equation for $\nu$ becomes the integro-differential one

$$\dot{\nu }+{\nu }^2+\mathrm{\Theta }\nu +{\displaystyle \frac{g^2}{c^2\left(t\right)}}{A}_3^2\left(t\right)+m_A^2{\displaystyle \frac{sin\mathrm{\Phi }\left(t\right)}{\mathrm{\Phi }\left(t\right)}}=0,$$

(63)

since the last term contains the function

$$\mathrm{\Phi }\left(t\right)=\mathrm{\Phi }\left(t_0\right)exp\left\{{\int }_{t_0}^{t}d\tau \nu \left(\tau \right)\right\},\mathrm{\Phi }\left(t_0\right)=\sqrt{{\varphi }^{\left(a\right)}\left(t_0\right){\varphi }_{\left(a\right)}\left(t_0\right)}.$$

(64)

c) When the potential is of the Higgs form, $V\left(\mathrm{\Phi }\right)={\displaystyle \frac{1}{2}}\rho \left[{\left({\mathrm{\Phi }}^2-{\mathrm{\Phi }}_{\ast }^2\right)}^2-{\mathrm{\Phi }}_{\ast }^4\right]$, the equation for $\nu \left(t\right)$ can be written as follows:

$$\dot{\nu }+{\nu }^2+\mathrm{\Theta }\nu +{\displaystyle \frac{g^2}{c^2}}{A}_3^2+\rho \left[{\mathrm{\Phi }}^2\left(t_0\right)exp\left\{2{\int }_{t_0}^{t}d\tau \nu \left(\tau \right)\right\}-{\mathrm{\Phi }}_{\ast }^2\right]=0.$$

(65)

In general case, $\nu \left(t\right)$ can be found only numerically, but below we consider the example of exact analytic solution.

3.5. Reduced Equations for the Gravitational Field

The obtained reduced master equations for the gravity field

$$R_q^p-{\displaystyle \frac{1}{2}}{\delta }_q^{p}R-\mathrm{\Lambda }{\delta }_q^p=(U^p{U}_q-{\delta }_q^p)C_2\dot{\mathrm{\Theta }}-{\displaystyle \frac{1}{2}}{\delta }_q^p{C}_2{\mathrm{\Theta }}^2+$$

(66)

$$+{\displaystyle \frac{\kappa }{c^2}}{\dot{A}}_3^2\left(-{\displaystyle \frac{1}{2}}{\delta }_q^p+{\delta }_3^p{\delta }_q^3+{\delta }_0^p{\delta }_q^0\right)+{\displaystyle \frac{1}{2}}{\delta }_q^p\kappa {\mathrm{\Psi }}_0^{2}V\left(\mathrm{\Phi }\right)+$$

$$+\kappa {\mathrm{\Psi }}_0^2\left[\left(U^p{U}_q-{\displaystyle \frac{1}{2}}{\delta }_q^p\right){\dot{\varphi }}^{\left(a\right)}{\dot{\varphi }}_{\left(a\right)}-{\displaystyle \frac{g^2}{c^2}}{A}_3^2\left({\delta }_3^p{\delta }_q^3-{\displaystyle \frac{1}{2}}{\delta }_q^p\right)f_{\left(c\right)\left(d\right)}^{\left(a\right)}{f}_{\left(a\right)\left(h\right)\left(e\right)}{\varphi }^{\left(d\right)}{\varphi }^{\left(e\right)}{q}^{\left(c\right)}{q}^{\left(h\right)}\right]$$

can be detailed for the Bianchi-I platform as follows:

$${\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{\dot{b}}{b}}+{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{\dot{c}}{c}}+{\displaystyle \frac{\dot{b}}{b}}{\displaystyle \frac{\dot{c}}{c}}-\mathrm{\Lambda }=-{\displaystyle \frac{1}{2}}{C}_2{\mathrm{\Theta }}^2+{\displaystyle \frac{\kappa }{2c^2}}{\dot{A}}_3^2+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^{2}V\left(\mathrm{\Phi }\right)+$$

(67)

$$+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^2{\dot{\varphi }}^{\left(a\right)}{\dot{\varphi }}_{\left(a\right)}+{\displaystyle \frac{g^2\kappa }{2c^2}}{\mathrm{\Psi }}_0^2{A}_3^2{f}_{\left(c\right)\left(d\right)}^{\left(a\right)}{f}_{\left(a\right)\left(h\right)\left(e\right)}{\varphi }^{\left(d\right)}{\varphi }^{\left(e\right)}{q}^{\left(c\right)}{q}^{\left(h\right)},$$

$${\displaystyle \frac{\ddot{b}}{b}}+{\displaystyle \frac{\ddot{c}}{c}}+{\displaystyle \frac{\dot{b}}{b}}{\displaystyle \frac{\dot{c}}{c}}-\mathrm{\Lambda }=-C_2\dot{\mathrm{\Theta }}-{\displaystyle \frac{1}{2}}{C}_2{\mathrm{\Theta }}^2-{\displaystyle \frac{\kappa }{2c^2}}{\dot{A}}_3^2+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^{2}V\left(\mathrm{\Phi }\right)-$$

(68)

$$-{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^2{\dot{\varphi }}^{\left(a\right)}{\dot{\varphi }}_{\left(a\right)}+{\displaystyle \frac{g^2\kappa }{2c^2}}{\mathrm{\Psi }}_0^2{A}_3^2{f}_{\left(c\right)\left(d\right)}^{\left(a\right)}{f}_{\left(a\right)\left(h\right)\left(e\right)}{\varphi }^{\left(d\right)}{\varphi }^{\left(e\right)}{q}^{\left(c\right)}{q}^{\left(h\right)},$$

$${\displaystyle \frac{\ddot{a}}{a}}+{\displaystyle \frac{\ddot{c}}{c}}+{\displaystyle \frac{\dot{a}}{a}}{\displaystyle \frac{\dot{c}}{c}}-\mathrm{\Lambda }=-C_2\dot{\mathrm{\Theta }}-{\displaystyle \frac{1}{2}}{C}_2{\mathrm{\Theta }}^2-{\displaystyle \frac{\kappa }{2c^2}}{\dot{A}}_3^2+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^{2}V\left(\mathrm{\Phi }\right)+$$

(69)

$$-{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^2{\dot{\varphi }}^{\left(a\right)}{\dot{\varphi }}_{\left(a\right)}+{\displaystyle \frac{g^2\kappa }{2c^2}}{\mathrm{\Psi }}_0^2{A}_3^2{f}_{\left(c\right)\left(d\right)}^{\left(a\right)}{f}_{\left(a\right)\left(h\right)\left(e\right)}{\varphi }^{\left(d\right)}{\varphi }^{\left(e\right)}{q}^{\left(c\right)}{q}^{\left(h\right)},$$

$${\displaystyle \frac{\ddot{b}}{b}}+{\displaystyle \frac{\ddot{a}}{a}}+{\displaystyle \frac{\dot{b}}{b}}{\displaystyle \frac{\dot{a}}{a}}-\mathrm{\Lambda }=-C_2\dot{\mathrm{\Theta }}-{\displaystyle \frac{1}{2}}{C}_2{\mathrm{\Theta }}^2+{\displaystyle \frac{\kappa }{2c^2}}{\dot{A}}_3^2+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^{2}V\left(\mathrm{\Phi }\right)+$$

(70)

$$-{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^2{\dot{\varphi }}^{\left(a\right)}{\dot{\varphi }}_{\left(a\right)}+{\displaystyle \frac{g^2\kappa }{2c^2}}{\mathrm{\Psi }}_0^2{A}_3^2{f}_{\left(c\right)\left(d\right)}^{\left(a\right)}{f}_{\left(a\right)\left(h\right)\left(e\right)}{\varphi }^{\left(d\right)}{\varphi }^{\left(e\right)}{q}^{\left(c\right)}{q}^{\left(h\right)}.$$

The sources in the right-hand sides of equations (68) and (69) coincide, and this fact allows us to assume that the spacetime possesses the local symmetry associated with the rotation around the $O{x}^3$ axis. This means that we can put $a\left(t\right)=b\left(t\right)$ and can omit, e.g., the equation (69) . In this situation one can work with two unknown functions $\mathrm{\Theta }\left(t\right)$ and $H\left(t\right)={\displaystyle \frac{\dot{a}}{a}}$ and make the replacement ${\displaystyle \frac{\dot{c}}{c}}\to \mathrm{\Theta }-2H$. In this terms the difference of equations (68) and (70) gives us the convenient consequence

$$\left(\dot{\mathrm{\Theta }}-3\dot{H}\right)+\mathrm{\Theta }\left(\mathrm{\Theta }-3H\right)=-{\displaystyle \frac{\kappa }{c^2\left(t\right)}}{\dot{A}}_3^2.$$

(71)

Similarly, the difference between the equations (70) and (67) gives the consequence

$$\dot{H}-H\left(\mathrm{\Theta }-3H\right)=-{\displaystyle \frac{1}{2}}{C}_2\dot{\mathrm{\Theta }}-{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^2{\nu }^2\left(t\right){\varphi }^{\left(a\right)}{\varphi }_{\left(a\right)}.$$

(72)

Since the Bianchi identities fulfill automatically on the solutions to the master equations for the vector and pseudoscalar fields, one of the equations (67)-(70) is a differential consequence of other ones. In fact, we can form the key set of master equations for the gravity field, if we take one of the equations (71) or (72) and add the equation

$$-3H^2+2H\mathrm{\Theta }-\mathrm{\Lambda }=-{\displaystyle \frac{1}{2}}{C}_2{\mathrm{\Theta }}^2+{\displaystyle \frac{\kappa }{2c^2}}{\dot{A}}_3^2+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^{2}V\left(\mathrm{\Phi }\right)+$$

(73)

$$+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^2\left[{\nu }^2\left(t\right)+{\displaystyle \frac{g^2}{c^2}}{A}_3^2\right]{\mathrm{\Phi }}^2\left(t_0\right)exp\left\{2{\int }_{t_0}^{t}d\tau \nu \left(\tau \right)\right\}.$$

4. Modeling of the Late-Time Evolution of the Multi-Component Pseudoscalar System

4.1. Key System of Equations

In order to obtain physically interesting results we formulate the so-called key system of equations with the following convenient notations and simplifications. First, we consider the quadratic potential of the pseudoscalar fields $V\left(\mathrm{\Phi }\right)=m_A^2{\mathrm{\Phi }}^2$, for which the term ${\displaystyle \frac{1}{2\mathrm{\Phi }}}{V}^{\prime }\left(\mathrm{\Phi }\right)$ converts into constant $m_A^2$, and thus, the equation for $\nu \left(t\right)$ is no longer integro-differential one:

$$\dot{\nu }+{\nu }^2+\mathrm{\Theta }\nu +g^2{A}^2+m_A^2=0.$$

(74)

Second, the new function defined as $A\left(t\right)={\displaystyle \frac{A_3}{c\left(t\right)}}$ satisfies the equation

$$\ddot{A}+\dot{A}\mathrm{\Theta }+A\left[\dot{\mathrm{\Theta }}-2\dot{H}+2H\mathrm{\Theta }-4H^2+g^2{\mathrm{\Psi }}_0^2{\mathrm{\Phi }}^2\left(t_0\right)exp\left\{2{\int }_{t_0}^{t}d\tau \nu \left(\tau \right)\right\}\right]=0.$$

(75)

Third, the functions $H\left(t\right)$ and $\mathrm{\Theta }\left(t\right)$ satisfy the pair of equations

$$\left(\dot{\mathrm{\Theta }}-3\dot{H}\right)+\mathrm{\Theta }\left(\mathrm{\Theta }-3H\right)=-\kappa {\left[\dot{A}+A(\mathrm{\Theta }-2H)\right]}^2,$$

(76)

$$-3H^2+2H\mathrm{\Theta }-\mathrm{\Lambda }=-{\displaystyle \frac{1}{2}}{C}_2{\mathrm{\Theta }}^2+{\displaystyle \frac{\kappa }{2}}{\left[\dot{A}+A(\mathrm{\Theta }-2H)\right]}^2+$$

(77)

$$+{\displaystyle \frac{1}{2}}\kappa {\mathrm{\Psi }}_0^2\left[m_A^2+{\nu }^2+g^2{A}^2\right]{\mathrm{\Phi }}^2\left(t_0\right)exp\left\{2{\int }_{t_0}^{t}d\tau \nu \left(\tau \right)\right\}.$$

Four self-consistent coupled integro-differential equations (74) - (77) for four required functions $\nu$, A, H, $\mathrm{\Theta }$ form the key system of equations.

4.2. Asymptotic Regime

We assume, first, that the cosmological constant is not vanishing, $\mathrm{\Lambda }\ne 0$; second, that at $t\to \infty$ the potential of the Yang-Mills and pseudoscalar fields tend to zero; third, the quantities H and $\mathrm{\Theta }$ tend to non-vanishing constants at $t\to \infty$. Clearly, we obtain from (76) that asymptotically

$$\mathrm{\Theta }-3H\to {\displaystyle \frac{\mathrm{const}}{a^{2}c}}\to 0,$$

(78)

and we deal with asymptotic isotropization of the Universe, ${\displaystyle \frac{\dot{c}}{c}}\to {\displaystyle \frac{\dot{a}}{a}}$. Taking into account (77) we obtain the asymptotic value of the Hubble function

$$H\to H_{\infty }=\sqrt{{\displaystyle \frac{\mathrm{\Lambda }}{3\left(1+\frac{3}{2}{C}_2\right)}}},$$

(79)

which includes the Jacobson constant $C_2$ and thus bears the imprint of the dynamic aether. The scale factors have the de Sitter form ${\displaystyle \frac{a\left(t\right)}{a\left(t_0\right)}}={\displaystyle \frac{c\left(t\right)}{c\left(t_0\right)}}=e^{H_{\infty }(t-t_0)}$, where $t_0$ is a certain reference point in time that marks the beginning of the asymptotic regime. In this case we can solve directly the transformed equation (74)

$${\displaystyle \frac{\dot{\nu }}{{\nu }^2+3H_{\infty }\nu +m_A^2}}=-1,$$

(80)

but the results depend, formally speaking, on the ratio between the parameters $m_A$ and $H_{\infty }$. We do not discuss here, what is the value of the ratio ${\displaystyle \frac{m_A}{H_{\infty }}}$ in reality, and list the results for all three standard situations.

4.2.1. Axion Mass Exceeds the Critical Value of the Hubble Constant $H_C={\displaystyle \frac{3}{2}}{H}_{\infty }$

When $m_A>{\displaystyle \frac{3}{2}}{H}_{\infty }$, the solution for $\nu \left(t\right)$ is

$$\nu \left(t\right)=-\left\{{\displaystyle \frac{3}{2}}{H}_{\infty }+\alpha tan\left[\alpha (t-t_{\ast })\right]\right\},$$

(81)

where we used the following notations

$$\alpha =\sqrt{m_A^2-{\displaystyle \frac{9}{4}}{H}_{\infty }^2},t_{\ast }=t_0+{\displaystyle \frac{1}{\alpha }}arctan\left[{\displaystyle \frac{\nu \left(t_0\right)+\frac{3}{2}{H}_{\infty }}{\alpha }}\right].$$

(82)

Respectively, the functions, which describe the multiplet of pseudoscalar fields, take the form

$${\varphi }^{\left(a\right)}\left(t\right)={\varphi }^{\left(a\right)}\left(t_0\right)e^{-{\displaystyle \frac{3}{2}}{H}_{\infty }(t-t_0)}\left[{\displaystyle \frac{cos\alpha (t-t_{\ast })}{cos\alpha (t_0-t_{\ast })}}\right].$$

(83)

Clearly, we deal with the regime of damped oscillations: the decrement of damping is equal to the Hubble critical parameter $H_C={\displaystyle \frac{3}{2}}{H}_{\infty }$, and the frequency of oscillations coincides with $\alpha$.

Formally speaking, the function $\nu \left(t\right)$ grows infinitely, when $t-t_{\ast }\to {\displaystyle \frac{\pi }{2}}$, however, the functions ${\varphi }^{\left(a\right)}\left(t\right)$ and their derivatives, containing $\nu \left(t\right)$

$${\dot{\varphi }}^{\left(a\right)}=-\left[{\displaystyle \frac{3}{2}}{H}_{\infty }+\alpha tan\left[\alpha (t-t_{\ast })\right]\right]{\varphi }^{\left(a\right)}\left(t_0\right)e^{-{\displaystyle \frac{3}{2}}{H}_{\infty }(t-t_0)}\left[{\displaystyle \frac{cos\alpha (t-t_{\ast })}{cos\alpha (t_0-t_{\ast })}}\right],$$

(84)

remain finite. When we calculate the stress-energy tensor of the pseudoscalar field

$$T_q^p={\mathrm{\Psi }}_0^2\left[{\displaystyle \frac{1}{2}}{\delta }_q^p{m}_A^2{\varphi }^{\left(a\right)}{\varphi }_{\left(a\right)}+\left(U^p{U}_q-{\displaystyle \frac{1}{2}}{\delta }_q^p\right){\dot{\varphi }}^{\left(a\right)}{\dot{\varphi }}_{\left(a\right)}\right],$$

(85)

we see that the corresponding scalars of energy density and pressure

$${\mathcal{W}}_A\left(t\right)={\displaystyle \frac{1}{2}}{\mathrm{\Psi }}_0^2{\mathrm{\Phi }}^2(m_A^2+{\nu }^2),\mathcal{P}\left(t\right)={\displaystyle \frac{1}{2}}{\mathrm{\Psi }}_0^2{\mathrm{\Phi }}^2(m_A^2-{\nu }^2)$$

(86)

also remain finite, when $t-t_{\ast }\to {\displaystyle \frac{\pi }{2}}$. Nevertheless, the growth of the function $\nu \left(t\right)$ can change the sign of pressure scalar. Indeed, if $m_A^2>{\nu }^2\left(t_0\right)$ and thus $\mathcal{P}\left(t_0\right)>0$, the moment of time $t=t_T$ exists

$$t_T=t_{\ast }+{\displaystyle \frac{1}{\sqrt{m_A^2-\frac{9}{4}{H}_{\infty }^2}}}arctan\sqrt{{\displaystyle \frac{m_A-\frac{3}{2}{H}_{\infty }}{m_A+\frac{3}{2}{H}_{\infty }}}},$$

(87)

when the pressure becomes equal to zero, and then remains negative at $t>t_T$.

In order to estimate the number density of pseudoscalar particles of the sort $\left(d\right)$ we can calculate the ratio of the energy density $W^{\left(d\right)}$ to the rest energy, i.e., to find the following quantities

$$N^{\left(d\right)}\left(t\right)={\displaystyle \frac{W^{\left(d\right)}\left(t\right)}{m_A}}=$$

(88)

$$={\displaystyle \frac{1}{2m_A}}{\mathrm{\Psi }}_0^2{\left[{\varphi }^{\left(d\right)}\left(t_0\right)\right]}^2{\displaystyle \frac{a^3\left(t_0\right)}{a^3\left(t\right)}}\left[m_A^2+{\left({\displaystyle \frac{3}{2}}{H}_{\infty }+\alpha tan\left[\alpha (t-t_{\ast })\right]\right)}^2\right]{\left[{\displaystyle \frac{cos\alpha (t-t_{\ast })}{cos\alpha (t_0-t_{\ast })}}\right]}^2.$$

These quantities decrease as $a^{-3}\left(t\right)$, are described by the oscillating functions and are predetermined by the initial values ${\varphi }^{\left(d\right)}\left(t_0\right)$.

4.2.2. Critical Value of the Hubble Constant Exceeds the Axion Mass

When $m_A<{\displaystyle \frac{3}{2}}{H}_{\infty }$ the solution for $\nu \left(t\right)$ takes the form

$$\nu \left(t\right)=-{\displaystyle \frac{3}{2}}{H}_{\infty }+\beta coth\beta (t-t_{\ast \ast }),$$

(89)

where $\beta =\sqrt{{\displaystyle \frac{9}{4}}{H}_{\infty }^2-m_A^2}$ and

$$t_{\ast \ast }=t_0+{\displaystyle \frac{1}{2\beta }}log\left[{\displaystyle \frac{\nu \left(t_0\right)+\frac{3}{2}{H}_{\infty }-\beta }{\nu \left(t_0\right)+\frac{3}{2}{H}_{\infty }+\beta }}\right].$$

(90)

The pseudoscalar fields can be now described by the regular functions

$${\varphi }^{\left(a\right)}\left(t\right)={\varphi }^{\left(a\right)}\left(t_0\right)e^{-{\displaystyle \frac{3}{2}}{H}_{\infty }(t-t_0)}\left[{\displaystyle \frac{sinh\beta (t-t_{\ast \ast })}{sinh\beta (t_0-t_{\ast \ast })}}\right],$$

(91)

which start to grow, reach maxima at $t_{\max }=t_{\ast \ast }+\mathrm{Arth}\sqrt{1-{\displaystyle \frac{4m_A^2}{9H_{\infty }^2}}}$, then decrease monotonically, and have the asymptotes

$${\varphi }^{\left(a\right)}(t\to \infty )\propto exp\left\{-{\displaystyle \frac{3}{2}}{H}_{\infty }t\left[1-\sqrt{1-{\displaystyle \frac{4m_A^2}{9H_{\infty }^2}}}\right]\right\}.$$

(92)

The asymptotic behavior of the function $N^{\left(d\right)}\left(t\right)$ is now described by the law

$$N^{\left(d\right)}\left(t\right)\propto {\left[a\left(t\right)\right]}^{-3\mu },\mu =\left[1-\sqrt{1-{\displaystyle \frac{4m_A^2}{9H_{\infty }^2}}}\right].$$

(93)

4.2.3. Critical Value of the Hubble Constant Coincides with the Axion Mass

When $m_A={\displaystyle \frac{3}{2}}{H}_{\infty }$, the results of integration can be written as

$$\nu \left(t\right)=-{\displaystyle \frac{3}{2}}{H}_{\infty }+{\displaystyle \frac{\left[\nu \left(t_0\right)+\frac{3}{2}{H}_{\infty }\right]}{1+\left[\nu \left(t_0\right)+\frac{3}{2}{H}_{\infty }\right](t-t_0)}},$$

(94)

$${\varphi }^{\left(a\right)}\left(t\right)={\varphi }^{\left(a\right)}\left(t_0\right)e^{-{\displaystyle \frac{3}{2}}{H}_{\infty }(t-t_0)}\left\{1+\left[\nu \left(t_0\right)+{\displaystyle \frac{3}{2}}{H}_{\infty }\right](t-t_0)\right\}.$$

(95)

The functions ${\varphi }^{\left(a\right)}\left(t\right)$ start to grow, reach maxima at $t_{\max }=t_0+{\displaystyle \frac{2\nu \left(t_0\right)}{3H_{\infty }\left(\nu \left(t_0\right)+\frac{3}{2}{H}_{\infty }\right)}}$ and then decrease monotonically.

The function $N^{\left(d\right)}\left(t\right)$ asymptotically behaves as

$$N^{\left(d\right)}\left(t\right)\propto a^{-3}{log}^{2}a.$$

(96)

4.3. Asymptotic Behavior of the Gauge Field

In the asymptotic regime the equation (75) transforms into

$$\ddot{A}+3H_{\infty }\dot{A}+2H_{\infty }^{2}A=0,$$

(97)

and its general solution containing two constants of integration $B_1$ and $B_2$ is

$$A\left(t\right)=B_1{e}^{-H_{\infty }t}+B_2{e}^{-2H_{\infty }t}.$$

(98)

Keeping in mind that $A_3=c\left(t\right)A$, and $c\left(t\right)=c\left(t_0\right)e^{H_{\infty }(t-t_0)}$, we obtain, that the physically motivated solution should contain only one term

$$A_3\left(t\right)=A_3\left(t_0\right)e^{-H_{\infty }(t-t_0)}.$$

(99)

The energy density of the gauge field behaves as follows:

$$W_{\left(G\right)}={\displaystyle \frac{{\dot{A}}_3^2}{2c^2\left(t\right)}}={\displaystyle \frac{A_3^2\left(t_0\right)H_{\infty }^2}{2c^2\left(t_0\right)}}{e}^{-4H_{\infty }(t-t_0)}.$$

(100)

The energy density is proportional to $a^{-4}\left(t\right)$; this behavior is typical for the electromagnetic field in the expanding Universe.

5. Conclusions and Outlook

In the introduction of this paper we formulated the ultimate goal of our work as the need to show that the evolution and decay of the color aether, that occurred in the early Universe, had an important "by-product" - the emergence of a SU(N) symmetric multiplet of pseudoscalar fields, which evolved in the late-time Universe into the multi-component dark matter. In other words, we intended to show that the color aether, which was presented by the SU(N) symmetric multiplet of vector fields $U^{\left(a\right)j}$, turned into the dynamic aether described by the unit vector field singlet $U^j$, and the pseudoscalar particles, which were previously designed in the form of an axionic multiplet, became a filler for several fractions of the cosmic dark matter. The first part of this task has been successfully solved: the complete mathematical model describing this scenario, is established and some physical conclusions are made. In particular, it was shown that the Yang-Mills field, that originally arose only to maintain the SU(N) symmetry of the color aether, formed the sources in the right-hand sides of the equations for the pseudoscalar field multiplet (25), so that a non-zero pseudoscalar field was forced to appear. Let us emphasize that the standard Peccei-Quinn source-term ${\displaystyle \frac{1}{4{\mathrm{\Psi }}_0^2}}{F}_{mn}^{\left(a\right)\ast }{F}_{\left(a\right)}^{mn}$ contains the complete convolutions of the Yang-Mills field strength and its dual tensor only, and thus it is colorless and is capable of generating the axion singlet only. However, the color aether "suggested" a supplementary tool for modeling, the vector in the group space, ${\omega }^{\left(a\right)}$, constructed using the covariant divergence of the vector fields (see (9) and (10)). Acquiring this new tool gives a possibility to activate principally new, the color sensitive sources in (25), which produce the whole axionic multiplet ${\varphi }^{\left(a\right)}$. If such a scenario works, the emergence of multi-axionic configuration can be realized, and next problem is connected with the evolution of this field configuration. We have described above one of the possible lines of further development of this story in the framework of the Bianchi-I cosmological model, based on the hypothesis of color polarization of the SU(N) symmetric system; although it is obvious that this story may have other outcomes.

We have left aside for now the second important question: how to explain the difference in masses of pseudoscalar particles that form different fractions of the cosmic dark matter? In the presented model all the particles of the pseudoscalar multiplet have the same mass $m_A$, and the information about its value is encoded in the potential $V\left(\mathrm{\Phi }\right)$. In the next work we plan to consider the potential $V\left(\mathrm{\Psi }\right)$, where ${\mathrm{\Psi }}^2={\varphi }^{\left(a\right)}{\varphi }^{\left(b\right)}\left[G_{\left(a\right)\left(b\right)}+\mu {\omega }_{\left(a\right)}{\omega }_{\left(b\right)}\right]$, and the parameter $\mu$ describes the difference of the particle masses. Clearly, the potential $V\left(\mathrm{\Psi }\right)$, containing the quantity ${\omega }^{\left(a\right)}$, depends now not only on the pseudoscalar fields ${\varphi }^{\left(a\right)}$, but also on the vector fields $U^{\left(b\right)j}$, on the Yang-Mills potentials $A_k^{\left(d\right)}$ and on the spacetime metric. This means that new source-terms will appear in the Jacobson, Yang-Mills, and Einstein equations significantly complicating their analysis. This task is beyond the scope of the presented work, but will be solved in the near future.