Cosmology from multimeasure multifield model

TL;DR

This paper explores a unified cosmological model that combines inflation, dark energy, dark matter, and the Higgs mechanism, demonstrating through numerical analysis that it can describe the entire Universe evolution with some model adjustments.

Contribution

It reformulates a recently proposed model using effective measures into a quintessential framework and shows its capability to model the full cosmic history with improved dynamics.

Findings

The model can describe the entire Universe evolution seamlessly.

Effective friction in the model leads to finite-time inflation with low slow-roll parameters.

Predicts a constant scalar field in the late Universe.

Abstract

We consider the cosmological application of a (variant of) relatively newly proposed model \cite{1609.06915} unifying inflation, dark energy, dark matter, and the Higgs mechanism. The model was originally defined using additional non-Riemannian measures but it can be reformulated into effective quintessential model unifying inflation, dark energy and dark matter. Here we demonstrate numerically that it is capable of describing the entire Universe evolution in a seamless way, but this requires some revision of the model setup. The main reason is that there is a strong effective friction in the model, a feature which has been neglected in the pioneer work. This improves the model potential for proper description of the Universe evolution, because the friction ensures a finite time inflation with dynamically maintained low-value slow-roll parameters in the realistic scenarios. In addition, the model predicts the existence of a constant scalar field in late Universe.

Equations

$$S=S_{darkon}+S_{inflaton}.$$

$$S_{darkon}$$

$$S_{inflaton}$$

$$L^{(0)}$$

$$L^{(1)}$$

$$L^{(2)}$$

$$\Phi (Z)=\frac{1}{3}{ϵ}^{\mu \nu \kappa \lambda }{\partial }_{\mu }{Z}_{\nu \kappa \lambda }.$$

$$0=\frac{\delta S^{\prime }}{\delta {\hat{z}}^{\mu }}=\frac{\delta \int d^{4}x({\partial }_{\nu }{\hat{z}}^{\nu })L^{\prime }}{\delta {\hat{z}}^{\mu }}\Rightarrow {\partial }_{\mu }{L}^{\prime }=0\forall \mu .$$

$$L^{(0)}$$

$$R+L^{(1)}$$

$$L^{(2)}+\frac{\Phi (\mathcal{H})}{\sqrt{-g}}$$

$$\frac{\Phi (\mathcal{B})}{\sqrt{-g}}$$

$${\Gamma }_{\lambda \kappa }^{\mu }=\frac{1}{2}{\tilde{g}}^{\mu \nu }({\partial }_{\nu }{\tilde{g}}_{\lambda \kappa }-{\partial }_{\lambda }{\tilde{g}}_{\nu \kappa }-{\partial }_{\kappa }{\tilde{g}}_{\nu \lambda })$$

$${\tilde{g}}_{\mu \nu }={\chi }_1{g}_{\mu \nu },$$

$${\chi }_1=\frac{\Phi (\mathcal{A})}{\sqrt{-g}}.$$

$$S^{(eff)}=\int d^{4}x\sqrt{-\tilde{g}}(\tilde{R}+L^{(eff)})$$

$${\tilde{R}}_{\mu \nu }-\frac{1}{2}{\tilde{g}}_{\mu \nu }\tilde{R}=\frac{1}{2}{T}_{\mu \nu }^{(eff)}$$

$$0$$

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }=\frac{1}{2}{T}_{\mu \nu }^{(1)}+\frac{{\chi }_2}{2{\chi }_1}{T}_{\mu \nu }^{(2)}+\frac{1+{\chi }_3}{2{\chi }_1}{T}_{\mu \nu }^{(0)}+\left(-\frac{M_1}{2}+\frac{{\chi }_2{M}_2}{2{\chi }_1}+\frac{{\chi }_3{M}_0}{{\chi }_1}\right)g_{\mu \nu }.$$

$${\chi }_1{L}^{(eff)}$$

$$u\to \tilde{u}$$

$$\tilde{Y}$$

$$\tilde{X}$$

$$\tilde{Y}=1/{\chi }_1.$$

$$L^{(eff)}=\tilde{X}-\tilde{Y}(V+M_1-{\chi }_2{b}_0{e}^{-\alpha \varphi }\tilde{X})+{\tilde{Y}}^2({\chi }_2(U+M_2)-2M_0).$$

$${\partial }_{\mu }\left(\sqrt{-\tilde{g}}{\tilde{g}}^{\mu \nu }{\partial }_{\nu }\tilde{u}\frac{\partial L^{(eff)}}{\partial \tilde{Y}}\right)=0.$$

$$\frac{1}{{\chi }_1}=\frac{2(L^{(1)}+M_1-T^{(1)}/2)}{{\chi }_2(T^{(2)}+4M_2)+(1+{\chi }_3)T^{(0)}+8{\chi }_3{M}_0}$$

$${\tilde{Y}}^{static}=\frac{V+M_1}{2{\chi }_2(U+M_2)-4M_0}.$$

$$U_{eff}=\frac{1}{4}\frac{(V+M_1{)}^2}{{\chi }_2(U+M_2)-2M_0}.$$

$$S^{(eff)}{|}_{FLRW}$$

$$v\left(V+M_1-\frac{1}{2}{\chi }_2{b}_0{e}^{-\alpha \varphi }{\dot{\varphi }}^2\right)+v^3({\chi }_2(U(\varphi )+M_2)-2M_0)=\frac{p_u}{a(t{)}^3}$$

$$\frac{d}{dt}\left(a(t{)}^3\dot{\varphi }(1+\frac{v^2}{2}{\chi }_2{b}_0{e}^{-\alpha \varphi })\right)+a(t{)}^3(\alpha \frac{v^2}{4}{\chi }_2{b}_0{e}^{-\alpha \varphi }{\dot{\varphi }}^2+\frac{v^2}{2}{V}_{\varphi }-{\chi }_2\frac{v^4}{4}{U}_{\varphi })=0$$

$$\frac{\dot{a}(t{)}^2}{a(t{)}^2}$$

$$\frac{\ddot{a}(t)}{a(t)}$$

$$p$$

Full text

1 Introduction

According to the current understanding, we live in an expanding Universe. The rate of this expansion, however, is a hot topic for the community due to the so called $H_{0}$ tension – i.e. the difference between the measurements of the Hubble constant locally (type IA supernova) and from the cosmic microwave background (the Planck mission). A recent study [2] has shown that the $H_{0}$ tension is even more significant than expected – the tension between the model-independent local measurements and the $\Lambda-CDM$ inferred Plank CMB measurements is about $4.4\sigma$. Such a discrepancy calls for a revision of the default $\Lambda-CDM$ model and opens the door for theories of modified gravity.

It is known that the Universe has passed trough different stages – radiation-dominated era, matter-dominated era and current, dark-energy dominated era. Those epochs can be described well using the $\Lambda-CDM$ model as different components in the energy density of the Universe. Because of this, any new theory of gravity, needs to describe those stages just as well. Additionally, however, any new theory of gravity, has to be able to describe the initial inflation of the Universe – the exponential expansion that is designed to solve the unresolved problems in cosmology – the horizon problem, the flatness problem, the missing monopols problem and the large-structures formation problem [3]). Numerous such theories have been proposed so far – chaotic inflation (one scalar field rolling in a potential), multi-field theories (more than one scalar fields), modified gravity (f(R), Brans-Dicke etc), for a review see [4, 5] , but they have seen their strongest limitation so far in the gravitational waves astronomy (for a review see [6]). Among the few surviving theories is the k-essence type of theories of modified gravity, which do not modify the speed of the gravitational waves and thus are not in violation with the observational data.

The considered in this article model is in the framework of multi-measure gravitational models. The basic idea is [7] that on the space-time manifold can be defined a metric-independent (non-Riemannian) volume form (or forms) in terms of an auxiliary antisymmetric gauge field(s) of maximal rank. Such forms can be used in the action integral to replace the usual integration measure (or to participate together with it). The approach was further developed in Refs.[8, 9, 10, 11] and successfully applied to define unified dark energy – dark matter model [12, 13], unified inflation – dark energy model [15, 14], and a model unifying inflation, dark energy, dark matter, and Higgs mechanism [1]. Some of these models are naturally connected to the $f(R)$ gravity and to the quintessential inflation proposed by Peebles and Vilenkin [16]. Effectively, the role of any such non-Riemannian measure is to produce very simple but powerful system of equations of motion for the quantity coupled to it in the Lagrangian. As a result on-shell the corresponding quantity is a (dynamically generated) constant, determined by the problem boundary conditions.

In Ref. [17] we have studied the cosmological aspects of a model [13] with a single scalar field (called darkon) coupled to two independent volume-forms. The model predicts both the existence of dark energy (as a dynamically generated cosmological constant) and of dark matter (as a dust contribution to the energy-momentum tensor). The model is capable to fit the Supernova Type Ia data still leaving a large freedom in its parameters. Here we consider an “upgrade” of the darkon model aiming to include initial inflation. The new model features two scalar fields (darkon and inflaton) and has been introduced in Ref.[1]. In this model, the inflaton moves over an effective step-like potential (which plays the role of cosmological constant) with two infinite plateaus connected with a steep slope.

It is advocated in Refs. [14, 1] that the left, higher plateau corresponds to the inflationary Universe, while the right, lower plateau – to the current expanding Universe. Our results show that for some particular values of the parameters we can obtain all epochs in the evolution of the Universe, including a graceful exit to current exponential expansion. However, the results demonstrate the existence and importance of a friction term in the inflaton equation of motion. As a consequences it is impossible to start the evolution from the left plateau and to reach the right plateau regardless of the initial inflaton velocity. The early inflation, instead, is generated on the slope of the effective potential after a brief ultra-relativistic period due to the initial singularity. While the large parameter-space has not been fully explored, the cases presented here are able to demonstrate the qualitatively different behaviors of the numerical solution.

2 The model

We work with a simplified version (without Higgs and gauge fields) of the model proposed in Ref. [1]. The model itself is based on earlier results on dark matter – dark energy unified model [12, 13] and inflation – dark energy unified model [15, 14]. We have two scalar fields — the darkon $u$ and inflaton $\phi$ and the action of the is the following:

[TABLE]

Here (in units $G_{N}=1/16\pi,\;c=1$):

[TABLE]

We use the notations

[TABLE]

where $Y=-\frac{1}{2}g^{\mu\nu}\partial_{\mu}u\partial_{\nu}u$ and $X=-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$ are the darkon and inflaton kinetic terms respectively, and $V(\phi)=f_{1}e^{-\alpha\phi},\;U(\phi)=f_{2}e^{-2\alpha\phi}$ and $W(u)$ are corresponding potentials. This particular form of the inflaton potentials ensures a global Weyl invariance of $S_{inflaton}$. On the other hand it turns out that the equations of motion do not depend on the form of the darkon potential $W(u)$ so we do not specify it. Note that we use $R=g^{\mu\nu}R_{\mu\nu}(\Gamma)$ written in Palatini formalism where $\Gamma^{\lambda}_{\mu\nu}$ are treated as a priori independent variables.

The quantity $\Phi(Z)$, where $Z=\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{H}$, which appears in Eqs. 2, 3 is an additional non-Riemannian generally covariant integration measure defined with the help of an auxiliary 3-index antisymmetric tensor gauge field $Z_{\nu\kappa\lambda}$:

[TABLE]

Thus the model depends on $4$ parameters — dimensionless $b_{0}$ and dimensional $f_{1},f_{2}$ and $\alpha$. The action in Eq.3 is exactly the same as that in Refs. [15, 14], while the action in Eq. 2 differs from that in [13] which is $S=\int{d^{4}x\;\sqrt{-g}(R-\epsilon R^{2})}+S_{darkon}$ because the first terms in the latter already present in $S_{inflaton}$ in a hidden form.

2.1 The model equations of motion

Our primary objective is to use the equations of motion for the auxiliary fields to eliminate them and to obtain an effective Lagrangian written down entirely in terms of metric and scalar fields. Probably the only unfamiliar equations of motion are those obtained by variation of the action with respect to the 3-forms $Z=\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{H}$. The easiest way to get them is to note that any such $Z$ has only 4 independent components, namely $\hat{z}^{\mu}=\epsilon^{\mu\nu\kappa\lambda}Z_{\nu\kappa\lambda}/3$. Written in terms of $\hat{z}^{\mu}$ the integration measure is $\Phi(Z)=\partial_{\mu}\hat{z}^{\mu}$ and now it is trivial to obtain the variation of the action $S^{\prime}=\int{d^{4}x\;\Phi(Z)L^{\prime}}$ (this is the most general form of the action with the non-Riemannian measures we use) with respect to independent components of $Z$:

[TABLE]

Thus, the solution of the $Z$-equation of motion is $L^{\prime}=\mathrm{constant}$. We have four auxiliary 3-forms in Eq.1, so we obtain four dynamically generated integration constants:

[TABLE]

Note that in the above equations only $\chi_{2}$ is dimensionless, while $M_{0},M_{1}$ and $M_{2}$ are with dimension $\mathrm{mass}^{4}$, the same as $f_{1}$ and $f_{2}$.

The full-fledged derivation of the equation of motion for $\Gamma^{\alpha}_{\beta_{\gamma}}$ and its solution can be found in Ref.[8] . The result is that $\Gamma$ is the Levi–Civita connection

[TABLE]

for the Weyl rescaled metric $\tilde{g}_{\mu\nu}$:

[TABLE]

where

[TABLE]

There is a simple way in two steps not to derive but to anticipate the result given in Eq. 13. First, note that from Eq.14 we get $\Phi(\mathcal{A})g^{\mu\nu}=\sqrt{-\tilde{g}}\tilde{g}^{\mu\nu}$, so the term $\int d^{4}x\;\Phi(\mathcal{A})R=\int d^{4}x\;\Phi(\mathcal{A})g^{\mu\nu}R_{\mu\nu}(\Gamma)$ in Eq.3 looks exactly as the Einstein–Hilbert action in Palatini formalism with metric $\tilde{g}_{\mu\nu}$. Second, under some general assumptions, the solution of Palatini formalism is the corresponding Levi-Civita connection Eq.13 Q.E.D.

Eq.13 signals that the effective Lagrangian we are looking for has to be written in terms of $\tilde{g}_{\mu\nu}$. More precisely, we are looking for $L^{(eff)}$ and effective action

[TABLE]

such that the corresponding Einstein equation

[TABLE]

is equivalent to the equation of motion obtained by varying action Eq. 1 with respect to $g^{\mu\nu}$, which is:

[TABLE]

Here $\tilde{R}_{\mu\nu}$ and $\tilde{R}=\tilde{g}^{\mu\nu}\tilde{R}_{\mu\nu}$ are the Ricci tensor and scalar respectively for metric $\tilde{g}_{\mu\nu}$, $T^{(eff)}_{\mu\nu}\equiv\tilde{g}_{\mu\nu}L^{(eff)}-2\partial L^{(eff)}/\partial\tilde{g}^{\mu\nu}$, and $T^{(i)}_{\mu\nu}\equiv g_{\mu\nu}L^{(i)}-2\partial L^{(i)}/\partial g^{\mu\nu},\;\;i=0,1,2$.

Using Eqs.10, 11, and the notation $\chi_{3}=\Phi(\mathcal{A})/\sqrt{-g}$, we put Eq.18 in the following form:

[TABLE]

We need two more relations, both following from Eq.15, in order to obtain $L^{(eff)}$. First, because by definition $\tilde{R}_{\mu\nu}=R_{\mu\nu}$, so, $\tilde{R}=R/\chi_{1}$ and therefore $\tilde{R}_{\mu\nu}-\frac{1}{2}\tilde{g}_{\mu\nu}\tilde{R}=R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}$. Second, $(\tilde{g}_{\mu\nu}-2\partial/\partial\tilde{g}^{\mu\nu})L^{(eff)}=\chi_{1}(g_{\mu\nu}-2\partial/\partial g^{\mu\nu})L^{(eff)}$. Then, it is straightforward to check that $L^{(eff)}$ written in terms of $g_{\mu\nu}$ is;

[TABLE]

where we have used Eq.9 to obtain the final form. In order to write $L^{(eff)}$ as a function depending on $\tilde{g}_{\mu\nu}$ we introduce the notations

[TABLE]

The change of variable Eq. 21 has been used in Refs.[12, 13] to demonstrate that the darkon model represents relativistic fluid. An important consequence of Eqs. 15, 21, 22 is that

[TABLE]

Taking into account all this, we obtain

[TABLE]

This Lagrangian is non-linear with respect to both scalar fields kinetic terms and thus can be classified to be of generalized k-essence type.

The effective Lagrangian possesses an obvious global symmetry $\tilde{u}\rightarrow\tilde{u}+w$ which produces the following current conservation

[TABLE]

which coincides with the equation of motion for $\tilde{u}$.

The variation with respect to the scalar field $\phi$ is easier done in a specific metric and we postpone it to subsection 2.3.

2.2 The static case

Taking the trace of Eq.18 we get:

[TABLE]

which in the case when the fields $u$ and $\phi$ are static and taking into account Eq.24 gives:

[TABLE]

Substituting the above equation into Eq.25 we see that the effective Lagrangian determines the following effective potential:

[TABLE]

This potential for certain parameters looks like two infinite plateaus with asymptotic values $U_{-}=U_{eff}|_{\phi\rightarrow-\infty}=f_{1}^{2}/4\chi_{2}f_{2}$ and $U_{+}=U_{eff}|_{\phi\rightarrow+\infty}=M_{1}^{2}/(4\chi_{2}M_{2}-8M_{0})$ and a steep slope connecting them. This means that the model naturally determines two different cosmological constants when $\phi\to\pm\infty$. However, there is a catch in the usage of the effective potential — the effective potential does not bring the kinetic energy in standard form. We write more on this topic at the end of next subsection.

2.3 The model in the Friedmann – Lemaître – Robertson – Walker (FLRW) metric

In the metric $\tilde{g}_{\mu\nu}=\mathrm{diag}\{-1,a(t)^{2},a(t)^{2},a(t)^{2}\}$ and with fields depending only on the time coordinate $(t)$, Eq.16 takes the form

[TABLE]

Here we use the notation $v=d\tilde{u}/dt$) and have performed the integration by parts when writing the term $\int{d^{4}x\;\sqrt{-\tilde{g}}\tilde{R}}$.

The equations of motion and their solutions are relatively simple in the FLRW metric. From Eq. 26 we obtain the following cubic equation with respect to $v$:

[TABLE]

where $p_{u}$ is an integration constant. The $\phi$-equation of motion is:

[TABLE]

where $U_{\phi}=\frac{\partial U}{\partial\phi}$, $V_{\phi}=\frac{\partial V}{\partial\phi}$.

The Friedmann equations for a perfect fluid with density $\rho$ and pressure $p$ are:

[TABLE]

where

[TABLE]

We will use only the first Friedmann equation (Eq. 33), because the second one (Eq. 34) is a consequence of Eqs. 31, 32 and 33.

Note that the model in FLRW metric has a strong friction term. To see this, one can decouple Eqs. 31–33 assuming that the system is in exponentially growing regime $a(t)=e^{Ht}$. In this case Eq. 32 can be put in the form:

[TABLE]

with a general solution

[TABLE]

We see that the friction is entirely determined by the Hubble parameter $H$. As a result one cannot start from the left plateau and reach the right plateau of the effective potential a feature we also observe numerically. The friction ensures that in the far future the kinetic energy of the $\phi$ field vanishes and we have a dynamically generated cosmological constant

[TABLE]

As it is seen form Eq. 37 the existence of friction is independent of the model parameters and it is a feature of FLRW metric.

It has been already mentioned the model possesses an effective potential depending only on the inflaton field (see Eq. 29). Unfortunately the effective potential does not bring the kinetic energy in standard form. For example in the limit of slow roll approximation (neglecting the terms $\sim\;\dot{\phi}^{2},\;\;\dot{\phi}^{3},\;\;\dot{\phi}^{4}$) the inflaton equation has the form:

[TABLE]

where $A=v(t)^{2}b\,\chi_{2}e^{-\alpha\phi}/2$ and $H=\dot{a}(t)/a(t)$. One can see that putting the kinetic term in standard form (i.e., dividing Eq. 39 by $(1+A)$) introduces a singularity at $a=0$ in the effective potential term (due to $v\sim p_{u}/a^{3}$). This changes significantly the behavior of the system in the region $a\sim 0$ as seen in Ref.[21] . Because of this problem, we will not use the effective potential, but will stick in our calculations to the system Eqs. 31–33.

3 Numerical results

Putting all things together, the algebraic–differential system we have to solve is:

[TABLE]

where the parameters of the cubic equation are:

$\mathbf{a}=-(V+M_{1}-\chi_{2}b_{0}e^{-\alpha\phi}\dot{\phi}^{2}/2)/3(\chi_{2}(U+M_{2})-2M_{0})$, $\mathbf{b}=-p_{u}/2a(t)^{3}(\chi_{2}(U+M_{2})-2M_{0})$.

Note that Eq.(40) always has at least one real root but there is no global continuous definition of it on $\{\mathbf{a},\mathbf{b}\}$ plane. The root we use in this article is:

[TABLE]

and a special care is taken that this root is real in the domain where we use it.

The system Eqs.(40–42) is of second order with respect to $\phi$ and we shall integrate it numerically111For the numerical solution of Eqs.(41–42) we use Fehlberg fourth-fifth order Runge-Kutta method with degree four interpolation implemented in Maple.. In order to put it on a computer we have to fix our units. We use $c=1$, $G_{N}=1/16\pi$ and $t_{u}=1$, where $c$ is the speed of light, $G_{N}$ is Newton constant and $t_{u}$ is the present day age of Universe ($t_{u}=13.810^{9}\;y$). Thus, our mass unit is $1.62\times 10^{59}M_{Pl}$ and the cosmological constant ($\Lambda=1.910^{-35}s^{-2}$) is $\Lambda\approx 3.6$. The modern-day Hubble constant in these units equals to $H_{0}\approx 0.95$. Also, in these units, the three periods of the evolution of the universe are : exit from early inflation: $t_{\mathrm{SD}}\sim 10^{-50}$ and start of present day acceleration $t_{\mathrm{AE}}\sim 0.71$.

Altogether in our system we have 4 free parameters, namely $\{\alpha,b_{0},f_{1},f_{2}\}$, 5 integration constants $\{M_{0},M_{1},M_{2},\chi_{2},p_{u}\}$ and three initial conditions $\{a(0),\phi(0),\dot{\phi}(0)\}$ 222Note that the differential equation Eq. 41 has a singular point in ${a}(0)=0$ due to the term $\sim 1/a(t)^{3}$ in $\rho$. So, a natural replacement of the initial condition $a(0)=0$ is the normalization condition $a(1)=1$, where $t=1$ corresponds to the current moment.. We take some steps to constrain the used parameters:

– In order to ensure that we are working with a real root of Eq.(40) , we choose $M_{0}<0$ (i.e. $\mathbf{b}<0$). Moreover, we use $M_{0}\sim-0.01$. We are free to fix $M_{0}$ in this way, because as we have shown in Ref. [17] in the asymptotic darkon-dominated late-time inflation universe, we have a nice fit of Supernovae Ia data for any $M_{0}$.

– In this article we work with $\dot{\phi}(0)=0$. In the general case, due to the strong friction, the system is not very sensitive to changes in $\dot{\phi}(0)$.

– We use $b_{0}>0$ to avoid problems with $\rho$, because $b_{0}<0$ can lead to negative kinetic term in $\rho$.

– It is assumed in Ref.[14] that the left plateau of the effective potential corresponds to the pre-inflationary Universe (Planck times) and the right plateau to the current and future accelerated expansion. Since we do not observe a solution for the inflaton connecting the two plateaus, the only restriction we apply to the model is that the left plateau is higher than the right one, i.e.

[TABLE]

Moreover, the estimates in Ref.[14] for the values of the ratio $f_{1}^{2}/(\chi_{2}f_{2})$ have to be considered as lower limit.

We have presented already [18] some results from the numerical investigation of the model. Here we focus our attention on some unexpected features of the model which may be useful for its further theoretical improvement. We observe that the solution of the system Eqs. 40–42 critically depends on the initial value of the inflaton field. If $\phi(0)$ is out of some very specific parameter dependent interval (which is always within the slope of effective potential), we obtain two-stage evolution of the Universe beginning with short deceleration followed by infinite acceleration. When $\phi(0)$ is in the interval we have four stage evolution which we refer to as“physically realistic”. We have a short first deceleration epoch (FD), early inflation (EI), second deceleration (SD) which we interpret as radiation and matter dominated epochs together and finally we have infinite accelerating expansion (AE). There are also some very specific ”pathological” parameter configurations with only deceleration. The physically realistic parameter configurations are just a small part of the parameter space, very sensitive to fine-tuning. Thus, a significant part of the numerical work is to find the points in the parameter space for which we observe the physically realistic behavior. This means that we have to find solutions for which the second derivative of the scale factor crosses the t-axis in three points ($t_{EI},t_{SD}$ and $t_{AE}$), turning the problem into root-finding problem. This can be done by plotting the function or if we require higher precision, it can be done by root-finding algorithms such as the 1-d Muller algorithm applied to the polynomial approximation of $\ddot{a}(t)$.

In order to illustrate the physically realistic evolution we consider two different and hopefully representative regions of the parameter-space. In the first case the dimensionless parameter $\chi_{2}\sim 1$ ($M_{2}<<|M_{0}|$) and in the second one $\chi_{2}<<1$ ($M_{2}\gg|M_{0}|$). The exact values of all other parameters are given in the Fig.1c caption. We have used the parameter $b_{0}$ to set $t_{\mathrm{AE}}\sim 0.71$ and parameters $f_{1}$ and $p_{u}$ to ensure $a(1)=1$.

The two different region of exponential expansion of the Universe are well visible on the logarithmic plot of $a(t)$ Fig 1c panel (a).

The other phases of the Universe evolution are easily traceable on Fig. 1c, panel (b) where the time dependence of $\ddot{a}(t)$ is plotted, and also from panel (d) where the equation of state (EOS) $w=p/\rho$ of the universe is presented. We can see that the Universe passes trough the following stages: First, beginning at $t=0$ we observe the EOS of ultra-relativistic matter with $w=1/3$. The existence of this phase doesn’t contradict to the observations, because currently we have information solely from the time after the initial inflation. Second, after the moment $t_{EI}$ we have inflation with EOS of dark energy $w\to-1$. Third, after the moment $t_{SD}$ the Universe is in matter domination stage where $w>-1/3$ and $w\to 0$. Finally, after the moment $t_{AE}$ we have slow accelerated expansion with $w<-1/3$. These numerical results coincide with the theoretically calculated asymptotic for $w$, namely $w|_{a(t)=0}\to 1/3$ and $w|_{a(t)=\infty}=-1$.

The friction in the model can be observed on Fig. 1c, panels (c), where it is seen that the inflaton scalar field $\phi$ tends to a constant. This means that the theory predicts the existence of a scalar field with in general nonzero average value in the late Universe.

On the final two panels (e) and (f) of Fig. 1c we present the logarithm of the Hubble parameter $H=\dot{a}(t)/a(t)$ and the “slow-roll” parameters of the model [20] defined as $\epsilon=-\dot{H}/H^{2}$ and $\eta=-\ddot{\phi}/H\dot{\phi}$. We observed that both $\eta$ and $\epsilon$ are $<|1|$ during inflation even though, they are not much smaller than $1$ as expected. The parameter $\eta$ is $\eta\sim 0.1$ during the inflation in both cases, but at the beginning of the accelerated expansion $\eta<1$ only in Case 2. We can attribute this to the fact that during the inflation in Case 1 we have $A(\phi)>>1$, while in Case 2 $A(\phi)\to 1$ (see Eq. 39). In Case 2, the parameter $\epsilon_{EI}\to 0$ and $\epsilon_{AE}<0.5$, while in Case 1, $\epsilon_{EI}\to 0.1$ and $\epsilon_{AE}\sim 1$. Since $\epsilon$ measures the ”exponentiality” of $a(t)$, this shows that in Case 2 the dynamics is dominated entirely by inflaton in nearly constant potential, while in Case 1 the dynamics of the inflaton and darkon are not decoupled. Physically, we attribute the observed difference mainly to the more complex movement of the inflaton in this case — it initially moves to smaller values, then at $t\sim 0.05$ it turns back and gradually increases its value (see Fig. 1c (c)). Thus for relatively long period of time $U_{eff}(\phi(t))$ does not change significantly which together with the small value of $p_{u}$ ensures almost pure exponential growth of $a(t)$.

The parameter $\alpha$ is crucial for the length of both inflation and matter domination epochs. We have investigated the model behavior with respect to small parameter variation as well. In both cases, higher value of $\alpha$ leads to a lower $t_{SD}$ and therefore longer matter domination. However, it turns out that it is impossible to obtain realistic $t_{\mathrm{SD}}$ solely by increasing $\alpha$. Two other parameters $f_{1}$ and $f_{2}$ affect significantly the value of $t_{SD}$. When they are increased we get smaller $t_{SD}$. As a result, we expect $\alpha,\;\;f_{1}$ and $f_{2}$ to have really big values for realistic $t_{SD}$. This is just the opposite to the estimates in Ref.[14] where it is asserted that in order for the theory to be compatible with Planck data [19], one needs $\alpha\to 0$ and in addition $f_{1}\sim f_{2}\ll 1$. From our experiments shown on the plots, the smaller $\alpha$ (Case 2) leads to significantly longer early inflation, thus leading to very short matter-domination epoch.

The e-folds parameter related to the inflation is defined as $N=\ln{(a_{SD}/a_{EI})}$ where the subscripts denotes the beginning ($(EI)$) and end ($(SD)$) of the inflation. The theoretical estimation for the number of e-folds needed to solve the horizon problem is model-dependent but is $N>70$. One such estimation can be found in Ref. [22] . In our cases, $t_{EI}=0.013,t_{SD}=0.349$ (Case 1) and $t_{EI}=0.014,t_{SD}=0.663$ (Case 2), and thus $N^{Case1}\approx 9$ and $N^{Case2}\approx 18$. The higher e-folds number in Case 2 is easily explained by the much longer inflation in this case. For the Hubble parameter, we obtain at $t=1$, in the two cases $H_{0}=\{1.2,2.4\}$ respectively. One can see that in the first case, the modern-day Hubble parameter is much closer to the theoretically expected value, due to the fact that in the first case, the inflaton is lower on the right plateau of the effective potential, thus $U_{eff}\to U_{+}$. In the second case, we are a bit higher on the potential, which affects the value of $H$. From the numerical experiments it seems that obtaining better Hubble and e-folds parameters is a matter of extensive numerical study of the full parameter space, which we have not performed yet.

4 Conclusions

We have presented the numerical study of the model of Guendelman-Nissimov-Pacheva aiming to describe inflation, dark energy, dark matter in a seamless way trough the evolution of two scalar fields – the darkon and the inflaton.

Our calculations have confirmed that qualitatively this model can successfully describe the evolution of Universe, since it naturally reproduces its main phases. The values of the parameters we use, however, differ significantly from the ones estimated in Refs. [14, 1], which poses important questions related to the nature of those parameters.

Furthermore, we have shown that there is a strong friction term in the system, due to which, the inflaton field tends to a constant. This means that the model predicts the existence of a non-zero averaged scalar field in the current Universe.

Also, we have shown that the inflationary region of the effective potential is its slope, while the evolution starting from higher plateau is nonphysical as it leads to eternal inflation. As a result, the duration of the inflation in the “realistic” Universe evolution is always finite, contrary to the case in Refs.[14, 1] where the inflation can be arbitrary long.

Finally, an interesting feature is that for certain choices of the parameters and the initial conditions of the system, the inflaton field increases its value in the beginning of the integration. This hints of a dynamical mechanism for enhancing the initial exponential expansion by holding the inflaton near its starting position on the effective potential.

Acknowledgments

It is a pleasure to thank E. Nissimov and S. Pacheva for the discussions. The work is supported by BAS contract DFNP – 49/21.04.2016, by Bulgarian NSF grant DN-18/1/10.12.2017 and by Bulgarian NSF grant 8-17.