Dynamical Analysis of an $n-H-T$ Cosmological Quintessence Real Gas Model with a General Equation of State

TL;DR

This paper analyzes the cosmological dynamics of a quintessence model based on real gas with a general equation of state, revealing stable cyclic solutions and conditions for inflation and deflation.

Contribution

It introduces a dynamical system framework for real gas quintessence models, finds global first integrals, and demonstrates the existence of cyclic universe solutions.

Findings

Existence of stable periodic (cyclic) solutions under certain conditions.

Identification of conserved quantities in the system.

Demonstration of inflation and deflation regions in cyclic solutions.

Abstract

The cosmological dynamics of a quintessence model based on real gas with general equation of state is presented within the framework of a three-dimensional dynamical system describing the time evolution of the number density, the Hubble parameter, and the temperature. Two global first integrals are found and examples for gas with virial expansion and van der Waals gas are presented. The van der Waals system is completely integrable. In addition to the unbounded trajectories, stemming from the presence of the conserved quantities, stable periodic solutions (closed orbits) also exist under certain conditions and these represent models of a cyclic Universe. The cyclic solutions exhibit regions characterised by inflation and deflation, while the open trajectories are characterised by inflation in a "fly-by" near an unstable critical point.

Equations

$$\displaystyle ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=c^{2}dt^{2}-a^{2}(t)\Bigl{[}\frac{dr^{2}}{1-kr^{2}}+r^{2}(d\theta^{2}+\sin^{2}\theta\,d\phi^{2})\Bigr{]},$$

$$G_{\mu \nu }+\Lambda g_{\mu \nu }=\kappa T_{\mu \nu },$$

$$T_{\mu \nu }=(\tilde{\rho }+\tilde{p})u_{\mu }{u}_{\nu }-\tilde{p}{g}_{\mu \nu }.$$

$$\ddot{a}$$

$${\dot{a}}^2$$

$$\frac{\ddot{a}}{a}$$

$$H^2$$

$$\dot{H}=-\frac{1}{2}({\rho }_d+\rho +p),$$

$$\frac{\partial n}{\partial t}+\vec{\nabla }\cdot \left(n\vec{v}\right)=0,$$

$$\displaystyle u^{\mu}=(1-v^{2})^{-\frac{1}{2}}\,\,\Bigl{(}1,{d\vec{x}\over dt}\Bigr{)}=(1-v^{2})^{-\frac{1}{2}}\,\,(1,\vec{v})=(1-v^{2})^{-\frac{1}{2}}\,\,\frac{dt}{d\tau}$$

$$\dot{n}+3Hn=0,$$

$$\dot{\rho }+3H(\rho +p)=0$$

$$\dot{{\rho }_d}+3H{\rho }_d=0.$$

$$\dot{H}=-\frac{3}{2}{H}^2-\frac{1}{2}p.$$

$$\displaystyle dS=\frac{1}{T}\,\,d\Bigl{(}{\rho\over n}\Bigr{)}+\frac{p}{T}\,\,d\Bigl{(}{1\over n}\Bigr{)}=-\left({\rho+p\over Tn^{2}}\right)\,\,dn+{1\over Tn}\,\,d\rho.$$

$$\displaystyle\biggl{[}\frac{\partial}{\partial n}\biggl{(}\frac{\partial S}{\partial\rho}\biggr{)}_{n}\biggr{]}_{\rho}=\biggl{[}\frac{\partial}{\partial\rho}\biggl{(}\frac{\partial S}{\partial n}\biggr{)}_{\rho}\biggr{]}_{n}\qquad\mbox{or}\qquad\biggl{[}\frac{\partial}{\partial n}\biggl{(}\frac{1}{Tn}\biggr{)}\biggr{]}_{\rho}=\biggl{[}\frac{\partial}{\partial\rho}\biggl{(}-\frac{\rho+p}{Tn^{2}}\biggr{)}\biggr{]}_{n}.$$

$$\displaystyle n\biggl{(}\frac{\partial T}{\partial n}\biggr{)}_{\rho}+(\rho+p)\biggl{(}\frac{\partial T}{\partial\rho}\biggr{)}_{n}=T\biggl{(}\frac{\partial p}{\partial\rho}\biggr{)}_{n}.$$

$$\displaystyle\rho+p=T\biggl{(}\frac{\partial p}{\partial T}\biggr{)}_{n}+n\biggl{(}\frac{\partial\rho}{\partial n}\biggr{)}_{T}.$$

$$\displaystyle\biggl{(}\frac{\partial T}{\partial n}\biggr{)}_{\rho}\dot{n}+\biggl{(}\frac{\partial T}{\partial\rho}\biggr{)}_{n}\,\,\dot{\rho}=-3HT\biggl{(}\frac{\partial p}{\partial\rho}\biggr{)}_{n},$$

$$\displaystyle\dot{T}=-3HT\biggl{(}\frac{\partial p}{\partial\rho}\biggr{)}_{n}.$$

$$\dot{n}$$

$$\dot{H}$$

$$\dot{T}$$

$$\frac{dT}{dn}=\frac{T}{n}{\left(\frac{\partial p}{\partial \rho }\right)}_n.$$

$${\rho }_d=Cn=3H^2-\rho (n,T).$$

$$J(n,H,T)=\frac{3H^2-\rho (n,T)}{n}=C=\text{mbox}const.$$

$$\displaystyle p=\frac{\mathcal{N}T}{V}\Bigl{[}1+\frac{\mathcal{N}}{V}F(T)+\Bigl{(}\frac{\mathcal{N}}{V}\Bigr{)}^{2}G(T)+\cdots\Bigr{]}.$$

$$p(n,T)=nT[1+nF(T)].$$

$$\displaystyle\biggr{[}\frac{\partial}{\partial n}\biggl{(}\frac{\rho}{n}\biggr{)}\biggr{]}_{T}=-T^{2}F^{\prime}(T).$$

$$\rho =n[\varphi (T)-n{T}^2{F}^{\prime }(T)],$$

$$\rho =n(m_0+\frac{3}{2}T)-n^2{T}^2{F}^{\prime }(T).$$

$$\dot{n}$$

Full text

Keywords: Dynamical systems, complete integrability, FRWL Cosmology, accelerated expansion, quintessence, phantom fields, real gas, cyclic universe, inflation.

1 Introduction

Supernova observations since 1998 have been showing that, over the last 5 billion years, the cosmic expansion has been accelerating [1, 2]. That is, for a Friedmann-Robertson-Walker-Lemaître (FRWL) Universe, both the scale factor $a(t)$ and its time derivative $\dot{a}(t)$ have been increasing. This has been happening at the same rate, so that Hubble parameter $H=\dot{a}/a$ has been constant. In other words, in the standard cosmological model, the Universe is evolving towards exponential expansion $a(t)\sim\exp(Ht).$ The current epoch has been labeled as dark energy dominated era after un unknown component called dark energy which makes most of the energy content of the Universe and provides the negative pressure, needed to overpower the gravitational pull and explain the resulting accelerated expansion. For a review and list of references, see [3]. Some doubts have been recently cast on whether this is, indeed, the case [4].

The cosmological principle — on a very large scale, the distribution of matter in the Universe is homogeneous and isotropic — most often leads to choosing perfect fluids to model the matter and energy in the universe. These are defined as fluids that are isotropic in their rest frame. Associated with a perfect fluid, there is its equation of state which is a relationship between the pressure $p$ of the fluid and the energy density $\rho$. The vast majority of cosmological models are based on relationship of the type $p=\omega\rho$, where $\omega$ is a constant, independent on time. The Friedmann equation $\ddot{a}/{a}=-(4\pi G/3)(\rho+3p)$ shows that, for cosmic acceleration, it is required that $\omega<-1/3$. On the other hand, the continuity equation (energy conservation equation) for the perfect fluid, $\dot{\rho}+3H(\rho+p)=0$ shows that $\rho+p$ must not be negative for a meaningful cosmological model in which the energy density of an expanding Universe decreases with time. This leads to $\omega\geq-1$. Dark energy could be defined as any physical field, which operates in the gap $-1\leq\omega<-1/3$ between the weak energy condition, $\rho\geq 0$ and $\rho+p\geq 0$, i.e. having positive energy density (to account for the necessary density to make the universe flat) and realistic cosmology, and the strong energy condition, $\rho+p\geq 0$ and $\rho+3p\geq 0$, for which the part $\rho+3p\geq 0$ is violated to account for the needed negative pressure. There are many forms in which the dark energy is sought. For an extensive review and list of references, see [5]. Observationally favoured is the cosmological $\Lambda$ cold dark matter model (corresponding to $\omega=-1$), in which the role of dark energy is played by the cosmological constant $\Lambda$. Quintessence models are based on a dynamical, evolving, spatially-inhomogeneous component with negative pressure [6, 7] and commonly considered quintessential cosmological models introduce a spatially-inhomogeneous slowly-evolving real scalar field rolling down a potential similar to the inflaton field in inflation theory. For a review and list of references, see [8]. Phantom cosmological models [5] violate all four energy conditions [9], $\omega\geq-1$, in particular. From quantum point of view, the phantom models are unstable, but it is not necessarily so from a classical perspective [9]. Such models exhibit Big Rip singularities (the scale factor $a(t)$ becoming singular over a finite time) and there are many proposed remedies for this — see [9] and the references therein. Caldwell introduced [10] the concept of phantom fields by constructing a toy model of a ”phantom” energy component with $\omega<-1$ and argued that it agrees, based on current data and understanding, with most classical tests of cosmology. Carroll et al. claim [11] that that temporary violation of $\omega\geq-1$ is not incompatible with a well-defined model. According to recent studies, see [12] and the references therein, phantom cosmological models are preferred over quintessence ones.

The quintessential scheme can also be achieved without the usually discussed scalar fields. Alternative quintessence models introduce real gas equations of state [13, 14, 15, 16, 17, 18]. General Relativity with a perfect fluid of any type can be recast equivalently as a modified gravity theory [19, 20]. See [21, 22] for a general review. Attempting to reconcile General Relativity and the observed cosmic acceleration, even more exotic dark matter sources are sought. The dark energy equation of state has been generalised in many ways: allowing it to change its structure/form during the universe evolution [23]; introducing dark fluid with a time-dependent equation of state leading to multiple de Sitter space [24]; considering non-linear inhomogeneous equations of state [25]; modifying the equation of state with an arbitrary function of the Hubble parameter $H$ and its derivatives [26]; abandoning the perfect fluid hypothesis [27] and studying the Redlich–Kwong, the modified Berthelot, the Dieterici, and the Peng–Robinson real gasses resulting in $a(t)$ not diverging in any finite time so that any Big Rip is avoided even if $\omega$ may lie today in the phantom regime and others.

The dark energy also provides for consistent models of cyclic cosmology, allowing the Universe to undergo infinitely many self-sustained cycles without failing the second law of thermodynamics, according to which the entropy can only increase, thus necessitating larger and larger successive cycles. The Steinhardt–Turok model [28, 29, 30, 31], based on the ekpyrotic scenario and M-theory, demonstrates a Big Bang–Big Crunch sequence with entropy removal in each cycle. The Baum–Frampton model [32, 33] is based on phantom cosmology. There are many other studies on cyclic Universes — see, for example, [17, 18, 26, 34].

In this study, a quintessence cosmological model, based on the most general real gas is presented as a nonlinear dynamical system of three variables — the number density, the Hubble parameter and the temperature. The nonlinear dynamics is simplified by the existence of two global first integrals. In addition, there are special (second) integrals, defined and conserved on hyper-surfaces in the three-dimensional phase space. These surfaces are invariant manifolds, which separate the phase space into subspaces where certain types of dynamic behaviour takes place.

Following a relatively general model description of real gas cosmology and introducing the three dynamical equations, a real gas with virial equation of state and a van der Waals gas are presented. For the latter, the system is completely integrable and the solutions are found. Various physically relevant possibilities for the critical points are identified and the corresponding dynamic behaviour is studied. Both periodic orbits and unbounded solutions are found in the model, depending on the initial conditions (and the level sets presented by the first integrals). The existing stable periodic solutions (closed orbits), which are models of a cyclic Universe, stem from the presence of the conserved quantities. Interestingly, periodic solutions exhibit regions characterised by inflation and deflation, while unbounded trajectories are characterised by inflation in a ”fly-by” near an unstable critical point.

2 The Set-up

In this model, the Universe is presented classically as a two-component mixture of dust, with energy density $\rho_{d}$ and pressure $p_{d}=0$, and a most general real gas with equation of state expressing the real gas pressure $p$ in terms of the real gas energy density $\rho$ and temperature $T$, i.e. $p=p(\rho,T)$. The energy density of the dust component, $\rho_{d}$, can be taken as positive (for example, one could think of the dust component as of ordinary baryonic matter in this case), zero (absence of dust component), or, to reveal more mathematical aspects of the model, negative. Dust with negative energy density is not a new feature — see [35, 36, 37, 38, 39, 40] and the references therein.

The setting for the analysis of the two-fraction Universe is the Friedmann-Robertson-Walker-Lemaître (FRWL) cosmology with metric:

[TABLE]

where $a(t)$ is the scale factor of the Universe and $k$ is the spatial curvature parameter.

Einstein’s equations are:

[TABLE]

where $\kappa=8\pi G/c^{4}$ and the matter energy-momentum tensor $T_{\mu\nu}$, representing the two fractions of the Universe, collectively modelled with a perfect fluid, is given by:

[TABLE]

Here $\tilde{\rho}=\rho_{d}+\rho$ and $\tilde{p}=p$ are, respectively, the cumulative density and pressure for both fractions and $u^{\mu}=dx^{\mu}/d\tau$ (with $\tau$ being the proper time) is the flow vector satisfying $g_{\mu\nu}u^{\mu}u^{\nu}=-1$.

Friedmann equations are [41]:

[TABLE]

Only the case of flat spatial three-sections ($k=0$) and without cosmological constant (i.e. $\Lambda=0$) will be of interest. Also, Planck units will be used, thus $\kappa=1,c=1,k_{B}=1$.

The Friedmann equations (4)–(5) then become:

[TABLE]

where $H(t)=\dot{a}(t)/a(t)$ is the Hubble parameter. It will be one of the three dynamical variables of the model. As $\ddot{a}/a=\dot{H}+H^{2}$, combining the Friedmann equations allows to express:

[TABLE]

which will be considered as one of three dynamical equations.

In absence of unbalanced particle creation or particle annihilation processes, the number of particles in the perfect fluid is conserved and this is manifested, in locally flat coordinates, by the continuity equation:

[TABLE]

where $n$ is the particle number density and the components of the velocity $\vec{v}$ are the spatial components of the four-velocity:

[TABLE]

(in locally flat coordinates).

Introducing the particle current $n^{\mu}=n\,u^{\mu}$, the particle conservation equation in covariant form can be written as $\nabla_{\mu}n^{\mu}=0$ or

[TABLE]

that is, the co-moving number of particles, $na^{3}$, is constant. The continuity equation for the real gas is:

[TABLE]

and that of the dust is:

[TABLE]

These continuity equations are energy conservation equations.

The final equation is the equation of state. The usual cosmological models are based on a barotropic equation of state $p=\omega\rho$, where $\omega$ is a parameter which is a constant. If $\omega$ is in the range $-1<\omega<-1/3$, the model is called quintessence. Negative pressure is needed to achieve this. The cosmological constant or vacuum energy is modelled by $\omega=-1$. Models with $\omega<-1$ are also characterised by negative pressure and are called phantom field models. For the observed cosmic acceleration, it is required that $\omega$ must be smaller than $-1/3$ — visible from the cosmic acceleration equation, (6). This amounts to $\rho+3p<0$ — a part-violation of the strong energy condition ($\rho+p\geq 0$ and $\rho+3p\geq 0$). Dark energy could be defined as any physical field with positive energy density (to account for the necessary density to make the universe flat) and negative pressure, violating the part $\rho+3p\geq 0$ of the strong energy condition. Phantom cosmological models, on the other hand, violate all four energy conditions, in particular, the weak energy condition: $\rho\geq 0$ and $\rho+p\geq 0$. The phantom field is unstable from a quantum field theory perspective, but not necessarily so from the perspective of classical cosmology.

It the part $\rho+p\geq 0$ of the weak energy condition is violated, then, according to equation (12), an expanding Universe ($H=\dot{a}/a>0$) would have growing density ($\dot{\rho}>0$) and vice versa: a contracting universe will be characterised by diminishing density. As it will be shown, when $\rho_{d}$ is taken as negative in the presented model, there will be initial conditions which lead to quintessence closed trajectories (but not to phantom ones, namely, the projections of the closed orbits in three-dimensional $\rho$–$H$–$T$ phase space onto the $\rho$–$H$ plane will be convex rather than concave, which are associated with phantom cosmologies — it is exactly over the ”indented” region where the weak energy condition is temporarily violated).

Not all of the above equations are independent. Differentiation of (7) with respect to time and substitution into it of: $\dot{H}$ from (6), $\dot{\rho}$ from (12), and $\dot{\rho_{d}}$ from (13), leads to an identity.

Using equation (7), the dust energy density can be expressed as $\rho_{d}=3H^{2}-\rho$ and eliminated from (8) so that the following dynamical equation holds [17, 18]:

[TABLE]

To find the dynamical equation for the temperature of the Universe, consider the Gibbs equation [42, 43]:

[TABLE]

Here $S$ is the specific entropy (entropy per particle) and, as it is a full differential, the following integrability condition must hold in thermodynamical variables $\rho$ and $n$:

[TABLE]

Thus the integrability condition becomes

[TABLE]

For any simple thermodynamical system, subject to the action of the generalised force $Z$, associated with the external parameter $\zeta$, the second initial proposition of thermodynamics leads to the existence of a thermic equation of state: $Z=Z(\zeta,T)$. Then, the following identity is valid: $(\partial Z/\partial\zeta)_{T}\,\,(\partial\zeta/\partial T)_{Z}\,\,(\partial T/\partial Z)_{\zeta}=-1$. With the help of this, the integrability condition results in the following thermodynamic identity:

[TABLE]

If the equation of state is substituted into this identity, a functional relationship between $T,\,\rho$ and $n$ stems.

Using the number conservation equation (11) to express $n$ as $-\dot{n}/(3H)$ and the energy conservation equation (12) to exptess $\rho+p$ as $-\dot{\rho}/(3H)$, the following temperature evolution law is valid [42, 43]:

[TABLE]

namely:

[TABLE]

This will be the second dynamical equation of the system (the temperature $T$ being another dynamical variable).

As third dynamical equation, the continuity equation (11) for the number density will be taken with the number density $n$ as the third dynamical variable. Therefore, the three-dimensional dynamical system is described by the equations:

[TABLE]

In (22), equation of state in the form $p=p(n,T)$ should be used as well as in (23) for the determination of $(\partial p/\partial T)_{n}$. If, however, the equation of state is given as $p=p(\rho,T)$, then the mass density $\rho$ should be excluded with the help of the identity (18) (the integrability condition) by expressing it as $\rho=\rho(n,T)$. From $\rho=\rho(n,T)$, one also calculates the derivative $(\partial\rho/\partial T)_{n}$ in the denominator in (23).

3 Conserved Quantities and Global Behaviour. Virial Gas Example

There is a symmetry in the model: if one divides (23) by (21), an expression independent of $H$ stems:

[TABLE]

Such, $H$-independent, first-order ordinary differential equation exists for the most general gas and leads to the existence of a first integral in the form $I=I(n,T)$, independent of $H$.

In view of their respective continuity equations, (11) and (13), the dust density $\rho_{d}$ and the number density $n$ are proportional: $\rho_{d}=Cn$, where $C$ is a constant having the same sign as $\rho_{d}$ (to ensure positive number density $n$).

In the absence of the dust component ($\rho_{d}=0$), the trajectory in the phase space always lies on the hyper-surface $\rho=3H^{2}$, given by equation (7) with $\rho_{d}=0$. Otherwise:

[TABLE]

This constant $C$ is of paramount importance in the analysis — it is the value of another first integrals present:

[TABLE]

(The fact that $J$ is a first integral can be easily seen by differentiating $J$ with respect to time and substituting the dynamical equations — an identity will follow.)

Absence of dust ($\rho_{d}=0$) corresponds to $C=0$.

If the first integral $I(n,T)=$ const is also known globally, then the presence of two first integrals means that the system is completely integrable.

To illustrate all this with an example, consider real gas whose pressure $p$ is related to the particle number $\mathcal{N}$, the temperature $T$, and the volume $V$ via the virial expansion:

[TABLE]

The term $F(T)$ represents the first correction to the ideal gas equation of state ($p=\mathcal{N}T/V$). Interactions involving three (and more) particles are described by the second correction $G(T)$ (and the following terms) and will not be considered. The temperature $T_{B}$ at which $F(T)$ vanishes (then the real virial gas resembles ideal gas mostly) is called Boyle temperature. For temperatures below $T_{B}$, the correction term $F(T)$ is negative**The virial expansion for a gas, based on Lennard-Jones potential, breaks down not only at high densities, but also at low temperatures: it is divergent as $T\to 0$. This reflects the fact that the attractive interactions lower the energy of the gas constituents and at low temperatures there is a condensation into a liquid phase. (For the van der Waals gas, at temperature $T=0$, the pressure is finite: $p=-Bn^{2}$). When a liquid is under tension, it pulls the confining surfaces. In the van der Waals case, for example, the liquid can withstand a maximum tension given by 27 times its critical pressure. Then the process of cavitation starts — formation of gas bubbles and two phases (liquid and gas) co-exist. This is characteristic to systems with Lennard–Jones potential [44]. The nucleation process can be homogeneous (for smooth and pure liquids) or not (for confined liquids or for liquids with impurities) and to trigger cavitation, the pressure in the fluid must fall below the saturated vapour pressure. Gagnon and Lesgourgues [45] consider $p_{DE}=\omega\rho_{DE}$ as the negative pressure that drives the cavitation process, where, in dark-energy domination, $-1<\omega<-1/3$ and the present-day dark energy density is $\rho_{DE}\sim 10^{-12}$ eV4*. In this case, the saturated vapour pressure is zero and, therefore, cavitation starts when the pressure becomes negative (which is the case if a dark energy component is present). The tensile strength of a fluid characterises its ability to counter-act the cavitation process: bubbles do not grow for as long as the force due to the intra-molecular interactions balances the outward pressure [45]. The tensile strength of a fluid is equal to the absolute value of the minimum negative pressure which the fluid can sustain without breaking apart and to achieve dark energy domination, while still maintaining a valid fluid description (the more bubbles form, the more the hydrodynamic description fails), the pressure should be negative but its absolute value should not grow above the tensile strength of the fluid [45]. For this reason, real virial gas description would not be applicable below a certain (very low) temperature, or the term $F(T)$ will have to be regularised..

For example, when $F(T)=A-B/T$, where $A$ and $B$ are two positive constants (with $A$ measuring the strength of the attractive force between the gas ingredients and $B$ being the volume per mole of substance), the van der Waals equation of state results.

In terms of the number density $n$, the virial equation of state is:

[TABLE]

Assuming $n$ and $T$ as the independent thermodynamical variables, finding $(\partial p/\partial T)_{n}$ from the above and substituting it, together with $p$, into (18), results in the following differential equation:

[TABLE]

This integrates directly into:

[TABLE]

where $\phi(T)$ can be determined as follows. One can take an ideal gas limit, i.e. set the co-efficient $F(T)$ of the second term of the virial expansion to zero. If this gas is monoatomic and has three translational degrees of freedom, the average kinetic energy of the particles is $(3/2)T$. Also, $n=(\mathcal{N}m)/(Vm)=(M/V)(1/m)=\rho/m$, where $M$ is the mass of the system and $m$ is the relativistic mass of a representative particle: $m=m_{0}+(1/2)m_{0}u^{2}+O(u^{4}).$ Here $m_{0}$ is the rest mass and $u$ — the speed of the particle. Thus, the mass density of an ideal gas can be approximately written as $\rho=n[m_{0}+(3/2)T]$. Therefore, $\phi(T)=m_{0}+(3/2)T$. The number density $n$, the mass density $\rho$ and the temperature $T$ of a real virial gas are related as follows:

[TABLE]

Thus the dynamical equations for the virial gas become:

[TABLE]

One also has to make sure that the mass density $\rho$ is not negative. This sets a range of allowed values of the number density $n$ as a function of the temperature $T$. In the case of real virial gas, requesting $\rho\geq 0$ in equation (31), leads to $n\leq[m_{0}+(3/2)T]/[T^{2}F^{\prime}(T)]$ for positive $F^{\prime}(T)$. If $F^{\prime}(T)$ is negative, then $\rho>0$ is always satisfied. In other words, one either has to have: (i) a ”standard” virial gas, for example, a van der Waals gas or gas based on Lennard–Jones type of potential, for which $F(T)=\alpha-\beta e^{\epsilon/T}$ and $F^{\prime}(T)>0$ always ($\alpha,\,\beta,$ and $\epsilon$ are all positive constants), combined with an upper limit on $n$ (which conforms with the virial expansion over the powers of small $n$), or (ii) an anomalous virial gas for which $F(T)$ is a function monotonically decreasing with the temperature ($n$ does not have to have an upper limit in this case).

For the real virial gas one has:

[TABLE]

As it has been made sure that $n$ is positive, one can keep track of when $\rho+p\geq 0$ (part of the weak energy condition) — it will be satisfied for as long as:

[TABLE]

For example, for the van der Waals gas, the mass density will be positive for as long as $n\leq[m_{0}+(3/2)T]/B$. For a gas of electrons ($m_{0}=511$ keV or $10^{9}$ K approximately), for quite high temperatures $T$, the term $T/m_{0}$ will be quite small and negligible and the upper limit on the number density could be considered as $m_{0}/B$. This is well above what is assumed for $n$ as the virial expansion is over the powers of small $n$. On the other hand, for a van der Waals gas, $\rho+p\geq 0$ will be satisfied for temperatures $T$ such that $T\geq(2Bn-m_{0})/(An+5/2)$. If the density $n$ is below $m_{0}/(2B)$, which is quite big, then $\rho+p\geq 0$ will always be satisfied.

The first integral $J(n,H,T)$ for the case of a virial gas is:

[TABLE]

It follows from (25) that $n=0$ leads to $H=0$. The plane $n=0$ is tangent to the hyper-surface $J(n,H,T)=C$. This hyper-surface will be entirely above the plane $T=-(2/3)(C+m_{0})$, as point $[n=0,H=0,T=-(2/3)(C+m_{0})]$ is its absolute minimum. For $C=-m_{0}$, the plane $T=0$ is tangent to the hyper-surface. For $C>-m_{0}$, the plane $T=0$ crosses the hyper-surface. Thus, the origin necessarily lies on the hyper-surface for any $C\geq-m_{0}$. As already mentioned, the signs of $\rho_{d}$ and $C$ are the same, while absence of dust ($\rho_{d}=0$) corresponds to $C=0$. As this case serves as a natural separator, it will be useful to introduce the separatrix hyper-surface $J(n,H,T)=0$ (see Figure 1a). The absolute minimum of this separatrix is at $T=-2m_{0}/3$. Thus, if a hyper-surface is above this separatrix, it will correspond to a negative $\rho_{d}$, while if it is below — to a positive $\rho_{d}$.

4 The Completely Integrable Van der Waals System

In the case of a van der Waals gas, the first integral $J$ is simply (see Figure 1b):

[TABLE]

For a van der Waals gas, the first integral $I(n,T)$ can be obtained explicitly in a straight-forward manner: using $F(T)=A-B/T$ and dividing (23) by (21) yields the ordinary differential equation in separate variables:

[TABLE]

The solution of (39) is:

[TABLE]

where $D$ is a positive constant, equal to the value of the other first integral for the van der Waals gas (see Figure 1b):

[TABLE]

Therefore, the van der Waals system is completely integrable.

Obviously, $n=0$ leads to $T=0$ for all $H$, including $H=0$. Thus, the hyper-surface $I(n,T)=D$ contains the origin for any value of the constant $D$.

The trajectories in the van der Waals phase space are obtained as intersections of the two hyper-surfaces given by the two first integrals $J_{vdw}(n,H,T)=3H^{2}/n-m_{0}-(3/2)T+Bn=C\mbox{ = const}$ and $I_{vdw}(n,T)=T\,n^{-2/3}e^{-2An/3}=D\mbox{ = const }>0$. All possible cases are given on Figure 2 (see also Figures 3 and 4).

From the first integrals, one can immediately find the solutions:

[TABLE]

where $n(t)$ is determined by separation of variables from $\dot{n}=-3Hn$:

[TABLE]

For small $n$, one finds easily that

[TABLE]

— behaviour similar to that in the $T=0$ plane. The origin is reachable in infinite time.

Formula (45) is valid for $C>-m_{0}$. When $C<-m_{0}$, the trajectory is either an open curve (Figures 2c and 2d) or a periodic one (Figure 2d), but, in either case, never over regions where $n$ approaches zero asymptotically, thus $C<-m_{0}$ is not applicable for small $n$. In the limiting (separatrix) case $C=-m_{0}$ and for small $n$, the relevant situations are the ones given by the open trajectories on Figure 2f and by the closed curves, terminating at the origin, on Figure 2g. Then one has $n(t)\sim(t-t_{0})^{-6/5}$ and, qualitatively, the same asymptotic behaviour — with a negative power of $t$ and trajectory extinguishing in the origin in infinite time.

As it is usually done, one can extend the validity of the model by allowing the consideration of large number densities. For the unbounded van der Waals trajectories and for large values of $n$ (see Figure 2 and Figure 3b), a blow-up is observed in finite time. In order to find the asymptotic behaviour near the time of the blow-up, retain the leading terms in (44). This yields the equation

[TABLE]

where $\sigma=\mathrm{sign}(H)$ and $t^{*}$ is an integration constant.

The integral on the left-hand side has, for $n\to\infty$, the asymptotic behaviour $-\frac{3}{A}n^{-11/6}\exp(-\frac{An}{3}),$ thus

[TABLE]

When $n\to\infty$, the left-hand side approaches zero and hence $t\to t^{*}.$ Therefore, $t^{*}$ is the blow-up time and the above formulae are valid for $t<t^{*}$ only. This is clearly possible only when $\sigma=-1$ and $H<0$. Hence $H\to-\infty$ and this blow-up represents a Big Crunch:

[TABLE]

The integration in (44) cannot be performed explicitly in the general case, but the qualitative behavior can be inferred from the previous (and the following) analysis.

5 Critical Points and Second Integrals

At an equilibrium point of the general dynamical system (21)–(23), all time derivatives on the left-hand sides of the equations vanish simultaneously.

Clearly, if $H=H^{*}=0$ and if the pressure vanishes, i.e. $p(n^{*},T^{*})=0$, then the set $\{[n^{*},H^{*}=0,T^{*}(n^{*})]\}$ is a continuum of equilibrium points. Separately, if $p(n,0)=0$, then there is another continuum of equilibrium points: the set $\{[n^{**},H^{**}=0,T^{**}=0]\}$. Finally, if $p(0,T)=0$, then there is a third continuum of equilibrium points: the set $\{[n^{***}=0,H^{***}=0,T^{***}]\}$.

As mentioned earlier, the system will choose one equilibrium point from the curve $T^{*}(n^{*})$. That is, there is another equation which, together with $p(n^{*},T^{*})=0$, determines $n^{*}$ and $T^{*}$. This can be seen as follows. $I(n,T)$ is constant, that is, for the initial condition $(n_{0},T_{0})$ at initial time $t=t_{0},\,\,$ $I(n,T)=I(n_{0},T_{0})=I_{1}=$ const. On the other hand, the curve $T^{*}=T^{*}(n^{*})$ with the equilibrium points (Figure 5) intersects $I(n,T)$ exactly at point with coordinates $(n^{\ast},T^{\ast})$, namely, the equation $I(n^{\ast},T^{\ast})=I(n_{0},T_{0})=I_{1}=$ const, together with the equation $p(n^{\ast},T^{*})=0$ are the two simultaneous equations giving the pair $n^{\ast}$ and $T^{\ast}$ — as dependent on the choice of initial conditions.

Not all trajectories in the phase space are affected by the equilibrium points $(n^{\ast},H^{*}=0,T^{\ast})$. This is due to the presence of second integrals in the phase space which fragment the phase space by non-crossable ”walls” with each piece being a centre manifold. A second integral $K(\vec{x})=0$ of an autonomous dynamical system of the type $\dot{\vec{x}}(t)=\vec{f}[\vec{x}(t)]$ is an invariant, but only on a restricted subset, given by its zero level set [46]. It is defined by $\dot{K}(\vec{x})=\mu(\vec{x})K(\vec{x})$. Second integrals reduce to first integrals when $\mu=0$ and to time-dependent first integrals when $\mu=$ const [46]. The surface $K_{1}$, defined by $3H^{2}-\rho(n,T)=0$, is a second integral and an invariant manifold, i.e. all trajectories originating from this invariant manifold remain there, and no trajectories originating from outside can penetrate this manifold, i.e. no trajectory can cross $3H^{2}-\rho(n,T)=0$. To see that the surface $K_{1}$ is a second integral, differentiate with respect to time: $(d/dt)[3H^{2}-\rho(n,T)]=6H\dot{H}-\dot{\rho}=-3H[3H^{2}-\rho(n,T)]$, as can be seen from the equations of motion. The hyper-surface $3H^{2}-\rho(n,T)=0$ corresponds to absence of dust, i.e. $C=0$ [this is the Friedmann equation (7) with $\rho_{d}=0$, namely, $C=0$]. Therefore, in the absence of dust, there will be a relationship between $n$, $H$, and $T$, stemming from $3H^{2}-\rho(n,T)=0$. In view of this, only two of the three dynamical equations (21)–(23) are independent and the motion in the three-dimensional phase space is along the hyper-surface with equation $3H^{2}-\rho(n,T)=0$.

No trajectory can cross the parabolic wall given by $K_{1}$. This can be seen from a physical point of view too — during the evolution, the dust component cannot disappear (or change the sign of $\rho_{d}$).

The surface $K_{2}$, defined by $n=0$, is another invariant manifold, i.e. $n=0$ is another second integral. For example, for the virial gas, the motion on the $K_{2}$ surface is given by $\dot{H}=-(3/2)H^{2}$ and $\dot{T}=-2HT$. These equations integrate easily and the solutions for initial data $(H_{0},T_{0})$ at $t=t_{0}$ are

[TABLE]

Thus, the long-time asymptotic is $(H,T)\to(0,0).$ The trajectories in $n=0$ converge to the origin over infinite time (bearing in mind that the origin is a critical point too). As the $T$ axis is a continuum of equilibrium points, should the equation of state satisfies $p(n,0)=0$, then the critical points $(0,0,T^{***})$ are not reachable due to the second integral $K_{2}$. Thus, a trajectory cannot end anywhere on the invariant plane $n=0$, except at the origin.

There is a third invariant manifold — the second integral $K_{3}$, defined as the surface $T=0$. The motion on this surface is given by the equations $\dot{H}=-(3/2)H^{2}$ (which is the same equation as the equation for $\dot{H}$ on the surface $K_{2}$) and $\dot{n}=-3Hn$. For initial data $(n_{0},H_{0})$ at $t=t_{0}$, the equations again integrate easily and the solutions are

[TABLE]

The long-time asymptotic is again towards the origin, $(n,H)\to(0,0),$ where the trajectories converge over infinite time. Similarly, if the equation of state satisfies $p(0,T)=0$, then the critical points $(n^{**},0,0)$ are not reachable due to the second integral $K_{3}$. Thus, a trajectory cannot end anywhere on the invariant plane $T=0$, except at the origin.

The equilibrium points of interest are the origin and the continuum set $\{[n^{*},H^{*}=0,T^{*}(n^{*})]\}$ where the pair $n^{*}$ and $T^{*}$ are determined by the simultaneous equations $p(n^{*},T^{*})=0$ and $I(n^{*},T^{*})=I(n_{0},T_{0})=I_{1}=$ const.

6 Eigenvalues, Eigenvectors and Trajectories. Linearization

Consider next the linearised form of the dynamical system near the equilibrium point $[n^{\ast},\,\,H^{\ast}=0,\,\,T^{\ast}(n^{*})]$:

[TABLE]

where the stars on the derivatives indicate that they are taken at the equilibrium point $(n^{\ast},0,T^{\ast})$.

In matrix form this can be written as:

[TABLE]

where:

[TABLE]

The stability matrix $L(n,H,T)$ at the equilibrium point $[n^{\ast},\,\,H^{\ast}=0,\,\,T^{\ast}]$ is:

[TABLE]

One of the eigenvalues of $L^{\ast}$ is zero ($\lambda_{1}=0$), while the other two satisfy

[TABLE]

If the expression on the right-hand side is negative and, say, equal to $-\omega^{2}$, then $\lambda_{2,3}=\pm i\omega$. If positive and equal to $q^{2}$, then $\lambda_{2,3}=\pm q$. If zero, then all three eigenvalues are zero††† A critical (fixed) point $P$ of a dynamical system is called stable if for any neighbourhood $U(P)$ there is a neighbourhood $V(P)$, such that every trajectory that starts in $V(P)$ does not leave $U(P)$ [47]. A periodic trajectory $\gamma(t)$ is called orbitally stable if for any neighbourhood $U(\gamma)$ there is a neighbourhood $V[\gamma(0)]$, such that every trajectory that starts in $V[\gamma(0)]$ does not leave $U(\gamma)$ [47].

If the two eigenvalues are imaginary ($\lambda_{2,3}=\pm i\omega$), given the existence of a globally defined first integral $J(n,H,T)$, and if the hyper-surface, defined by the other first integral $I(n,T)$, is well behaved around $(n^{*},0,T^{*})$ (which is to be expected from the physical context), then there are periodic trajectories which are always confined by the integral surfaces and are orbitally stable. They are orbiting around the critical point $(n^{*},0,T^{*})$ and this critical point is stable. If the eigenvalues are real, $\lambda_{2,3}=\pm q$ (they are with opposite signs), then the critical point $(n^{*},0,T^{*})$ is unstable.

The eigenvectors $u_{i}$, corresponding to eigenvalues $\lambda_{i},\,\,i=1,2,3,\,$ are:

[TABLE]

where the choice $\pm i\omega$ or $\pm q$ depends on the sign of the right-hand side of (74).

It is interesting to study the closed trajectories which correspond to $\lambda_{2,3}=\pm i\omega$ (see Figure 6).

Replacing $L(n^{\ast},H^{\ast},T^{\ast})$ in (56) with $M\Lambda M^{-1}$, where $M$ is the matrix whose columns are the eigenvectors $u_{1,2,3}$ and $\Lambda$ is the diagonal matrix with the eigenvalues $\lambda_{1,2,3}$ along its main diagonal, multiplying (56) from the left with $M^{-1}$, noting that $M^{-1}=$ const, and introducing the column-vector $Z(t)=[z_{1}(t),z_{2}(t),z_{3}(t)]=M^{-1}X(t)$, allows to write the dynamical system in diagonalised form: $\dot{z}_{i}(t)=\lambda_{i}z_{i}(t)$ with $i=1,2,3.$ Thus $z_{1}=\alpha_{1}=$ const, $z_{2,3}(t)=\alpha_{2,3}e^{\pm i\omega t},$ where $\alpha_{2,3}=$ const.

Using $X(t)=MZ(t)$ and requesting reality for the dynamical variables — which, in turn, necessitates that $z_{2}(t)=\bar{z}_{3}(t)$, that is, $\alpha_{2}=\bar{\alpha}_{3}\equiv\alpha$ — allows to find:

[TABLE]

One can immediately deduce from here that:

[TABLE]

The right-hand side is a constant, which is not necessarily small. In the linear approximation however, $n$ is very close to $n^{*}$ and $T$ is very close to $T^{*}$. Thus, this constant must be set equal to zero, namely, either $\alpha_{1}$ should be zero or the bracketed term on the right-hand side should be zero. The latter cannot be taken as zero because this would impose an extra relationship on $p$, $n$, $\rho$, and $T$. Thus one must have $\alpha_{1}$ equal to zero.

Without loss of generality, $\alpha$ can be taken as positive.

The trajectory in the three-dimensional phase space is a closed curve and the equation of its projection (in the linear approximation) onto the $n$–$H$ plane is an ellipse (see Figure 6c):

[TABLE]

while the equation of its projection onto the $T$–$H$ plane is another ellipse (Figure 6d):

[TABLE]

These are two ellipses with semi-axes proportional to $\alpha$ in each of the cases. The constant $\alpha$, itself, is fixed by the initial conditions.

Equation (88) is the equation of a plane and this plane is the tangent plane at point $(\rho^{\ast},0,T^{\ast})$ (the centre) to the surface along which the dynamical system evolves in the three-dimensional phase space. This plane (88), intersected with the plane $H=0$, yields the straight line (Figure 6b):

[TABLE]

in the linear approximation.

To illustrate the above linear approximation with an example, start with the eigenvalues for a real virial gas system with equation of state given by (28). These are $\lambda_{1}=0$ and $\lambda_{2,3}$ satisfying:

[TABLE]

If one again considers a van der Waals gas [for which $F(T)=A-B/T$], then:

[TABLE]

These are purely imaginary if $n^{*}<(\sqrt{3/2}-1)/A$. Given that $n^{*}=T^{*}/(B-AT^{*})$ for the van der Waals gas, one gets that closed trajectories exist for $T^{*}<(1-\sqrt{3/2})T_{B}\approx 0.1835\,T_{B}$ (recall that $T_{B}=B/A$). If this is the case, then the trajectory in the linear approximation is an ellipse and the projection of this ellipse onto the the $n$–$H$ plane is the same as (89), while its projection onto the $T$–$H$ plane is

[TABLE]

From these two, or from (91), one gets that the projection of the trajectory onto the $n$–$T$ plane is given by

[TABLE]

7 Periods of Inflation and Deflation for the Linearized Theory

To achieve solution describing inflation [by definition, inflation is equivalent to $\ddot{a}(t)>0$ and $\dot{a}(t)>0$], it is obvious from (14) that negative pressure (28) is needed. Inflation occurs when

[TABLE]

since the scale factor $a(t)$ is always strictly positive. That is, when $p<-H^{2}$.

The inflation parameter $\varepsilon_{\mbox{\tiny I}}$, given by:

[TABLE]

must be less than 1. Here $dN(t)=d\ln a(t)=H(t)dt$ measures the number $N$ of $e$-folds of inflationary expansion.

For the most general gas which admits negative pressure, let $G(n,T)$ denote the region(s) in the $n$–$T$ plane for which $p(n,T)<-H^{2}$. If inflation occurs, it will be over $G(n,T)$.

To analyse the situation in the linear approximation (Figure 6), it is most convenient to firstly find the scale factor $a(t)$. When $H(t)$ is given by (86), i.e. $H(t)=2\alpha\omega\sin\omega t$ (purely imaginary eigenvalues):

[TABLE]

Recalling that $\alpha>0$, positivity of $\dot{a}(t)$ means that $0<\omega t<\pi$, while positivity of $\ddot{a}(t)$ means:

[TABLE]

Note that $1/(4\alpha)+[1+(4\alpha)^{-2}]^{1/2}$ is greater than 1 for all $\alpha$, thus the second inequality is always satisfied.

On the other hand, $-1<1/(4\alpha)-[1+(4\alpha)^{-2}]^{1/2}<0$ for all $\alpha$. It is convenient to introduce $t_{\alpha}$ via

[TABLE]

so that $\pi/2<\omega t_{\alpha}<\pi.$ Therefore, inflation occurs for times $t$ satisfying

[TABLE]

The upper limit depends on the initial conditions via $\alpha$. Therefore, inflation occurs for $\omega t$ between 0 and, at least, $\pi/2$.

One should note that the model also includes deflation — occurring for $\omega t$ between $\pi$ and, at least $3\pi/2$ (this can be seen in analogous manner).

To see how much inflation is generated, consider the number $N$ of $e$-folds of inflationary expansion: $dN(t)=H(t)dt$. This integrates easily to give $N=2\alpha$ and this depends on the initial data through $\alpha$.

If the eigenvalues are real, then the trajectories are open curves and $H(t)=2\alpha q\sinh qt$. Thus:

[TABLE]

Since $\alpha$, $q$, and time $t$ are all positive, then $\dot{a}(t)$ will also be positive. The acceleration $\ddot{a}(t)$ is also positive. Therefore, a “fly-by” near an unstable critical point (i.e. in the neighbourhood of the unstable critical point, defined as the region where the linear approximation is applicable) is always characterised by inflation.

For a van der Waals gas at very high densities and temperatures, the leading term in $\dot{H}$ from (33) is $-\frac{A}{2}n^{2}T$, while the leading term in $H^{2}$ from (38) is $\frac{3}{2}nT$. Thus

[TABLE]

which disagrees with the inflation condition. Therefore, near the blow-up time inflation is not observed.

8 Hamiltonian Formulation

The two integrals of motion completely determine the global behaviour of the full system — see Figure 2. In the case of van der Waals gas, using $J_{vdw}(n,H,T)$ from (38) and with the help of (39), allows to exclude the temperature from the picture: $T(n)=Dn^{2/3}\exp(2An/3)$. The dynamics is then governed by:

[TABLE]

where $p(n)$ is the pressure expressed in terms of the number density $n$ and for the van der Waals gas it is given by:

[TABLE]

The critical points of the underlying two-dimensional system occur where the Hubble parameter and the pressure $p(n)$ vanish, namely, these are the points $(n^{\ast}=0,H^{\ast}=0)$ and $(n^{**},H^{**}=0)$. The zero $n^{**}$ of (110) depends on the van der Waals gas parameters $A$ and $B$ and on the initial conditions via $D$.

For the most general equation of state, solving (24) leads to the appearance of the implicit integral $I(n,T)=$ const and, again, $T=T(n)$ can be implicitly excluded, resulting in the same dynamical system: (108) and (109), with the relevant $p(n)$.

Introducing new variables $u(n)=2/(3\sqrt{n})$ and $v=H/\sqrt{n}$, the system becomes

[TABLE]

where $\varphi(u)=-p(n)/(2\sqrt{n}).$

It can be written in terms of the canonical variables $u$ and $v$ with Hamiltonian (conserved quantity)

[TABLE]

and dynamics:

[TABLE]

Even in the case of a most general equation of state, the underlying two-dimensional system is Hamiltonian. As it is well known, for Hamiltonian systems the only allowed critical points are centres and saddles. Therefore, in the case of imaginary eigenvalues, one is dealing with a centre-type behaviour.

The system corresponds to one-dimensional motion of a particle in a potential field $\mathscr{V}(u)=-\int_{0}^{u}\varphi(\widetilde{u})d\widetilde{u}$ — see [48] and [49]. The centre-type oscillations are near the local minima of $\mathscr{V}(u)$.

It is not difficult to check that the centres of the canonical system transform back to centres of the original system (108), (109) in terms of $n$ and $H$.

9 Discussion

The type of the dynamical behaviour of the system is influenced by two groups of parameters — the initial data of the system $(n_{0},H_{0},T_{0})$ and the gas parameters from the equation of state. The integral curve is an intersection of the level sets of two global first integrals and the possible integral curves are classified in dependence on the aforementioned parameters. A distinct feature of the system are the stable periodic orbits around the critical point $[n^{*},0,T^{*}(n^{*})]$. They occur in the case of phantom dust and for imaginary eigenvalues — see (74) and (92) for the case of virial gas. Such scenario depends on the gas parameters as well, which can be seen from formula (93) for the van der Waals gas. Inflationary and deflationary periods in the cyclic regime are observed (103). Alternatively, when the eigenvalues at $[n^{*},0,T^{*}(n^{*})]$ are real (positive and negative), there are no closed orbits, but all “fly-by” trajectories near $[n^{*},0,T^{*}(n^{*})]$ are characterised by inflation (104)–(106).

In the case of positive dust density, the trajectories end up asymptotically either at $(0,0,0)$ or are unbounded (Big Chrunch). The special explicit solutions in the invariant planes $n=0$ (49)–(50) and $T=0$ (51)–(52), as well as the exact van der Waals solution in this asymptotic regime (45), indicate that the trajectories terminating at the origin $(n,H,T)\to(0,0,0)$ behave asymptotically as negative powers of $t$ and reach the origin in infinite time. This excludes possible inflation near the origin. The unbounded van der Waals trajectories (Figure 3b) represent a blow-up in finite time (Big Crunch) and are not inflatory.