Interaction of the Cosmic Dark Fluid with Dynamic Aether: Parametric Mechanism of Axion Generation in the Early Universe

Abstract

We consider an isotropic homogeneous cosmological model with five interacting elements: first, the dynamic aether presented by a unit timelike vector field; second, the pseudoscalar field describing an axionic component of the dark matter; third, the cosmic dark energy, described by a rheologic fluid; fourth, the non-axionic dark matter coupled to the dark energy; and fifth, the gravity field. We show that the early evolution of the Universe described by this model can include two specific epochs: the first one can be characterized as a super-inflation epoch; the second epoch is associated with an oscillatory regime. The dynamic aether carries out a regulatory mission; the rheologic dark fluid provides the specific features of the spacetime evolution. The oscillations of the scale factor and of the Hubble function are shown to switch on the parametric (Floquet-type) mechanism of the axion number growth.

1. Introduction

The paper is dedicated to the memory of Steven Weinberg (AB: half a century ago, the excellent book [1] predetermined my scientific life).

1.1. Motivation of the Work

1.1.1. Inflation vs. Super-Inflation

We live and work in a unique situation, when we obtain new sensational results of observations made on the James Webb Space Telescope (JWST) practically in real time (for information see, e.g., the official website webb.nasa.gov (accessed on 23 August 2023)). These new data concern, in particular, the discovery of Mothra, an extremely magnified monster star; estimations of the masses of warm dark matter particles and of axion dark matter particles [2]; the abundance of carbon-containing molecules [3], etc. In addition, starting from the discovery of gravitational waves in 2015 [4], when two black holes with masses 36 $M_{\mathrm{Sun}}$ and 29 $M_{\mathrm{Sun}}$ collided, more than a hundred events of this type have been recorded. These events are associated with the merger and collision of black holes, whose masses are in the range from 20 to 90 solar masses.

Many of the data obtained by JWST and the LIGO–VIRGO collaboration are so unexpected that they make us think about the revision of the models of the early evolution of our Universe. What is, from our point of view, the main element of such a revision? We think that the starting period of the Universe evolution, when the size of the Universe does not exceed the size of the region allowing macroscopic causal processes, has to be longer than we think (of course, for comparison, we have to use the time scale that is predetermined by the corresponding value of the effective Hubble function H associated with that epoch).

Such a possibility appears, for instance, if we consider super-inflation. The term super-inflation has been already used, e.g., in the model of Loop Quantum Cosmology [5], in order to mark a specific episode of the Universe inflationary evolution, when the kinetic energy of the scalar field is much more than the potential energy. Our goal is to consider not approximate but exact solutions describing the super-inflation as an alternative to standard inflation. Of course, the inflation scenario has already explained many details of the early Universe evolution, and when we pose a question about a super-inflation, we have to motivate this step. Keeping in mind this simple argument, we would like to attract attention to one new fact only. The LIGO–VIRGO collaboration has proved that black holes with intermediate masses (from 20 to 90 solar masses) do exist. Astrophysicists are ready to theoretically explain the presence of black holes with masses of several solar masses obtained in the scenario of a star collapse; also, one can explain the existence of super-massive black holes. However, now there is no adequate theory for the formation of the medium-sized black holes and super-massive stars. But if the causal period of the early Universe evolution lasted longer than the inflation theory predicts, we have a natural opportunity to explain the observed set of the black hole masses. We hope that a corresponding model for the formation of medium-sized black holes will be formulated in the near future: for example, similar to how the problem of the causal limit of the neutron star maximum mass was solved in [6].

Mathematically, one can explain the mentioned idea by comparing two functions $y_1=e^{H(t-t_0)}$ and $y_2=exp[\alpha sinhH(t-t_0)]$. Both functions start with the same values $y_1\left(t_0\right)=1=y_2\left(t_0\right)$. The first function describes the standard inflation, and we can require that at the moment $t=t_{*}$, the scale factor increased ${10}^{26}$ times, i.e., $\frac{a\left(t_{*}\right)}{a\left(t_0\right)}=e^{H(t_{*}-t_0)}={10}^{26}$. We obtain in this case the well-known 60 e-folds as follows: $H(t_{*}-t_0)=26ln10\approx 60$. If we use the second function, which describes the so-called super-inflation, and again require that $\frac{a\left(t_{*}\right)}{a\left(t_0\right)}=exp[\alpha sinhH(t_{*}-t_0)]={10}^{26}$, we have to assume that $\alpha \approx 120·{10}^{-26}$. Thus, for small values of time $H(t-t_0)\to 0$, we can decompose the mentioned functions as $y_1\approx 1+H(t-t_0)$ and $y_2\approx 1+120·{10}^{-26}H(t-t_0)$. Finally, we compare the time moments $t_1$ and $t_2$ for which the sizes of the expanding Universe become of the same order as the causal domain size $L_{\mathrm{causal}}$. We obtain that $t_1-t_0=120·{10}^{-26}(t_2-t_0)$, or equivalently, $t_2-t_0\propto {10}^{24}(t_1-t_0)$. This means that the super-exponential growth, on the one hand, guarantees that at $t=t_{*}$ the Universe expanded ${10}^{26}$ times as due to the standard inflation; on the other hand, the causal period of the Universe evolution lasts much longer, ensuring the development of causal phenomena, which could be associated, for example, with the formation of proto-galaxies, proto-stars, medium-sized black holes, etc.

Super-inflation can be described, in particular, as the exact solution of the set of master equations in the framework of models with a dark fluid of rheological type. In the works [7,8], we obtained the super-inflationary solution assuming that the equation of state of dark energy contains the convective derivative of the pressure; we also assumed that the dynamics of the dark matter particles is under control of the Archimedean-type force induced by dark energy. In the work [9], the super-inflationary solution appears as the exact solution of the dark fluid model with the kernel of contact interaction of the integral Volterra type. In the work [10], the mentioned solution appeared when we used the integral representation of the equations of state of the dark energy and dark matter. In other words, the solutions of the super-inflationary type seem to be typical ones for the rheological models of the cosmic dark fluid.

In this work, we consider again the rheological models for the dark fluid; however, now, we assume that the dark matter is a multi-component substratum [11]. We separate the axionic component of the dark matter [12] and consider axions on the language of field theory as a pseudoscalar field with modified periodic potential. Other components of the dark matter (WIMPs, ALPs, warm and hot dark matter parts, etc.) are unified and described as a dark medium with rheological properties.

1.1.2. Dynamic Aether as a Guiding Element of the Cosmic Evolution

The concept of the dynamic aether [13,14,15,16] gives the theorist a unique tool for modeling the process of controlling cosmic expansion. The dynamic aether is described by the unit-time-like global vector field, which is associated with the aether velocity four-vector $U^j$ and realizes the idea of a privileged frame of reference [17,18].

In addition to the geometric aspects of control, the dynamic aether can control the rhythm of life in the Universe. The fact is that the scalar of expansion defined as the divergence of the velocity four-vector $\Theta ={\nabla }_k{U}^k$ coincides with the tripled Hubble function $\Theta =3H$, when the Universe is isotropic and homogeneous. Thus, for the spacetimes of the FLRW type, the scalar $\mathcal{T}=\frac{3}{\Theta }$ predetermines the typical time scale of the Universe evolution. Traditional history of the Universe is written using energy units and the equivalent temperature: the main stages of the Universe evolution are tied to some milestones associated with breaking of symmetry of fundamental interactions. However, for many purposes, we need to link the energy scale (or temperature scale) with the appropriate time scale. Clearly, requirements of the covariant approach do not give us possibility to use time-dependent parameters of equations of state or the time-dependent cosmological constant directly. Of course, one needs to introduce appropriate scalars, associated with some field, to solve the dynamic equation for this field and then to reconstruct the required guiding scalar. In this sense, the unit vector field $U^j$ is well suited for this role, since the evolution of this field is well described by the Jacobson equations [13,14,15], and the basic scalar $\Theta ={\nabla }_j{U}^j$ can be associated with the cosmological time scale.

Keeping in mind this idea we introduced in [19], the mechanism of the aetheric control on the axion field evolution is through introduction of the guiding function ${\Phi }_{*}(\Theta )$ into the axion field potential $V(\varphi ,{\Phi }_{*}(\Theta ))$. Such a modification of the periodic axion potential led to the emergence of a concept of equilibrium states in the axion containing systems, for which the axion field $\varphi$ takes the values $\varphi =n{\Phi }_{*}$ with an integer n [20]. In this context, the value of the expansion scalar $\Theta$ predetermines the position and depth of minima of the axion field potential. The extension of this approach has shown that for cosmological and astrophysical applications, it would be interesting to enlarge the number of guiding functions, which could be constructed using the covariant derivative of the aether velocity four-vector. For instance, when we deal with cosmological models of the Bianchi types, we cannot ignore the fact that in addition to the expansion scalar, we can use the non-vanishing symmetric traceless shear tensor describing the aether flow, ${\sigma }_{mn}$, in the theory. Correspondingly, the new geometric aspects of the Einstein–Maxwell–aether theory can be associated with the term ${\sigma }^{mn}{F}_{mp}{F}_n^p$ in the modified Lagrangian; the new geometric aspect of the Einstein–aether–axion theory can be connected with additional part of the axion kinetic energy ${\sigma }^{mn}{\nabla }_m\varphi {\nabla }_n\varphi$. As for the description of the aetheric control over the axion system, one can add the square of the shear tensor ${\sigma }^2={\sigma }_{mn}{\sigma }^{mn}$ as an argument of the guiding function ${\Phi }_{*}(\Theta ,{\sigma }^2)$. For the static spherically symmetric model, the square of the acceleration four-vector $a^2$ may be in demand; the square of the antisymmetric vorticity tensor ${\omega }^2={\omega }_{mn}{\omega }^{mn}$ could appear as an argument of the guiding function in the model of Gödel type, describing the rotating Universe.

1.1.3. Interaction of the Dynamic Aether with Dark Fluid

All the mentioned extensions of the Einstein–aether theory are formulated in the language of field theory. The description of the dark fluid is made in terms of relativistic phenomenological hydrodynamics of two-component fluids. In order to realize the idea of aetheric control over the dark fluid, we suggest to include the scalars $\Theta$, ${\sigma }^2$, $a^2$, and ${\omega }^2$ into the Lagrangian $L_{\left(\mathrm{DF}\right)}$ of the dark fluid. Such an approach can be indicated as the semi-phenomenological one. In this work, we restrict ourselves by the ansatz that the function $L_{\left(\mathrm{DF}\right)}$ has a multi-step structure. This representation is based on the Heaviside step functions, arguments of which contain the expansion scalar $\Theta$. Using this approach, we assume that the aether divides the history of the Universe evolution into episodes, which can be indicated as inflation, super-inflation, oscillatory stage, etc. In this division into episodes, some critical values of the expansion scalar appear, ${\Theta }_{*}^{\left(1\right)}$, ${\Theta }_{*}^{\left(2\right)}$, ${\Theta }_{*}^{\left(3\right)}$, etc., which are assumed to play roles analogous to the roles of critical temperatures $T_{*}^{\left(1\right)}$, $T_{*}^{\left(2\right)}$, and $T_{*}^{\left(3\right)}$ in the series of phase transitions associated with the Universe restructuring.

1.1.4. The Role of Axionic Dark Matter in Our Approach

We support the point of view that the axionic component is the key constructive element of the multi-component dark matter. The history of investigations (theoretical and experimental) of the axionic dark matter phenomenon (see, e.g., [21,22,23,24,25,26,27,28,29]) hints to us that the anomalous growth of the axion number in the early Universe could take place (the so-called “axionization” of the Universe); now, these relic axions form the basic part of cold dark matter. Following this idea, we consider two mechanisms of instability in the axion system provoked by the dark fluid controlled by the dynamic aether. The first mechanism can be realized in the scheme of super-inflationary expansion of the Universe. The second mechanism can be switched on at the oscillatory stage of the Universe evolution; it can be associated with the parametric instability described by the Floquet theorem for the Hill equation with periodic coefficients.

1.2. The Structure of the Work

The paper is organized as follows. In Section 2, we describe the mathematical formalism and derive the master equations of the presented theory. In Section 3, we consider applications of this theory to the model of evolution of the isotropic homogeneous Universe filled with the dynamic aether, two-component dark fluid, and axion field. In Section 3.1, we reduce the basic master equations to the chosen spacetime symmetry and find the exact solution to the equations for the unit vector field (aether velocity). In Section 3.2, we focus on the super-inflationary scenario of the Universe evolution; we find exact solutions for the scale factor and Hubble function, and we reconstruct the state functions for the rheologically active dark energy and non-axionic dark matter; in Section 3.2.3, we consider the problem of instability in the axion system controlled by the dynamic aether. In Section 3.3, we study the oscillatory episode of the Universe evolution; we again find the exact solution for the geometric quantities and for the dark fluid state functions; in Section 3.3.3, we focus on the analysis of the parametric mechanism of the axion generation. Section 4 includes discussion and conclusions.

2. The Formalism

2.1. The Total Action Functional

We consider the extension of the Einstein–aether–axion model and add the term associated with the dark fluid $L_{\left(\mathrm{DF}\right)}$ to the total Lagrangian. The total action functional

$$-S_{\left(\mathrm{EA}\right)}=\int d^{4}x\sqrt{-g}\left\{\frac{1}{2\kappa }\left[R+2\Lambda +\lambda \left(g_{mn}{U}^m{U}^n-1\right)+K_{mn}^{ab}{\nabla }_a{U}^m{\nabla }_b{U}^n\right]+\frac{1}{2}{\Psi }_0^2\left(V-{\nabla }_k\varphi {\nabla }^k\varphi \right)+L_{\left(\mathrm{DF}\right)}\right\}$$

(1)

contains the standard Einstein–Hilbert term with the determinant of the metric g, the covariant derivative ${\nabla }_k$, the Ricci scalar R, the cosmological constant $\Lambda$, and the Einstein constant $\kappa =\frac{8\pi G}{c^4}$. The unit-time-like vector field $U^j$ is associated with the velocity four-vector of the dynamic aether (see, e.g., [13,14,15,16] for history, mathematical details, and definitions); the term with the Lagrange multiplier $\lambda$ in front is introduced to provide the normalization of the vector field, $U^k{U}_k=1$. The constitutive tensor

$$K_{mn}^{ab}=C_1{g}^{ab}{g}_{mn}+C_2{\delta }_m^a{\delta }_n^b+C_3{\delta }_n^a{\delta }_m^b+C_4{U}^a{U}^b{g}_{mn}$$

(2)

contains four phenomenological constants $C_1$, $C_2$, $C_3$, and $C_4$. The pseudoscalar field $\varphi$ stands to describe the axionic part of the dark matter; the term V describes the potential of the axion field; the parameter ${\Psi }_0$ relates to the coupling constant of the axion–photon interaction $g_{A\gamma \gamma }$, $\frac{1}{{\Psi }_0}=g_{A\gamma \gamma }$. The potential of the axion field

$$V(\varphi ,{\Phi }_{*})=\frac{m_A^2{\Phi }_{*}^2}{2{\pi }^2}\left[1-cos\left(\frac{2\pi \varphi }{{\Phi }_{*}}\right)\right]$$

(3)

is periodic since it inherits the discrete symmetry $\frac{2\pi \varphi }{{\Phi }_{*}}\to \frac{2\pi \varphi }{{\Phi }_{*}}+2\pi n$ (n is an integer). Also, this potential can be indicated as modified periodic potential since it contains the guiding function ${\Phi }_{*}$, which describes an averaged value of the pseudoscalar field; generally, this guiding function depends on coordinates via the scalars associated with the model as a whole (see, e.g., [19]). This periodic potential has its minima at $\varphi =n{\Phi }_{*}$; for these values of the axion field, the potential and its first derivative vanish, $V_{|\varphi =n{\Phi }_{*}}=0$, ${\frac{\partial V}{\partial \varphi }}_{|\varphi =n{\Phi }_{*}}=0$. We indicate the states $\varphi =n{\Phi }_{*}$ as the equilibrium state. Near the minima, when $\varphi \to n{\Phi }_{*}+\psi$ and $\left|\frac{2\pi \psi }{{\Phi }_{*}}\right|$ is small, the potential takes the standard form $V\to m_A^2{\psi }^2$, where $m_A$ is the axion rest mass. The integer n describes the level on which the axion field is fixed; when $n=1$, we deal with the basic level; when $n=0$, we assume that the axions are absent.

If the guiding function ${\Phi }_{*}$ is constant, the more convenient term can be used for this quantity, namely, the vacuum expectation value. In this regard, it is important to mention that, unlike the axion theory, the standard theory of the dynamic aether does not contain an appropriate vector field potential. In this sense, the vacuum expectation value of the vector field does not appear in the standard version of this theory; thus, there is no fixed direction in the space, and the spatial isotropy violation can not take place.

The term $L_{\left(\mathrm{DF}\right)}$ describes the dark fluid, which indicates two interacting cosmic substrates: the dark energy and the non-axionic dark matter.

2.2. Auxiliary Elements of Analysis

Based on the velocity four-vector of the aether $U^j$, we decompose all the tensor quantities highlighting the longitudinal and transversal components. In particular, the covariant derivative can be decomposed as follows

$${\nabla }_k=U_{k}D+{\stackrel{\perp }{\nabla }}_k,D=U^s{\nabla }_s,{\stackrel{\perp }{\nabla }}_k={\Delta }_k^j{\nabla }_j,{\Delta }_k^j={\delta }_k^j-U^j{U}_k.$$

(4)

${\Delta }_k^j$ is the projector. The covariant derivative ${\nabla }_k{U}_j$ can be decomposed in the standard sum

$${\nabla }_k{U}_j=U_{k}D{U}_j+{\sigma }_{kj}+{\omega }_{kj}+\frac{1}{3}{\Delta }_{kj}\Theta ,$$

(5)

where the acceleration four-vector $D{U}_j$, the symmetric traceless shear tensor ${\sigma }_{kj}$, the skew-symmetric vorticity tensor ${\omega }_{kj}$ and the expansion scalar $\Theta$ are presented by the well-known formulas

$$D{U}_j=U^s{\nabla }_s{U}_j,{\sigma }_{kj}=\frac{1}{2}\left({\stackrel{\perp }{\nabla }}_k{U}_j+{\stackrel{\perp }{\nabla }}_j{U}_k\right)-\frac{1}{3}{\Delta }_{kj}\Theta ,{\omega }_{kj}=\frac{1}{2}\left({\stackrel{\perp }{\nabla }}_k{U}_j-{\stackrel{\perp }{\nabla }}_j{U}_k\right),\Theta ={\nabla }_k{U}^k.$$

(6)

Using the decomposition (5) we form one linear and three quadratic scalars

$$\Theta ={\nabla }_k{U}^k,a^2=D{U}_{k}D{U}^k,{\sigma }^2={\sigma }_{mn}{\sigma }^{mn},{\omega }^2={\omega }_{mn}{\omega }^{mn}.$$

(7)

Generally, the guiding function ${\Phi }_{*}$ is assumed to be the function of all four scalars, ${\Phi }_{*}(\Theta ,a^2,{\sigma }^2,{\omega }^2)$, and we prepare the basic formulas for the general case. But below, we consider the application of the theory to the spatially isotropic homogeneous model, for which $D{U}^j=0$, ${\sigma }_{mn}=0$, ${\omega }_{mn}=0$; naturally, we will focus on the case when the guiding function depends on the expansion scalar $\Theta$ only, i.e., ${\Phi }_{*}={\Phi }_{*}(\Theta )$. The scalar $\mathcal{K}\equiv K^{abmn}\left({\nabla }_a{U}_m\right)\left({\nabla }_b{U}_n\right)$ can also be expressed via the mentioned invariants as follows:

$$\mathcal{K}=(C_1+C_4)D{U}_{k}D{U}^k+(C_1+C_3){\sigma }_{ik}{\sigma }^{ik}+(C_1-C_3){\omega }_{ik}{\omega }^{ik}+\frac{1}{3}\left(C_1+3C_2+C_3\right){\Theta }^2.$$

(8)

Clearly, for the spatially isotropic and homogeneous model, this term reduces to $\mathcal{K}=\frac{1}{3}\left(C_1+3C_2+C_3\right){\Theta }^2$.

2.3. Master Equations of the Model

2.3.1. Master Equations for the Unit Vector Field

Variations of the action functional (1) with respect to the Lagrange multiplier $\lambda$ and to the four-vector $U^i$ give the following set of equations, respectively:

$$g_{mn}{U}^m{U}^n=1,$$

(9)

$${\nabla }_a{\mathcal{J}}_j^a=\lambda U_j+C_{4}D{U}_m{\nabla }_j{U}^m+\frac{1}{2}\kappa {\Psi }_0^2\frac{\delta V}{\delta U^j}+\kappa \frac{\delta L_{\left(\mathrm{DF}\right)}}{\delta U^j}.$$

(10)

The tensor ${\mathcal{J}}_j^a$ is of the form

$${\mathcal{J}}_j^a=K_{jn}^{ab}{\nabla }_b{U}^n=C_1{\nabla }^a{U}_j+C_2{\delta }_j^a\Theta +C_3{\nabla }_j{U}^a+C_4{U}^{a}D{U}_j.$$

(11)

The variational derivatives in (10) can be represented as follows:

$$\begin{array}{c}\hfill \frac{1}{2}\kappa {\Psi }_0^2\frac{\delta V}{\delta U^j}\Rightarrow -{\nabla }_j\left(\Omega \frac{\partial {\Phi }_{*}}{\partial \Theta }\right)-2D{U}_{j}D\left(\Omega \frac{\partial {\Phi }_{*}}{\partial a^2}\right)-{\nabla }^n\left(2\Omega \frac{\partial {\Phi }_{*}}{\partial {\sigma }^2}{\sigma }_{jn}\right)+{\nabla }^n\left(2\Omega \frac{\partial {\Phi }_{*}}{\partial {\omega }^2}{\omega }_{jn}\right)+\\ \hfill +2\Omega \frac{\partial {\Phi }_{*}}{\partial a^2}\left[D{U}_k{\nabla }_j{U}^k-\Theta D{U}_j-D^2{U}_j\right]-2\Omega \frac{\partial {\Phi }_{*}}{\partial {\sigma }^2}D{U}^n{\sigma }_{jn}-2\Omega \frac{\partial {\Phi }_{*}}{\partial {\omega }^2}D{U}^n{\omega }_{jn},\end{array}$$

(12)

$$\begin{array}{c}\hfill \kappa \frac{\delta L_{\left(\mathrm{DF}\right)}}{\delta U^j}\Rightarrow \kappa \left\{-{\nabla }_j\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right)-2D{U}_{j}D\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial a^2}\right)-{\nabla }^n\left(2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial {\sigma }^2}{\sigma }_{jn}\right)+{\nabla }^n\left(2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial {\omega }^2}{\omega }_{jn}\right)+\right.\\ \hfill \left.+2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial a^2}\left[D{U}_k{\nabla }_j{U}^k-\Theta D{U}_j-D^2{U}_j\right]-2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial {\sigma }^2}D{U}^n{\sigma }_{jn}-2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial {\omega }^2}D{U}^n{\omega }_{jn}\right\}.\end{array}$$

(13)

Here, we used the convenient definition

$$\Omega =\frac{1}{2}\kappa {\Psi }_0^2\frac{\partial V}{\partial {\Phi }_{*}}=\frac{\kappa {\Psi }_0^2{m}_A^2}{2{\pi }^2}\left\{{\Phi }_{*}\left[1-cos\left(\frac{2\pi \varphi }{{\Phi }_{*}}\right)\right]-\pi \varphi sin\left(\frac{2\pi \varphi }{{\Phi }_{*}}\right)\right\}.$$

(14)

The Lagrange multiplier can be formally presented as the sum $\lambda ={\lambda }_{\left(\mathrm{U}\right)}+{\lambda }_{\left(\mathrm{V}\right)}+{\lambda }_{\left(\mathrm{DF}\right)}$, where

$${\lambda }_{\left(\mathrm{U}\right)}=U^j{\nabla }_a{\mathcal{J}}_j^a-C_{4}D{U}_{m}D{U}^m,{\lambda }_{\left(\mathrm{V}\right)}=-\frac{1}{2}\kappa {\Psi }_0^2{U}^j\frac{\delta V}{\delta U^j},{\lambda }_{\left(\mathrm{DF}\right)}=-\kappa U^j\frac{\delta L_{\left(\mathrm{DF}\right)}}{\delta U^j}.$$

(15)

Using (12) and (13), we obtain immediately

$${\lambda }_{\left(\mathrm{V}\right)}=D\left(\Omega \frac{\partial {\Phi }_{*}}{\partial \Theta }\right)-2\Omega {\sigma }^2\frac{\partial {\Phi }_{*}}{\partial {\sigma }^2}-2\Omega {\omega }^2\frac{\partial {\Phi }_{*}}{\partial {\omega }^2}-4\Omega a^2\frac{\partial {\Phi }_{*}}{\partial a^2},$$

(16)

$${\lambda }_{\left(\mathrm{DF}\right)}=\kappa \left\{D\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right)-2{\sigma }^2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial {\sigma }^2}-2{\omega }^2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial {\omega }^2}-4a^2\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial a^2}\right\}.$$

(17)

When we deal with the spatially isotropic homogeneous model, only the first terms in the right-hand sides of the Equations (16) and (17) are non-vanishing. In the state of the axionic equilibrium, when $\varphi =n{\Phi }_{*}$, we obtain that $\Omega =0$.

2.3.2. Master Equation for the Axion Field

Variation of the total action functional with respect to the axion field yields

$${\nabla }^k{\nabla }_k\varphi +\frac{m_A^2{\Phi }_{*}}{2\pi }sin\left(\frac{2\pi \varphi }{{\Phi }_{*}}\right)=0.$$

(18)

If we deal with the basic equilibrium state of the axion field and put $\varphi ={\Phi }_{*}$ into (18), we can conclude that the guiding function ${\Phi }_{*}$ satisfies the massless Klein–Gordon equation

$${\nabla }^k{\nabla }_k{\Phi }_{*}=0.$$

(19)

2.3.3. Master Equations for the Gravitational Field

Variation with respect to the metric gives the gravity field equations:

$$R_{ik}-\frac{1}{2}R{g}_{ik}-\Lambda g_{ik}=T_{ik}^{\left(\mathrm{U}\right)}+\kappa T_{ik}^{\left(\mathrm{A}\right)}+\kappa T_{ik}^{\left(\mathrm{DF}\right)}.$$

(20)

The first term in the right-hand side of (20) is the standard stress–energy tensor associated with the aether flow [13,14,15]:

$$\begin{gathered} T_{ik}^{\left(\mathrm{U}\right)}=\frac{1}{2}{g}_{ik}{K}^{abmn}{\nabla }_a{U}_m{\nabla }_b{U}_n+{\lambda }_{\left(\mathrm{U}\right)}{U}_i{U}_k+ \\ +{\nabla }^m\left[U_{(i}{\mathcal{J}}_{k)m}-{\mathcal{J}}_{m(i}{U}_{k)}-{\mathcal{J}}_{\left(ik\right)}{U}_m\right]+ \\ +C_1\left[\left({\nabla }_m{U}_i\right)\left({\nabla }^m{U}_k\right)-\left({\nabla }_i{U}_m\right)\left({\nabla }_k{U}^m\right)\right]+C_{4}D{U}_{i}D{U}_k. \end{gathered}$$

(21)

The parentheses symbolize the symmetrization of indices. The term ${\lambda }_{\left(\mathrm{U}\right)}{U}_i{U}_k$ is included into $T_{ik}^{\left(\mathrm{U}\right)}$ since ${\lambda }_{\left(\mathrm{U}\right)}$ (see (15)) contains the aether velocity and its derivatives only. The second term

$$\begin{gathered} \kappa T_{ik}^{\left(\mathrm{A}\right)}=\kappa {\Psi }_0^2\left[{\nabla }_i\varphi {\nabla }_k\varphi +\frac{1}{2}{g}_{ik}\left(V-{\nabla }_s\varphi {\nabla }^s\varphi \right)\right]+{\lambda }_{\left(\mathrm{V}\right)}{U}_i{U}_k-g_{ik}(D+\Theta )\left(\Omega \frac{\partial {\Phi }_{*}}{\partial \Theta }\right)+ \\ +2\Omega \frac{\partial {\Phi }_{*}}{\partial a^2}D{U}_{i}D{U}_k+2{\nabla }_s\left\{\Omega \frac{\partial {\Phi }_{*}}{\partial a^2}\left[D{U}^s{U}_i{U}_k-2U^s{U}_{(i}D{U}_{k)}\right]\right\}- \\ -2{\nabla }_s\left[\Omega \frac{\partial {\Phi }_{*}}{\partial {\sigma }^2}{U}^s{\sigma }_{ik}\right]-4\Omega \frac{\partial {\Phi }_{*}}{\partial {\sigma }^2}\left[D{U}_n{\sigma }_{(i}^n{U}_{k)}+{\sigma }_{(i}^n{\omega }_{k)n}\right]+ \\ +4{\nabla }_s\left[\Omega \frac{\partial {\Phi }_{*}}{\partial {\omega }^2}{U}_{(i}{\omega }_{k)·}^s\right]-4\Omega \frac{\partial {\Phi }_{*}}{\partial {\omega }^2}\left[D{U}^n{U}_{(i}{\omega }_{k)n}+{\sigma }_{(i}^n{\omega }_{k)n}\right] \end{gathered}$$

(22)

describes the contribution of the axion field coupled to the aether. The term ${\lambda }_{\left(\mathrm{V}\right)}{U}_i{U}_k$ is included into this part of the total stress–energy tensor, since it describes the interaction of the aether with the axion field. The stress–energy tensor of the dark fluid $T_{ik}^{\left(\mathrm{DF}\right)}$ is presented by two terms:

$$T_{ik}^{\left(\mathrm{DF}\right)}=\frac{1}{\kappa }{\lambda }_{\left(\mathrm{DF}\right)}{U}_i{U}_k+\frac{(-2)}{\sqrt{-g}}\frac{\delta }{\delta g^{ik}}\left[\sqrt{-g}{L}_{\left(\mathrm{DF}\right)}\right].$$

(23)

The first term $\frac{1}{\kappa }{\lambda }_{\left(\mathrm{DF}\right)}{U}_i{U}_k$ (see (15)) is the last element of the total term $\lambda U_i{U}_k$, and the second one is given by the standard variation formula. This last construction can be decomposed algebraically using the aether velocity four-vector $U^j$,s and we obtain

$$T_{ik}^{\left(\mathrm{DF}\right)}=\left[\frac{1}{\kappa }{\lambda }_{\left(\mathrm{DF}\right)}+\mathcal{W}\right]U_i{U}_k+q_i{U}_k+q_k{U}_i+{\mathcal{P}}_{ik}.$$

(24)

Here, $\mathcal{W}$ denotes the total energy density of the dark fluid; $q_j$ is the heat-flux four-vector, and ${\mathcal{P}}_{ik}$ is the pressure tensor; these quantities are defined in the standard way:

$$\frac{1}{\kappa }{\lambda }_{\left(\mathrm{DF}\right)}+\mathcal{W}=U^i{T}_{ik}^{\left(\mathrm{DF}\right)}{U}^k,q_j={\Delta }_j^i{T}_{ik}^{\left(\mathrm{DF}\right)}{U}^k,{\mathcal{P}}_{ik}={\Delta }_i^p{T}_{pq}^{\left(\mathrm{DF}\right)}{\Delta }_k^q.$$

(25)

When we deal with the isotropic homogeneous cosmological model, we have to put $q_j=0$ and ${\mathcal{P}}_{ik}=-\mathcal{P}{\Delta }_{ik}$. We assume that $\mathcal{W}=W+E$ is the sum of the energy densities of the dark energy and of the non-axionic dark matter; similarly, $\mathcal{P}=P+\Pi$ is the sum of the corresponding scalars of pressure.

2.3.4. Bianchi Identity and Conservation Law

The Bianchi identity requires that

$${\nabla }^k\left[T_{ik}^{\left(\mathrm{U}\right)}+\kappa T_{ik}^{\left(\mathrm{A}\right)}+\kappa T_{ik}^{\left(\mathrm{DF}\right)}\right]=0.$$

(26)

We would like to emphasize two important details. First, since we included the term ${\lambda }_{\left(\mathrm{U}\right)}{U}_i{U}_k$ into the stress–energy tensor $T_{ik}^{\left(\mathrm{U}\right)}$, we have guaranteed that ${\nabla }^k{T}_{ik}^{\left(\mathrm{U}\right)}=0$ on the solutions of the master Equation (10). Second, the presence of the term ${\lambda }_{\left(\mathrm{V}\right)}{U}_i{U}_k$ in the stress–energy tensor $T_{ik}^{\left(\mathrm{A}\right)}$ guarantees that ${\nabla }^k{T}_{ik}^{\left(\mathrm{A}\right)}=0$ on the solutions of the master Equation (18). Thus, the equality ${\nabla }^k{T}_{ik}^{\left(\mathrm{DF}\right)}=0$ is the consequence of the Bianchi identity on the one hand, and it gives the evolutionary equation for the state functions of the dark fluid on the other hand.

3. Application to the Spatially Isotropic Homogeneous Cosmological Model

3.1. The Spacetime Platform and Reduced Master Equations

We work below with the spacetime of the FLRW type with the metric

$$d{s}^2=d{t}^2-a^2\left(t\right)\left(d{x}^2+d{y}^2+d{z}^2\right).$$

(27)

The symmetry of the model hints that one can search for the velocity four-vector of the aether in the form $U^j={\delta }_0^j$. It is well known that for this metric, the covariant derivative of the vector field

$${\nabla }_k{U}_i=\frac{1}{2}{\dot{g}}_{ik}$$

(28)

is symmetric, and we see explicitly that $D{U}_j=0$, ${\sigma }_{mn}=0$, ${\omega }_{mn}=0$. The expansion scalar is non-vanishing

$$\Theta =3\frac{\dot{a}}{a}=3H\left(t\right)\Rightarrow {\nabla }_k{U}_i={\Delta }_{ik}H\left(t\right).$$

(29)

Here, $H\left(t\right)=\frac{\dot{a}}{a}$ is the Hubble function, and the dot denotes, as usual, the derivative with respect to the cosmological time t (here and below $c=1$).

3.1.1. Solution to the Equations of the Vector Field

When we start to analyze the reduced Equations (10) and (11), we have to emphasize the following circumstance. The sum of the coupling constants $C_1+C_3$ can be estimated as $-6\times {10}^{-15}<C_1+C_3<1.4\times {10}^{-15}$, and below, we put $C_1+C_3=0$. This estimation is based on the result of observation of the binary neutron star merger (the events GW170817 and GRB 170817A [30]), which has shown that the ratio of the velocities of the gravitational and electromagnetic waves satisfies the inequalities $1-3\times {10}^{-15}<\frac{v_{\mathrm{gw}}}{c}<1+7\times {10}^{-16}$), while according to [15], the square of the velocity of the tensorial aether mode is equal to $S_{\left(2\right)}^2=\frac{1}{1-(C_1+C_3)}$. Keeping in mind that $D{U}_j=0$, ${\sigma }_{mn}=0$, and ${\omega }_{mn}=0$, we find that (11) converts into

$$J_j^a=C_2\Theta {\delta }_j^a,$$

(30)

and the equations for the unit vector field (10) now takes the form

$$C_2{\nabla }_j\Theta =\lambda U_j-{\nabla }_j\left(\Omega \frac{\partial {\Phi }_{*}}{\partial \Theta }\right)-\kappa {\nabla }_j\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right).$$

(31)

Clearly, the set of Equation (31) contains only one non-trivial equation, which gives, in fact, the solution for the Lagrange multiplier $\lambda$:

$$D\left(C_2\Theta +\Omega \frac{\partial {\Phi }_{*}}{\partial \Theta }+\kappa \frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right)=\lambda .$$

(32)

3.1.2. Reduced Equation for the Axion Field

For the fixed spacetime symmetry, the Equation (18) takes the form

$$\ddot{\varphi }+3H\dot{\varphi }+\frac{m_A^2{\Phi }_{*}}{2\pi }sin\left(\frac{2\pi \varphi }{{\Phi }_{*}}\right)=0.$$

(33)

Let us imagine that the axion field $\varphi$ is frozen in the first minimum of the axion potential; it coincides with the guiding function, and thus, it follows the variations of ${\Phi }_{*}\left(t\right)$. Then, we put $\varphi ={\Phi }_{*}$ into (33), and we have to admit that the axion field remains in the first minimum if the guiding function ${\Phi }_{*}$ satisfies the equation

$${\ddot{\Phi }}_{*}+3H{\dot{\Phi }}_{*}=0.$$

(34)

Our ansatz is that (34) defines the missing master equation for the guiding function ${\Phi }_{*}$ in the framework of the established model.

Clearly, the Equation (34) admits the first integral

$${\dot{\Phi }}_{*}\left(t\right)=\frac{\mathrm{const}}{a^3\left(t\right)}={\dot{\Phi }}_{*}\left(t_0\right){\left[\frac{a\left(t_0\right)}{a\left(t\right)}\right]}^3.$$

(35)

Here and below, the parameter $t_0$ describes the time moment, starting from which our semi-phenomenological model is valid. We assume that $t_0>t_{\mathrm{quantum}}$, so that the results of evolution on the stage associated with the quantum description of the Universe are fixed in the initial data, e.g., in ${\dot{\Phi }}_{*}\left(t_0\right)$.

3.1.3. Key Equations for the Gravity Field

Solving the Einstein equations for the isotropic homogeneous cosmological model, we follow the standard scheme: we consider the gravity field Equation (20) for the indices $i=0$, $k=0$, and the consequence of the Bianchi identity (26) for $i=0$, thus obtaining the first and second key equations of the model. Equation (18) for the axion field forms the third key equation. The first key equation can be written as

$$3H^2\left(1+\frac{3}{2}{C}_2\right)-\Lambda =\kappa {\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right)}^{·}+\frac{1}{2}\kappa {\Psi }_0^2\left(V+{\dot{\varphi }}^2\right)-\Omega \Theta \frac{\partial {\Phi }_{*}}{\partial \Theta }+\kappa (W+E).$$

(36)

The second key equation is presented by the balance equation for the dark fluid

$$\dot{W}+\dot{E}+3H(W+E+P+\Pi )=-\kappa \left[{\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right)}^{··}+3H{\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right)}^{·}\right].$$

(37)

3.2. Super-Inflationary Scenario of Early Universe Evolution

The novelty of our approach is associated with a multi-step-like structure of the function $L_{\left(\mathrm{DF}\right)}$. Let us explain this idea for the one-step function. We assume that

$$L_{\left(\mathrm{DF}\right)}=\eta ({\Theta }_{*}-\Theta )L_{\left(\mathrm{DF}\right)}^{\left(1\right)}+\eta (\Theta -{\Theta }_{*})L_{\left(\mathrm{DF}\right)}^{\left(2\right)}.$$

(38)

Here, $\eta \left(F\right)$ is the Heaviside function, which is equal to one if the argument F is positive, $F>0$, and is equal to zero when $F<0$. The constant ${\Theta }_{*}\equiv \Theta \left(t_{*}\right)$ is the value of the expansion scalar at some fixed time moment $t_{*}>t_0$. We assume that both functions $L_{\left(\mathrm{DF}\right)}^{\left(1\right)}$ and $L_{\left(\mathrm{DF}\right)}^{\left(2\right)}$ depend on time, but do not include $\Theta$. Generally, $L_{\left(\mathrm{DF}\right)}^{\left(1\right)}\left(t\right)\ne L_{\left(\mathrm{DF}\right)}^{\left(2\right)}\left(t\right)$; however, the continuity condition is met, $L_{\left(\mathrm{DF}\right)}^{\left(1\right)}\left(t_{*}\right)=L_{\left(\mathrm{DF}\right)}^{\left(2\right)}\left(t_{*}\right)$. This means that we have to solve the key equations of the model for two intervals: $t_0<t<t_{*}$ and $t>t_{*}$ and sew the found state functions at $t=t_{*}$. In fact, we guarantee that the derivative ${\left(\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }\right)}^{·}$ vanishes in both mentioned intervals, as well as in the point $t=t_{*}$. To be more precise, we obtain

$$\frac{\partial L_{\left(\mathrm{DF}\right)}}{\partial \Theta }=\delta (\Theta -{\Theta }_{*})\left[L_{\left(\mathrm{DF}\right)}^{\left(2\right)}-L_{\left(\mathrm{DF}\right)}^{\left(1\right)}\right]=0,$$

(39)

if the function $\Theta \left(t\right)$ is monotonic and the equation $\Theta \left(t_{*}\right)={\Theta }_{*}$ has only one solution $t_{*}$.

3.2.1. Epoch of the Dark Fluid Domination

Let us assume that during the first episode of the Universe evolution, when $t_0<t<t_{*}$, the axion field is absent, $\varphi =0$; thus, $\Omega =0$ and $V=0$, so the key equations of the model as a whole are reduced to the following two equations:

$$H^2=\frac{{\Lambda }_{*}}{3}+\frac{{\kappa }_{*}}{3}(W+E),$$

(40)

$$\dot{W}+\dot{E}+3H(W+E+P+\Pi )=0.$$

(41)

Here, we used the definitions ${\Lambda }_{*}=\frac{\Lambda }{1+\frac{3}{2}{C}_2}$ and ${\kappa }_{*}=\frac{\kappa }{1+\frac{3}{2}{C}_2}$ and assumed that $C_2>-\frac{2}{3}$. The next convenient step is to introduce the new variable $x=\frac{a\left(t\right)}{a\left(t_0\right)}$, providing the following relationships:

$$\dot{W}=xH\frac{dW}{dx},t-t_0={\int }_1^{\frac{a\left(t\right)}{a\left(t_0\right)}}\frac{dx}{xH\left(x\right)}.$$

(42)

In these terms, Equation (41) can be rewritten as

$$x{W}^{\prime }+3(W+P)=-\left[x{E}^{\prime }+3(E+\Pi )\right].$$

(43)

Following the idea of modeling of the so-called kernel of interaction Q (see, e.g., [31,32]), we can rewrite (43) as a pair of equations

$$x{W}^{\prime }+3(W+P)=Q\left(x\right),x{E}^{\prime }+3(E+\Pi )=-Q\left(x\right).$$

(44)

The next step is to propose the constitutive equations for dark energy and for non-axionic dark matter and to formulate the structure of the kernel $Q\left(x\right)$. We assume that during the first episode of early Universe evolution, the following relationships are appropriate:

$$P=(\Gamma -1)W,\Pi =(\gamma -1)E,Q\left(x\right)=K_0{\int }_1^x\frac{dy}{y}[E\left(y\right)-W\left(y\right)].$$

(45)

Here, $\Gamma$, $\gamma$, and $K_0$ are some phenomenological constants. The representation of the kernel Q in the integral form was motivated in the work [9]. Using (45), we obtain the integral equation of the Volterra type

$$x{W}^{\prime }+3\Gamma W=K_0{\int }_1^x\frac{dy}{y}(E\left(y\right)-W\left(y\right)),$$

(46)

from which one can extract the energy density of the non-axionic dark matter

$$E\left(x\right)=W\left(x\right)+\frac{1}{K_0}\left[x^2{W}^{″}+x{W}^{\prime }(1+3\Gamma )\right].$$

(47)

Then, we put the function $E\left(x\right)$ into (43) and obtain the Euler equation of the third order for the function $W\left(x\right)$:

$$x^3{W}^{\prime \prime \prime }+3x^2{W}^{″}(1+\Gamma +\gamma )+x{W}^{\prime }\left[2K_0+(1+3\Gamma )(1+3\gamma )\right]+3K_0(\Gamma +\gamma )W=0.$$

(48)

The characteristic equation for this differential Euler equation is of the form

$${\sigma }^3+3(\Gamma +\gamma ){\sigma }^2+\sigma (2K_0+9\Gamma \gamma )+3K_0(\Gamma +\gamma )=0.$$

(49)

We are especially interested in the analysis of the case with the threefold degenerate root ${\sigma }_1={\sigma }_2={\sigma }_3=0$ of this characteristic equation. Clearly, such a situation takes place when

$$\Gamma +\gamma =0,K_0=\frac{9}{2}{\gamma }^2.$$

(50)

For this special case, the solutions for $W\left(x\right)$ and $E\left(x\right)$ are of the form

$$W\left(x\right)=W\left(1\right)\left(1-3\Gamma logx\right)+\frac{9}{4}{\Gamma }^2\left[W\left(1\right)+E\left(1\right)\right]{log}^{2}x,$$

(51)

$$E\left(x\right)=E\left(1\right)\left(1+3\Gamma logx\right)+\frac{9}{4}{\Gamma }^2\left[W\left(1\right)+E\left(1\right)\right]{log}^{2}x.$$

(52)

Here, we took into account the consequence of the integral relationship (46), which gives $W^{\prime }\left(1\right)=-3\Gamma W\left(1\right)$ (similarly, we see that $E^{\prime }\left(1\right)=3\Gamma E\left(1\right)$). For the square of the Hubble function, we obtain from (40)

$$\frac{3}{{\kappa }_{*}}\left[H^2\left(x\right)-H^2\left(1\right)\right]=3\Gamma \left[E\left(1\right)-W\left(1\right)\right]logx+\frac{9}{2}{\gamma }^2\left[W\left(1\right)+E\left(1\right)\right]{log}^{2}x,$$

(53)

where the following definition was used:

$$H^2\left(1\right)=\frac{{\Lambda }_{*}}{3}+\frac{{\kappa }_{*}}{3}\left[W\left(1\right)+E\left(1\right)\right].$$

(54)

3.2.2. Exact Explicit Solutions in the Model of Super-Inflation

One can obtain now the analytic formula which links the scale factor and the cosmological time for arbitrary initial data $W\left(1\right)$ and $E\left(1\right)$. For illustration only, we put $E\left(1\right)=W\left(1\right)$ and obtain the integral

$$t-t_0=\pm {\int }_1^{\frac{a\left(t\right)}{a\left(t_0\right)}}\frac{dx}{x\sqrt{H^2\left(t_0\right)+3{\kappa }_{*}{\gamma }^{2}W\left(t_0\right){log}^{2}x}}.$$

(55)

We choose the sign plus in this formula, and direct integration gives the super-exponential law of evolution of the scale factor

$$log\left[\frac{a\left(t\right)}{a\left(t_0\right)}\right]=\frac{H\left(t_0\right)}{\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}}sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right].$$

(56)

When $W\left(t_0\right)$ tends to zero, i.e., dark fluid is absent, this formula gives the de Sitter-type law

$$\underset{W\left(t_0\right)\to 0}{lim}\left[\frac{a\left(t\right)}{a\left(t_0\right)}\right]=e^{\sqrt{\frac{{\Lambda }_{*}}{3}}(t-t_0)}$$

(57)

with the aetherically modified cosmological constant ${\Lambda }_{*}=\frac{\Lambda }{1+\frac{3}{2}{C}_2}$ (see Figure 1 for illustration).

In terms of the cosmological time, the Hubble function $H\left(t\right)$ is presented by the function

$$H\left(t\right)=H\left(t_0\right)cosh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right],$$

(58)

which grows monotonically. The acceleration parameter is of the form

$$-q\left(t\right)\equiv \frac{\ddot{a}}{a{H}^2}=1+\frac{\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}}{H\left(t_0\right)}\left\{\frac{sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]}{{cosh}^2\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]}\right\}.$$

(59)

The energy density scalars for the dark energy and for the non-axionic dark matter take the form, respectively,

$$W\left(t\right)=W\left(t_0\right)\left\{1+\frac{\sqrt{3}H\left(t_0\right)}{\sqrt{{\kappa }_{*}W\left(t_0\right)}}sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]+\frac{3H^2\left(t_0\right)}{2{\kappa }_{*}W\left(t_0\right)}{sinh}^2\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]\right\},$$

(60)

$$E\left(t\right)=W\left(t_0\right)\left\{1-\frac{\sqrt{3}H\left(t_0\right)}{\sqrt{{\kappa }_{*}W\left(t_0\right)}}sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]+\frac{3H^2\left(t_0\right)}{2{\kappa }_{*}W\left(t_0\right)}{sinh}^2\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]\right\}.$$

(61)

3.2.3. On the Stability of the Axion Field Configuration

We assumed above that the axion field is in the equilibrium state with $n=0$. Since the scale factor is found, we can calculate guiding function ${\Phi }_{*}\left(t\right)$ using (35), obtaining the formal integral

$${\Phi }_{*}\left(t\right)={\Phi }_{*}\left(t_0\right)+{\dot{\Phi }}_{*}\left(t_0\right){\int }_{t_0}^{t}d\tau exp\left\{-\frac{\sqrt{3}H\left(t_0\right)}{\gamma \sqrt{{\kappa }_{*}W\left(t_0\right)}}sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(\tau -t_0)\right]\right\}.$$

(62)

This integral cannot be presented via elementary functions, but can be expressed via incomplete modified Bessel functions [33]. Clearly, if we consider ${\dot{\Phi }}_{*}\left(t_0\right)=0$ and ${\Phi }_{*}\left(t_0\right)=2\pi$, we obtain the potential $V\left(\varphi \right)=2m_A^2\left[1-cos\varphi \right]$, which does not admit control over the axion evolution carried out by the dynamic aether.

Now, we admit that a fluctuation occurs $\varphi \to \psi \ne 0$ with $|\frac{2\pi \psi }{{\Phi }_{*}}|<<1$, so that the master equation for the pseudoscalar field converts into

$$\ddot{\psi }+3H\left(t\right)\dot{\psi }+m_A^2\psi =0,$$

(63)

where $H\left(t\right)$ is given by (58). Using the replacement $\psi \to {\left[\frac{a\left(t\right)}{a\left(t_0\right)}\right]}^{-\frac{3}{2}}\Psi$, we obtain the so-called modified Hill equation

$$\begin{gathered} \ddot{\Psi }=\mathcal{J}\left(t\right)\Psi ,\mathcal{J}\left(t\right)=-m_A^2+\frac{9}{4}{H}^2\left(t\right)+\frac{3}{2}\dot{H}= \\ =\left[\frac{9}{8}{H}^2\left(t_0\right)-m_A^2\right]+\frac{9}{8}{H}^2\left(t_0\right)cosh\left[2\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]+ \\ +\frac{3}{2}H\left(t_0\right)\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]. \end{gathered}$$

(64)

This terminology is based on the idea that if we replace the hyperbolic functions cosh and sinh by the trigonometric cos and sin, we obtain the standard Hill equation with periodic Hill potential. The Hill potential $\mathcal{J}\left(t\right)$ (64) is monotonic; thus, there exists a moment $t=t_{+}$, so that $\mathcal{J}(t>t_{+})>0$. In other words, when $t>t_{+}$, the Hill potential is positive, meaning that the term $\frac{\ddot{\Psi }}{\Psi }$ is also positive (see (64)). This means that if $\Psi \left(t_{+}\right)>0$, the function $\Psi$ grows exponentially, while the factor ${\left[\frac{a\left(t\right)}{a\left(t_0\right)}\right]}^{-\frac{3}{2}}$ means that the function $\psi =\Psi {\left[\frac{a\left(t\right)}{a\left(t_0\right)}\right]}^{-\frac{3}{2}}$ decreases super-exponentially. Clearly, the function $\psi \left(t\right)$ grows; then, it reaches the maximum at $t=t_{\max }$ and then decreases inevitably under the influence of the Universe expansion. One can use the following estimation of the behavior of this perturbation:

$$\psi \left(t\right)\propto cosh\left[2\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]exp\left\{-\frac{\sqrt{3}H\left(t_0\right)}{2\gamma \sqrt{{\kappa }_{*}W\left(t_0\right)}}sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t-t_0)\right]\right\}.$$

(65)

We deal with the instability in the axion field behavior on the interval $t_{+}<t<t_{\max }$. Indeed, when $W\left(t_0\right)\ne 0$, the function $\psi \left(t\right)$ starts to grow and reaches the maximal value, with the height of the maximum determined by the initial value $W\left(t_0\right)$ (see Figure 2 for illustration).

The special case $W\left(t_0\right)=0$ relates to the case of constant Hubble function $H>0$, for which Equation (63) converts into the differential equation with constant coefficients, and the set of solutions of the corresponding characteristic equation

$$k_{1,2}=-\frac{3}{2}H\pm \sqrt{\frac{9}{4}{H}^2-m^2}$$

(66)

do not contain real positive roots.

3.3. Periodic Episode of the Early Universe Evolution

3.3.1. Modifications of the Model

Now, we assume that on the quantum stage of the Universe evolution, which corresponds to the time interval $0<t<t_0$, standard inflation already led to the expansion characterized by 60 e-folds. We expect that the oscillatory regime of the Universe evolution will be switched on; however, in itself, such an event is unlikely. But, the aether, which controls the behavior of the dark fluid, is able to organize a fine-tuning to provide such an opportunity. Indeed, let us consider the following composite dark fluid model: first, during the interval $t_0<t<t_{*}$, the dark fluid is described in the same way as in the previous section; but, at the moment $t=t_{*}$, the aether organizes the restructuring of the equations of state for the dark energy and non-axionic dark matter. Now, we deal with $\Gamma =0$ for typical dark energy, and $\gamma =1$ for typical dust. Also, we assume that the kernel of the contact interaction in the dark fluid is of the local type, but the dark energy now possesses self-interaction of the rheological type (see [10] for details). Mathematically, this means that we use the replacements

$$Q\left(x\right)=K_0{\int }_{x_{*}}^x\frac{dy}{y}[E\left(y\right)-W\left(y\right)]\Rightarrow \tilde{Q}\left(x\right)={\omega }_0[E\left(x\right)-W\left(x\right)],$$

(67)

$$P\left(x\right)=(\Gamma -1)W\left(x\right)\Rightarrow \tilde{P}\left(x\right)=-W+K_1{\int }_{x_{*}}^x\frac{dy}{y}{\left(\frac{y}{x}\right)}^{\nu }W\left(y\right).$$

(68)

Here, $x_{*}=\frac{a\left(t_{*}\right)}{a\left(t_0\right)}$. Then, instead of (44) with (45), we consider the modified coupled balance equations

$$x{W}^{\prime }+3K_1{\int }_{x_{*}}^x\frac{dy}{y}{\left(\frac{y}{x}\right)}^{\nu }W\left(y\right)={\omega }_0(E-W),x{E}^{\prime }+3E={\omega }_0(W-E).$$

(69)

One can derive the key differential equations of the third order for $E\left(x\right)$ from the pair of the integro-differential Equations (69) if we extract $W\left(x\right)$ from the second Equation (69), put it into the first one, and fulfill the appropriate differentiation. The key equation has the form of the Euler equation

$$x^3{E}^{\prime \prime \prime }+x^2{E}^{″}\left(6+\nu +2{\omega }_0\right)+x{E}^{\prime }\left(3K_1+4+5{\omega }_0+4\nu +2\nu {\omega }_0\right)+3E\left[K_1\left(3+{\omega }_0\right)+\nu {\omega }_0\right]=0.$$

(70)

The corresponding characteristic equation

$${\sigma }^3+{\sigma }^2\left(3+\nu +2{\omega }_0\right)+\sigma \left(3K_1+3{\omega }_0+3\nu +2\nu {\omega }_0\right)+3\left[K_1\left(3+{\omega }_0\right)+\nu {\omega }_0\right]=0$$

(71)

has three coinciding roots ${\sigma }_1={\sigma }_2={\sigma }_3=0$, when

$$\nu =-(3+2{\omega }_0),K_1=\frac{{\omega }_0(3+2{\omega }_0)}{(3+{\omega }_0)},$$

(72)

and ${\omega }_0$ satisfies the following cubic equation:

$$4{\omega }_0^3+15{\omega }_0^2+27{\omega }_0+27=0.$$

(73)

This equation has only one real root ${\omega }_0\approx -2.06$. This unique set of parameters gives us the submodel with the following functions $E\left(x\right)$ and $W\left(x\right)$:

$$E\left(x\right)=E\left(x_{*}\right)+x_{*}{E}^{\prime }\left(x_{*}\right)log\left(\frac{x}{x_{*}}\right)+{\tilde{C}}_3{log}^2\left(\frac{x}{x_{*}}\right),$$

(74)

$$W\left(x\right)=W\left(x_{*}\right)+x_{*}{W}^{\prime }\left(x_{*}\right)log\left(\frac{x}{x_{*}}\right)+{\tilde{C}}_3\frac{(3+{\omega }_0)}{{\omega }_0}{log}^2\left(\frac{x}{x_{*}}\right).$$

(75)

The modified initial data at $t=t_{*}$ ($x=x_{*}$) can be obtained from (69); they now have the form:

$$x_{*}{W}^{\prime }\left(x_{*}\right)={\omega }_0[E\left(x_{*}\right)-W\left(x_{*}\right)],x_{*}{E}^{\prime }\left(x_{*}\right)={\omega }_{0}W\left(x_{*}\right)-(3+{\omega }_0)E\left(x_{*}\right),$$

(76)

so that the constant ${\tilde{C}}_3$ can be presented in two equivalent versions:

$$2{\tilde{C}}_3={\omega }_0{x}_{*}{W}^{\prime }\left(x_{*}\right)-(3+{\omega }_0)x_{*}{E}^{\prime }\left(x_{*}\right)=E\left(x_{*}\right)\left[{\omega }_0^2+{(3+{\omega }_0)}^2\right]-{\omega }_{0}W\left(x_{*}\right)(3+2{\omega }_0).$$

(77)

The square of the Hubble function can be presented as follows:

$$\begin{gathered} \frac{3}{{\kappa }_{*}}\left[H^2\left(x\right)-H^2\left(x_{*}\right)\right]= \\ -3E\left(x_{*}\right)log\left(\frac{x}{x_{*}}\right)+\left(\frac{3+2{\omega }_0}{2{\omega }_0}\right)\left[E\left(x_{*}\right)\left(2{\omega }_0^2+6{\omega }_0+9\right)-W\left(x_{*}\right){\omega }_0(3+2{\omega }_0)\right]{log}^2\left(\frac{x}{x_{*}}\right), \end{gathered}$$

(78)

where we used the following definition:

$$H^2\left(x_{*}\right)=\frac{{\Lambda }_{*}}{3}+\frac{{\kappa }_{*}}{3}\left[W\left(x_{*}\right)+E\left(x_{*}\right)\right].$$

(79)

The principal detail of our model is that the aether fixes the moment $t_{*}$ (and thus the value $x_{*}$) so that the coefficient in front of ${log}^2\left(\frac{x}{x_{*}}\right)$ in Formula (78) is negative. Keeping in mind that $E\left(x_{*}\right)$ and $W\left(x_{*}\right)$ are fixed by the Formulas (61) and (60) with $t=t_{*}$, respectively, we see that this is possible when

$$0.41<\frac{\sqrt{3}H\left(t_0\right)}{\sqrt{{\kappa }_{*}W\left(t_0\right)}}sinh\left[\gamma \sqrt{3{\kappa }_{*}W\left(t_0\right)}(t_{*}-t_0)\right]<4.91.$$

(80)

(Here, we used ${\omega }_0\approx -2.06$). Since the hyperbolic sinus is the monotonic function, we can find the appropriate value $t_{*}>t_0$ for every set of parameters $H\left(t_0\right)>0$, ${\kappa }_{*}>0$, $W\left(t_0\right)>0$, and $\gamma >0$.

3.3.2. Exact Solutions in the Periodic Model

In order to find the scale factor and the Hubble function for the mentioned case, we can use the simplified method. We introduce the auxiliary definitions

$${\alpha }^2=-{\kappa }_{*}\left(\frac{3+2{\omega }_0}{6{\omega }_0}\right)\left[E\left(x_{*}\right)\left(2{\omega }_0^2+6{\omega }_0+9\right)-W\left(x_{*}\right){\omega }_0(3+2{\omega }_0)\right],2\beta ={\kappa }_{*}E\left(x_{*}\right),$$

(81)

and search for the $a\left(t\right)$ and $H\left(t\right)$ in the form

$$a\left(t\right)=a_{X}exp\left\{h_{X}sin\left[\alpha (t-t_{*})+\phi \right]\right\},a_X=a\left(x_{*}\right)exp\left(-\frac{\beta }{{\alpha }^2}\right),$$

(82)

$$H\left(t\right)=\frac{\dot{a}}{a}=\alpha h_{X}cos\left[\alpha (t-t_{*})+\phi \right].$$

(83)

Then, we put these quantities into (78) and obtain the relationship between the introduced parameters

$${\alpha }^2{h}_X^2=H^2\left(t_{*}\right)+\frac{{\beta }^2}{{\alpha }^2}.$$

(84)

The last point is to find the auxiliary parameter $\phi$ which plays the role of a phase of the metric oscillations. We require that at $t=t_{*}$, the following initial condition takes place

$$a\left(t_{*}\right)=a_{X}exp\left[h_{X}sin\phi \right],H\left(t_{*}\right)=\alpha h_{X}cos\phi ,$$

(85)

and thus

$$sin\phi =\frac{\beta }{h_X{\alpha }^2}\le 1,cos\phi =\frac{H\left(t_{*}\right)}{\alpha h_X}\le 1\Rightarrow tan\phi =\frac{\beta }{H\left(t_{*}\right)\alpha }.$$

(86)

Illustration of the behavior of the scale factor is presented on Figure 3.

This exact solution demonstrates that the Hubble function is periodic with a period $T=\frac{2\pi }{\alpha }$; the parameter $\alpha$ plays the role of the frequency of the oscillations of the spacetime. The energy densities of the non-axionic dark matter and dark energy are also periodic:

$$E\left(t\right)=E_X+\frac{6{\omega }_0{\alpha }^2{h}_X}{{\kappa }_{*}{(3+2{\omega }_0)}^2}sin\left[\alpha (t-t_{*})+\phi \right]-\frac{3{\omega }_0{\alpha }^2{h}_X^2}{{\kappa }_{*}(3+2{\omega }_0)}{sin}^2\left[\alpha (t-t_{*})+\phi \right],$$

(87)

$$W\left(t\right)=W_X-\frac{6{\omega }_0{\alpha }^2{h}_X}{{\kappa }_{*}{(3+2{\omega }_0)}^2}sin\left[\alpha (t-t_{*})+\phi \right]-\frac{3{\alpha }^2{h}_X^2(3+{\omega }_0)}{{\kappa }_{*}(3+2{\omega }_0)}{sin}^2\left[\alpha (t-t_{*})+\phi \right].$$

(88)

Here, for the sake of convenience, we marked the values of the energies at the moment when the argument of sinus takes a zero value by $E_X$ and $W_X$. The energy densities of the dark energy (88) and of the non-axionic dark matter (87) describe oscillations not only with the frequency $\alpha$, but also with the double frequency $2\alpha$.

Finally, one can find the guiding function ${\Phi }_{*}\left(t\right)$ in the integral form

$${\Phi }_{*}\left(t\right)={\Phi }_{*}\left(t_{*}\right)+{\dot{\Phi }}_{*}\left(t_{*}\right)\frac{1}{\alpha }{e}^{\frac{3\beta }{2{\alpha }^2}}{\int }_{{\tau }_1}^{{\tau }_2}d\tau e^{-\frac{3}{2}{h}_{X}cos\tau },$$

(89)

where

$${\tau }_2={\tau }_1+\alpha (t-t_{*}),{\tau }_1=arctan\left[\frac{\beta }{\alpha H\left(t_{*}\right)}\right]-\frac{\pi }{2}.$$

(90)

Keeping in mind the representation of the complete Bessel function of the imaginary argument (see Section 6.22 (4) of the book [33]).

$$I_0\left(z\right)=\frac{1}{\pi }{\int }_0^{\pi }d\xi e^{zcos\xi },$$

(91)

we could rewrite the integral in (89) as

$${\int }_{{\tau }_1}^{{\tau }_2}d\tau e^{-\frac{3}{2}{h}_{X}cos\tau }={\tilde{I}}_0\left(-\frac{3}{2}{h}_X,{\tau }_2\right)-{\tilde{I}}_0\left(-\frac{3}{2}{h}_X,{\tau }_1\right),{\tilde{I}}_0(z,\tau )=\frac{1}{\pi }{\int }_0^{\tau }d\xi e^{zcos\xi },$$

(92)

using the incomplete Bessel function ${\tilde{I}}_0(z,\tau )$. However, these mathematical details are out from the frame of this article.

3.3.3. Parametric Generation of the Axion Field

Let us return to the perturbed equations for the axion field evolution (see (63) and (64)), but now, we use the periodic Hubble function (83) and the transformation $\psi \to {\left[\frac{a\left(t\right)}{a\left(t_{*}\right)}\right]}^{-\frac{3}{2}}\Psi$ based on the scale factor (82). Now, we deal with the standard Hill equation

$$\ddot{\Psi }=\mathcal{J}\left(t\right)\Psi ,$$

(93)

with the periodic Hill potential

$$\mathcal{J}\left(t\right)=\left[\frac{9}{8}{H}^2\left(t_X\right)-m_A^2\right]-\frac{3}{2}\alpha H\left(t_X\right)sin\left[\alpha (t-t_{*})+\phi \right]+\frac{9}{8}{H}^2\left(t_X\right)cos\left[2\alpha (t-t_{*})+2\phi \right].$$

(94)

The Hill potential is characterized by the period $T=\frac{2\pi }{\alpha }$. The Wronsky determinant for the Equation (93) is equal to a constant. As usual, we work with the reduced dimensionless time variable $\tau =\alpha (t-t_{*})+\phi$ and consider the starting point to be $\tau =0$. In these terms, the Hill potential possesses the period $2\pi$. It is convenient to consider the fundamental solutions ${\Psi }_1\left(\tau \right)$ and ${\Psi }_2\left(\tau \right)$ which satisfy the initial data

$${\Psi }_1\left(0\right)=1,{\Psi }_1^{\prime }\left(0\right)=0,{\Psi }_2\left(0\right)=0,{\Psi }_2^{\prime }\left(0\right)=1,$$

(95)

where the prime denotes the derivative with respect to $\tau$. For such fundamental solutions, the Wronsky determinant is equal to 1. If ${\Psi }_1\left(\tau \right)$ and ${\Psi }_2\left(\tau \right)$ satisfy the Hill equations, the functions ${\Psi }_1(\tau +2\pi )$ and ${\Psi }_2(\tau +2\pi )$ are also the solutions to the Hill equation with the Hill potential possessing the property $\mathcal{J}(\tau +2\pi )=\mathcal{J}\left(\tau \right)$. According to the Floquet theorem [34], we know that for the periodic Hill potential, there exists the so-called normal solution ${\Psi }_{\mathrm{N}}$, which satisfies the rule ${\Psi }_{\mathrm{N}}(\tau +2\pi )=\sigma {\Psi }_{\mathrm{N}}\left(\tau \right)$. The parameter $\sigma$ is the solution to the characteristic equation

$${\sigma }^2-2\mathcal{A}\sigma +1=0\Rightarrow \frac{1}{2}\left(\sigma +\frac{1}{\sigma }\right)=\mathcal{A},$$

(96)

where $2\mathcal{A}$ denotes the trace of the Wronsky matrix calculated at the end of the first period

$$\mathcal{A}=\frac{1}{2}\left[{\Psi }_1\left(2\pi \right)+{\Psi }_2^{\prime }\left(2\pi \right)\right].$$

(97)

Then, we represent the characteristic number $\sigma$ as an exponent $\sigma =e^{2\pi \mu }$ providing the characteristic equation to have the form $cosh2\pi \mu =\mathcal{A}$. Clearly, when $\mathcal{A}>1$, there are two real solutions to the characteristic equation

$${\sigma }_1=\mathcal{A}+\sqrt{{\mathcal{A}}^2-1}=e^{2\pi \mu },{\sigma }_2=\mathcal{A}-\sqrt{{\mathcal{A}}^2-1}=e^{-2\pi \mu },{\sigma }_1{\sigma }_2=1,$$

(98)

$$2\pi \mu =log(\mathcal{A}+\sqrt{{\mathcal{A}}^2-1}).$$

(99)

Finally, the equation ${\Psi }_{\mathrm{N}}(\tau +2\pi )=e^{2\pi \mu }{\Psi }_{\mathrm{N}}\left(\tau \right)$ can be rewritten as

$$e^{-\mu (\tau +2\pi )}{\Psi }_{\mathrm{N}}(\tau +2\pi )=e^{-\mu \tau }{\Psi }_{\mathrm{N}}\left(\tau \right),$$

(100)

and we see that $e^{-\mu \tau }{\Psi }_{\mathrm{N}}\left(\tau \right)$ is a periodic function, and thus ${\Psi }_{\mathrm{N}}\left(\tau \right)=e^{\mu \tau }Y\left(\tau \right)$, where $Y\left(\tau \right)$ denotes a periodic function. Thus, following the idea of the Floquet theorem, we obtain that one of the solutions to the Hill equation grows exponentially, if the characteristic Equation (96) has two different real roots. The periodic function $Y\left(\tau \right)$ can be obtained as follows. We put the function ${\Psi }_{\mathrm{N}}\left(\tau \right)=e^{\mu \tau }Y\left(\tau \right)$ into the Equation (93) and obtain the equation for $Y\left(\tau \right)$ in the form

$$Y^{″}+2\mu Y^{\prime }=\left(\frac{\mathcal{J}}{{\alpha }^2}-{\mu }^2\right)Y.$$

(101)

Then, we search for periodic function $Y\left(\tau \right)$ in the form of trigonometric series

$$Y\left(\tau \right)=f_0+\sum _{k=1}^{\infty }\left[f_{k}sink\tau +g_{k}cosk\tau \right].$$

(102)

It is convenient to use now the auxiliary terms

$$F_0=\frac{9H^2\left(t_X\right)}{8{\alpha }^2}-\frac{m_A^2}{{\alpha }^2}-{\mu }^2,F_1=-\frac{3H\left(t_X\right)}{2\alpha },F_2=\frac{9H^2\left(t_X\right)}{8{\alpha }^2},$$

(103)

in order to formulate the recurrent relations between the coefficients $f_0$, $f_k$, and $g_k$. These recurrent relations contain four subgroups: for $k=0$, $k=1$, $k=2$, and $k\ge 3$:

$$k=0\Rightarrow F_0{f}_0+\frac{1}{2}{F}_1{f}_1+\frac{1}{2}{F}_2{g}_2=0,$$

(104)

$$k=1\Rightarrow f_0{F}_1+f_1\left(1+F_0-\frac{1}{2}{F}_2\right)-2\mu g_1-\frac{1}{2}{F}_1{g}_2+\frac{1}{2}{F}_2{f}_3=0,$$

(105)

$$k=1\Rightarrow g_1\left(1+F_0+\frac{1}{2}{F}_2\right)+2\mu f_1+\frac{1}{2}{F}_1{f}_2+\frac{1}{2}{F}_2{g}_3=0,$$

(106)

$$k=2\Rightarrow f_2\left(4+F_0\right)-4\mu g_2+\frac{1}{2}{F}_1(g_1-g_3)+\frac{1}{2}{F}_2{f}_4=0,$$

(107)

$$k=2\Rightarrow g_2\left(4+F_0\right)+4\mu f_2+F_2{f}_0+\frac{1}{2}{F}_1(f_3-f_1)+\frac{1}{2}{F}_2{g}_4=0,$$

(108)

$$k\ge 3\Rightarrow (k^2+F_0)f_k+2\mu k{g}_k+\frac{1}{2}{F}_1\left(g_{k-1}-g_{k+1}\right)+\frac{1}{2}{F}_2\left(f_{k+2}+f_{k-2}\right)=0,$$

(109)

$$k\ge 3\Rightarrow (k^2+F_0)g_k-2\mu k{f}_k+\frac{1}{2}{F}_1\left(f_{k+1}-f_{k-1}\right)+\frac{1}{2}{F}_2\left(g_{k+2}+g_{k-2}\right)=0.$$

(110)

Clearly, there is the consistent sequential procedure: we extract $g_{k+2}$ from (110), $f_{k+2}$ from (109), $g_4$ from (108), $f_4$ from (107), $g_3$ from (106), $f_3$ from (105), and $g_2$ from (104), thus allowing the coefficients $f_2$, $g_1$, $f_1$, and $f_0$ to remain arbitrary. As an illustration of the instable solution, we can put $g_1=0$, $f_2=0$ and express all the coefficients via the parameters $f_0$ and $f_1$. The corresponding decomposition has the form

$$\begin{gathered} {\Psi }_{\mathrm{N}}\left(\tau \right)=e^{\mu \tau }\left\{f_0+f_{1}sin\tau -\frac{1}{F_2}\left(2F_0{f}_0+F_1{f}_1\right)cos2\tau -\right. \\ \left.-\frac{2}{F_2}\left[F_1{f}_0\left(1+\frac{F_0}{F_2}\right)+f_1\left(2+F_0-\frac{1}{2}{F}_2\right)\right]sin3\tau -\frac{4\mu }{F_2}{f}_{1}cos3\tau +\dots \right\}. \end{gathered}$$

(111)

The coefficients $f_0$ and $f_1$ can be formally found from the initial data

$${\Psi }_{\mathrm{N}}\left(0\right)=f_0+\sum _{k=1}^{\infty }{g}_k,{\Psi }_{\mathrm{N}}^{\prime }\left(0\right)=\mu {\Psi }_{\mathrm{N}}\left(0\right)+\sum _{k=1}^{\infty }k{f}_k.$$

(112)

To conclude, during the episode of the Universe evolution, when the scale factor is the periodic function (82), one of the solutions to the perturbed equation for the axion field behaves as

$$\psi \left(\tau \right)=Y\left(\tau \right)exp\left\{\mu \tau -\frac{3}{2}{h}_{X}sin\tau \right\}.$$

(113)

The envelope function $exp\left\{\mu \tau -\frac{3}{2}{h}_{X}sin\tau \right\}$ is the sinusoidally deformed exponent which grows during the periodic episode of Universe evolution (see Figure 4 for illustration). We deal with the instability of the axion field, which allows us to speak about an “axionization” of the Universe in analogy with the process indicated as scalarization in the work [35].

3.4. Remark about the Post-Inflationary Stage of the Universe Evolution

When we consider the model of aetheric control over the dark fluid evolution, we assume that there exists a critical value of the expansion scalar, say, ${\Theta }_{\ast \ast }=3H\left(t_{\ast \ast }\right)$, so that for $t>t_{\ast \ast }$, the rheological interactions in the dark fluid happen to be switched off, and we deal with the post-inflationary epoch. For this epoch, the dark energy and non-axionic dark matter are assumed to be characterized by the standard equations of state $P=-W$ and $\Pi =0$, respectively. Also, ${\omega }_0=0$, and the dark matter interacts with dark energy via the gravity field only. The conservation laws for the dark energy and non-axionic dark matter become decoupled, and we now have

$$\dot{W}=0,\dot{E}+3HE=0\Rightarrow W=const=W\left(t_{\ast \ast }\right),E\left(t\right)=E\left(t_{\ast \ast }\right){\left[\frac{a\left(t_{\ast \ast }\right)}{a\left(t\right)}\right]}^3.$$

(114)

As for the axion field, we assume that during the period of instability, this field reached the maximal value at $t=t_{\left(\max \right)}$, and at the moment $t=t_{\ast \ast }>t_{\left(\max \right)}$, two parameters, $\varphi \left(t_{\ast \ast }\right)$ and $\dot{\varphi }\left(t_{\ast \ast }\right)$, start to play the roles of the initial data for the master equations in the post-inflationary period. Clearly, we can recover all the specific details of the standard post-inflationary epochs, but here, we mention two interesting cases of the axion field behavior at $t>t_{\ast \ast }$.

3.4.1. Frozen Axion Field with Non-Vanishing Effective Cosmological Constant

If the value $\varphi \left(t_{\ast \ast }\right)$ happens to be close to one of the minima of the axion field potential, i.e., $\varphi \left(t_{\ast \ast }\right)\to n_{\ast \ast }{\Phi }_{*}$ with some integer $n_{\ast \ast }$ and thus $V\to 0$, we obtain the situation analogous to the one described in [36], when the large kinetic energy of axions dominates its potential energy, so the axion field evolves as a stiff matter. Now, one can see that the contribution of the derivative of the potential also is small, and the equation of axion dynamics (33) admits the integral

$$\dot{\varphi }=\dot{\varphi }\left(t_{\ast \ast }\right){\left[\frac{a\left(t_{\ast \ast }\right)}{a\left(t\right)}\right]}^3.$$

(115)

Thus, the key equation of the gravity field takes the form

$$3H^2\left(1+\frac{3}{2}{C}_2\right)=\Lambda +\kappa W\left(t_{\ast \ast }\right)+\frac{1}{2}\kappa {\Psi }_0^2{\left(\dot{\varphi }\left(t_{\ast \ast }\right)\right)}^2{\left[\frac{a\left(t_{\ast \ast }\right)}{a\left(t\right)}\right]}^6+\kappa E\left(t_{\ast \ast }\right){\left[\frac{a\left(t_{\ast \ast }\right)}{a\left(t\right)}\right]}^3.$$

(116)

If the effective cosmological constant $\Lambda +\kappa W\left(t_{\ast \ast }\right)$ is non-vanishing, we can rewrite this equation in terms of auxiliary variable $x=\frac{a\left(t\right)}{a\left(t_{\ast \ast }\right)}$ as

$$H^2=H_{\infty }^2\left[1+2{\mathcal{M}}_1{x}^{-3}+{\mathcal{M}}_2{x}^{-6}\right],$$

(117)

where the following auxiliary parameters are introduced:

$$H_{\infty }^2=\frac{\Lambda +\kappa W\left(t_{\ast \ast }\right)}{3\left(1+\frac{3}{2}{C}_2\right)},2{\mathcal{M}}_1{H}_{\infty }^2=\frac{\kappa E\left(t_{\ast \ast }\right)}{3\left(1+\frac{3}{2}{C}_2\right)},{\mathcal{M}}_2{H}_{\infty }^2=\frac{\kappa {\Psi }_0^2{\left(\dot{\varphi }\left(t_{\ast \ast }\right)\right)}^2}{6\left(1+\frac{3}{2}{C}_2\right)}.$$

(118)

The solution to the equation

$$H_{\infty }(t-t_{\ast \ast })={\int }_1^{\frac{a\left(t\right)}{a\left(t_{\ast \ast }\right)}}\frac{dx}{x\sqrt{1+2{\mathcal{M}}_1{x}^{-3}+{\mathcal{M}}_2{x}^{-6}}}$$

(119)

can be presented in the form

$$\frac{a\left(t\right)}{a\left(t_{\ast \ast }\right)}={\left\{(1+{\mathcal{M}}_1)cosh\left[3H_{\infty }(t-t_{\ast \ast })\right]-{\mathcal{M}}_1)+\sqrt{1+2{\mathcal{M}}_1+{\mathcal{M}}_2)}sinh\left[3H_{\infty }(t-t_{\ast \ast })\right]\right\}}^{\frac{1}{3}}.$$

(120)

In the asymptotic regime $t\to \infty$, we obtain that $H\to H_{\infty }$ and $a\left(t\right)\propto e^{H_{\infty }t}$, i.e., we deal with the de Sitter type asymptote.

3.4.2. Frozen Axion Field with Vanishing Effective Cosmological Constant

When $\Lambda +\kappa W\left(t_{\ast \ast }\right)=0$, we use another re-parametrization of the key equation for the gravity field:

$$H^2={\rho }_1{x}^{-3}+{\rho }_2{x}^{-6},{\rho }_1=\frac{\kappa E\left(t_{\ast \ast }\right)}{3\left(1+\frac{3}{2}{C}_2\right)},{\rho }_2=\frac{\kappa {\Psi }_0^2{\left(\dot{\varphi }\left(t_{\ast \ast }\right)\right)}^2}{6\left(1+\frac{3}{2}{C}_2\right)}.$$

(121)

The solution for the scale factor is now of the form

$$\frac{a\left(t\right)}{a\left(t_{\ast \ast }\right)}={\left[1+3\sqrt{{\rho }_1+{\rho }_2}(t-t_{\ast \ast })+\frac{9}{4}{\rho }_1{(t-t_{\ast \ast })}^2\right]}^{\frac{1}{3}}.$$

(122)

In the asymptotic regime, we obtain the power-law behavior of the scale factor, $a\left(t\right)\propto t^{\frac{2}{3}}$.

4. Discussion and Conclusions

Starting from 1996, cosmologists use specific terminology which is associated with the spontaneous growth of fields and particle numbers in the early Universe. These terms are spontaneous scalarization (see, e.g., [35,37,38]), vectorization [39,40], tensorization [41], and spinorization [42]. Thinking along this line, we pose the question: when and why did the axionization of the Universe occur? In other words, we are interested to elaborate upon the models of spontaneous growth of axions which form the main part of the cosmic dark matter. There are two mechanisms of such spontaneous growth: the first one is associated with axion production due to the decay of gauge fields (see, e.g., [43]); the second (the most known) mechanism is connected with the instability of the axion system. In this work, we consider models of the second type and assume that the origin of the axion instability is connected with the rate of the Universe evolution. We can explain the idea using a few examples. The first example relates to standard inflation, when the scale factor is described by the de Sitter law $a\left(t\right)\propto e^{H_0(t-t_0)}$ with $H\left(t\right)=H_0=const$. Equation (63) is now a differential equation of second order with constant coefficients, and we conclude that both fundamental solutions to this equation do not increase; the model is stable with respect to perturbations. The second example is connected with the power-law behavior of the scale factor $a\left(t\right)\propto t^{\xi }$ (it can be obtained, e.g., for the Universe filled with matter with the equation of state $P=(\Gamma -1)W$). Equation (63) can be now reduced to the Bessel equation, and if $\Gamma \le 2$, the fundamental solutions to this equation also do not increase, so the model is stable. The third example is the super-inflationary solution obtained in this work. As we see from Formula (65) and from Figure 2, there exists the interval of the cosmological time when the axion field increases, such that the axion system should be considered instable. The fourth example is the model with periodic behavior of the scale factor (see (82)). Now, we see from (113) and from Figure 4 that the perturbations of the axion field grow during the episode of periodic evolution of the Universe. Again, we deal with axionic instability, and thus, we can speak about Universe axionization.

To conclude, we have to say that at present, we do not know certainly that the super-inflationary and/or oscillatory periods took place in the history of the early Universe. Systematic analysis of observational data is necessary to clarify this question, and we hope to take part in the detailed discussions. As a first step, one could try to think about imprints that were left by the oscillations of the scale factor in the early-time and late-time stages of Universe evolution. For sure, the main features relate to the axion field growth in the early Universe and formation of the cold dark matter from these relic axions in our epoch. Of course, it would be important to consider a model which unifies these features along the lines elaborated, e.g., in [44,45,46,47]; we hope to present this work in our next papers.

The mechanism of the axion field growth described by the periodic model can be indicated as the parametric mechanism of the Floquet type; this mechanism is well known in nonlinear mechanics. In the work [48], the parametric generation of the scalar field was studied in the case when the oscillations of the scalar field are intrinsic, i.e., they are connected with the properties of the scalar field potential. Our approach is based on the idea that just the oscillations of the gravity field are the cause of parametric instability of the pseudoscalar axion field. As for the model of super-inflation studied in this work, the mechanism of axion instability can be indicated by a generalized parametric model, since instead of the standard Hill equation, we use the generalized Hill equation for the analysis, replacing the trigonometric functions with hyperbolic ones.

We have to emphasize that in all the stages of the Universe evolution, Peccei–Quinn symmetry is assumed to be non-violated; the discrete symmetry $\varphi \to \varphi +m{\Phi }_{*}$ with integer m is supported in our model due to the specific periodic structure of the axion field potential (3).

The last point of our discussion is the status of two exact solutions to the master equations of the presented model: the super-inflationary and the periodic solutions. These exact solutions appear in the model with rheologically active dark fluid consisting of coupled dark energy and non-axionic dark matter. These solutions are marked in the catalogs presented in [9,10] as the ones corresponding to the triple degenerate roots of the characteristic equation. In this paper, we considered both solutions in detail and studied a new composite model in which the super-inflationary behavior of the Universe performs a fine-tuning step, allowing the start of the periodic regime. In this context, we have to emphasize that such a behavior of the expanding Universe is possible due to the aetheric guidance, which is described in our model by inclusion of the expansion scalar $\Theta$ into the axion field potential and into the Lagrangian of the dark fluid. Also, we believe that namely the dynamic aether predetermines the milestone values of the cosmological time $t_0$, $t_{*}$, as well as the time moment $t_{\ast \ast }$, when the periodic episode of the Universe evolution finishes and the standard cosmological expansion starts.