Van der Waals Universe with Adiabatic Matter Creation

TL;DR

This paper explores a cosmological model with van der Waals fluid and dust, incorporating adiabatic particle creation, analyzing its global dynamics and inflationary regimes through a reduced dynamical system.

Contribution

It introduces a novel analysis of a three-component cosmological model with particle creation, reducing the system using a global first integral, and studies specific cases numerically.

Findings

Existence of trajectories with temporary inflationary regimes.

Reduction of the dynamical system to two components using a first integral.

Numerical examples illustrating special cases of particle creation rates.

Abstract

A FRWL cosmological model with perfect fluid comprising of van der Waals gas and dust has been studied in the context of dynamical analysis of a three-component autonomous non-linear dynamical system for the particle number density $n$, the Hubble parameter $H$, and the temperature $T$. Perfect fluid isentropic particle creation at rate proportional to an integer power $\alpha$ of $H$ has been incorporated. The existence of a global first integral allows the determination of the temperature evolution law and hence the reduction of the dynamical system to a two-component one. Special attention is paid to the cases of $\alpha = 2$ and $\alpha = 4$ and these are illustrated with numerical examples. The global dynamics is comprehensively studied for different choices of the values of the physical parameters of the model. Trajectories in the $(n, H)$ phase space are identified for which temporary inflationary regime exists.

Tables

Critical Points	Parameters	$\alpha =2$					$\alpha =4$
$\left(n^{\ast }\ne 0,H^{\ast }=0\right)$		$\beta$						$n^{\ast }>{\tilde{n}}^{\ast }$
					$\beta >\frac{\sqrt{8B}}{m_0}$
				$n^{\ast }>{\tilde{n}}^{\ast }$		$n^{\ast }>{\tilde{n}}^{\ast }$
		Always existing: either unstable spiral in — Fig. 2a, 2b, 4a, 5a, 5d, or unstable node in , — Fig 2a, 2b		Always existing saddle, Fig. 2a, 2b, 4c,5c, 5e, 6e	Unstable node, Fig. 2a, 2c, 6a, 6d	Saddle, Fig. 2a, 2c, 6b, 6e	Always existing unstable spiral with purely imaginary eigenvalues (centre manifold theory) — Fig. 7a, 7c	Always existing saddle — Fig. 7c
	$\tau >\tilde{\tau }$	Does not exist, Fig. 2a, 4c, 5f, 6f					Does not exist
$(0,0)$		Attracts trajectories from the upper half-plane $H>0$
		Repels trajectories from the lower half-plane

Critical Points	Parameters	$T\left(n^{\ast \ast }\right)=T^{\ast \ast }\left(n^{\ast \ast }\right)$ at:	$\alpha =2$ and $\alpha =4$
			$\beta >{\beta }_0$	$\beta ={\beta }_0$
$\left(n^{\ast \ast },H^{\ast \ast }={\beta }^{\frac{1}{\alpha -1}}\right)$	$\tau >\hat{\tau }\left(\beta \right)$	$n^{\ast \ast }={\hat{\nu }}_0^{\ast \ast }$	Saddle, Fig. 3a, 3b	Saddle, Fig. 3a, 3b	Saddle, Fig. 3a, 3b	Saddle, Fig. 3a, 3b
		$n^{\ast \ast }={\hat{\nu }}_0^{\ast \ast }$	Saddle, Fig3b, 3c, 7a, 7b, 7c	Saddle, Fig3b, 3c, 7a, 7b, 7c	Saddle, Fig3b, 3c, 7a, 7b, 7c	Saddle, Fig3b, 3c, 7a, 7b, 7c
		$n^{\ast \ast }={\hat{\nu }}_1^{\ast \ast }$	Stable, , Fig. 3b, 3c, 7d	Stable, , Fig. 3b, 3c, 7d	Stable, , Fig. 3b, 3c, 7d	Saddle for all $n^{\ast \ast }$
		$n^{\ast \ast }={\hat{\nu }}_2^{\ast \ast }$	Saddle, ${\hat{\nu }}_2^{\ast \ast }>{\sigma }_2^{\ast \ast }$, Fig. 3b, 3c	Saddle, ${\hat{\nu }}_2^{\ast \ast }>{\chi }_2^{\ast \ast }$, Fig. 3b, 3c	Saddle, ${\hat{\nu }}_2^{\ast \ast }>{\xi }_2^{\ast \ast }$, Fig. 3b, 3c	Saddle for all $n^{\ast \ast }$

Equations

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

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=d{t}^2-a^2(t)[d{r}^2+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2)],$$

$$T_{\mu \nu }=(\rho +{\rho }_d+p+\Pi )u_{\mu }{u}_{\nu }-(p+\Pi )g_{\mu \nu },$$

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

$$H^2$$

$$\dot{H}=-\frac{3}{2}{H}^2-\frac{1}{2}(p+\Pi ).$$

$$\Gamma =3\beta H^{\alpha },$$

$$\dot{n}=-3nH+n\Gamma =-3nH(1-\beta H^{\alpha -1}).$$

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

$$\dot{{\rho }_d}+3H{\rho }_d$$

$$\Pi =-\frac{\Gamma (\rho +p)}{3H}=-\beta (\rho +p)H^{\alpha -1}.$$

$$\dot{H}=-\frac{3}{2}{H}^2-\frac{1}{2}\left[p(n,T)(1-\beta H^{\alpha -1})-\beta \rho (n,T)H^{\alpha -1}\right].$$

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

$$\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\dot{\rho}(n,T)=\biggl{(}\frac{\partial\rho}{\partial n}\biggr{)}_{T}\dot{n}+\biggl{(}\frac{\partial\rho}{\partial T}\biggr{)}_{n}\,\,\dot{T}.$$

$$\dot{\rho }(n,T)=(\rho +p)(\Gamma -3H)=-3(\rho +p)H(1-\beta H^{\alpha -1}).$$

$$\displaystyle-3(\rho+p)H(1-\beta H^{\alpha-1})=-3nH(1-\beta H^{\alpha-1})\,\biggl{(}\frac{\partial\rho}{\partial n}\biggr{)}_{T}+\biggl{(}\frac{\partial\rho}{\partial T}\biggr{)}_{n}\,\,\dot{T}.$$

$$\displaystyle\dot{T}=-3H(1-\beta H^{\alpha-1})T\biggl{(}\frac{\partial p}{\partial\rho}\biggr{)}_{n}=-3H(1-\beta H^{\alpha-1})T\frac{\left(\frac{\partial p}{\partial T}\right)_{n}}{\Bigl{(}\frac{\partial\rho}{\partial T}\Bigr{)}_{n}}$$

$$p(n,T)=nT(1+An)-B{n}^2,$$

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

$$\rho =n[\varphi (T)-Bn].$$

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

$$\dot{n}$$

$$\dot{H}$$

$$\dot{T}$$

$$\displaystyle\frac{dT}{dn}=\frac{2T(1+An)}{3n}>0\,\mbox{ as $n>0$,}$$

$$T(n)=\tau n^{\frac{2}{3}}{e}^{\frac{2An}{3}},$$

$$\dot{n}$$

$$\dot{H}$$

$$I(n,T)=T{n}^{-\frac{2}{3}}{e}^{-\frac{2An}{3}}=\tau \text{mbox}=const>0.$$

$$L_{11}$$

$$L_{12}$$

Full text

Keywords: Dynamical systems, FRWL cosmology, accelerated expansion, real gas.

1 Introduction

The acceleration of the cosmic expansion and observational data (Supernovæ  Type Ia, Cosmic Microwave Background, Baryon Acoustic Oscillations) are fit best by the current concordance model — the $\Lambda$CDM model which incorporates Dark Energy, modelled by the cosmological constant $\Lambda$, and cold (pressureless) Dark Matter. There are open issues in relation to such model — the so called Cosmological Coincidence Problem: it is known observationally that the present values of the densities of dark energy and dark matter are of the same order of magnitude while, under the $\Lambda$CDM model, the dark-energy density is constant and the dark-matter density is proportional to the inverse third power of the scale factor with the ratio of the two densities varying in time from infinity to zero. There are numerous alternative models, not without open issues on their own, which accommodate acceleration of the cosmic expansion: modified gravity theories, inhomogeneous cosmologies, gravitationally induced particle creation models. In the literature, special attention has been gathered by the adiabatic, or isentropic, production [1, 2, 3, 4] of perfect fluid particles in which the specific entropy (entropy per particle) is conserved (with “isentropic” referring to this). There is overall entropy production due to the enlargement of the phase space of the system as the particle number increases. The imposed condition of conserved specific entropy during the production of perfect fluid particles leads to a simple relationship between the particle production rate and particle “creation” pressure. Zimdahl [5] studies cosmological particle production with production rate which depends quadratically on the Hubble rate $H$ and confirms the existence of solutions which describe a smooth transition from inflationary to non-inflationary behavior. The present work falls in this category and offers a full dynamical analysis of isentropic perfect fluid particle production rate that depends on $H^{\alpha}$ with $\alpha$ being a positive integer. Special attention is paid to the cases of $\alpha=2$ and $\alpha=4$, but the analysis can be easily extended to any other integer positive values of $\alpha$, including odd values — due to the second law of thermodynamics, these work in the regime of expansion only [6]. The setting of the proposed model is a flat FRWL Universe with perfect fluid comprising of two fractions: real gas wit van der Waals equation of state and dust and the tools used are those of dynamical systems, see for example, [7, 8], and as those used in the study of $n$–$H$–$T$ (where $n$ is the particle number density, and $T$ is the temperature) dynamical analysis of cosmological quintessence real gas model with a general equation of state [9]. The dynamical variables are again $n$, $H$, and $T$, but due to the existence of a global first integral (in addition to second integrals), the temperature evolution law has been easily determined and the dynamical system reduced to a two component one over the $(n,H)$ phase space. Inflationary regime with exit from the inflationary behaviour has been identified, both for $\alpha=2$ and for $\alpha=4$, and full classification of the possible phase-space trajectories, subject to the variation of the several physical parameters of the model, has been provided.

2 The Model

This paper studies a Universe modelled classically as a fluid comprising of a binary mixture of dust with energy density $\rho_{d}$ and pressure $p_{d}=0$ and a van der Waals gas with equation of state

[TABLE]

where $p$ is the pressure, $T$ is the temperature, $n=N/V$ — the number of particles $N$ per unit volume $V$ — is the particle number density and $F(T)$ is the term describing two-particle interaction: $F(T)=A-B/T$, where $A$ and $B$ are positive constants***To aid the analysis, a numerical example is presented in this work. It is for van der Waals gas, the parameters of which are $A=1/100$ and $B=10$..

The Universe is described, using Planck units, by the flat Friedmann–Robertson–Walker–Lemaître metric:

[TABLE]

where $a(t)$ is the scale factor of the Universe.

The particle number is not conserved due to a process of particle creation and annihilation [1, 2]. This process manifests itself, geometro-thermodynamically [3, 4], through the appearance of “creation pressure” $\Pi$ in the cumulative energy-momentum tensor $T_{\mu\nu}$ [5, 10, 11]:

[TABLE]

where $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$.

The Friedmann equations are:

[TABLE]

where $H(t)=\dot{a}(t)/a(t)$, the Hubble parameter, will be considered as one of three dynamical variables of a three-component autonomous dynamical system, also involving the particle number density $n(t)$ and the temperature $T(t)$.

Combining (4) and (5), yields:

[TABLE]

The continuity equation for the particles of the perfect fluid is $N^{\mu}_{\phantom{\mu};\mu}=n\Gamma,$ where $N^{\mu}=nu^{\mu}$ is the particle flow vector and $\Gamma$, the particle production rate, is an input quantity in the phenomenological description [5]. In this work, the dynamics of a model with particle production rate [12]:

[TABLE]

where $\beta$ is a constant, will be studied. As will be shown shortly, due to the second law of thermodynamics, one must have $\Gamma>0$ so that the entropy is never decreasing.

With such particle production rate, the particle conservation equation reads off as

[TABLE]

This equation will be further used as one of the evolution equations of the dynamical system.

The energy conservation equation for the van der Waals gas and for the dust are

[TABLE]

respectively.

The separate conservation laws stipulate that there would be no exchange between the two components of the Universe.

The “creation pressure” $\Pi$, in the case of conserved specific entropy $s$ (i.e. entropy per particle, $s=S/N$, where $S$ is the total entropy), is given by [2]:

[TABLE]

Note that the total entropy $S$ is not conserved due to the enlargement of the phase space resulting from the particle production [2].

On the issue of the equivalence of bulk viscosity and matter creation, Calvaõ et al. [3] and Lima et al. [4] argue that the matter creation process, as described by Prigogine [2], can generate the same dynamic behavior as a FRWL universe with bulk viscosity, while the models being quite different from a thermodynamic point of view. Brevik et al. [13] conclude that creation and viscosity concepts do not describe one and the same physical process — it is shown that viscous and creation universes can develop dynamically in the same manner, but the thermodynamic requirement for their identification is violated. The dynamic pressure $\Pi$ in case of bulk viscosity is given by $\Pi=-3\zeta H$, where $\zeta$ is the bulk viscosity co-efficient [3, 4, 13], while in the case of matter creation processes, similar arguments lead to $\Pi=-\alpha n\Gamma/(3H)$, where $\alpha$ is a phenomenological co-efficient, called creation co-efficient, and it is closely related to the creation process — see [3, 4, 13] and the references therein. The adiabaticity of the fluid, namely, the conservation of the specific entropy, $\dot{s}=0$, leads to the dependence on time of the creation co-efficient $\alpha$: one gets $\alpha=(\rho+p)/n$ — see [13] — and with this, $\Pi=-\alpha n\Gamma/(3H)$ becomes the same as (11).

Substituting (11) into (6) gives:

[TABLE]

This equation describes the dynamical evolution of the Hubble parameter and will be the second equation of the dynamical system.

The dynamical equation (8), multiplied across by $a^{3}$, reads off as $dN/dt=(d/dt)(a^{3}n)=a^{3}n\Gamma$. Differentiating separately $N=na^{3}$ with respect to time, using $\dot{a}=aH$ and (8) gives $\dot{N}=3\beta NH^{\alpha}$. Also, from $s=S/N=$ const, one gets $\dot{S}/S=\dot{N}/N=3\beta H^{\alpha}$. Thus, the constant $\beta$ will be taken as positive and $\alpha$ will be taken as a positive even integer. In the analysis, $\alpha$ and $\beta$ will be considered as parameters of the model.

The integrability condition for the Gibbs equation

[TABLE]

is

[TABLE]

The latter can be written as the following thermodynamic identity:

[TABLE]

In thermodynamical variables $n$ and $T$, the time evolution of the energy density is:

[TABLE]

On the other hand, the energy conservation equation for the van der Waals gas can be written as:

[TABLE]

Using the number conservation equation (8) in (16) and equating to (17) gives:

[TABLE]

Expressing $\rho+p$ from (15) and substituting in the above gives the temperature evolution law:

[TABLE]

and third dynamical equation of the system.

In the absence of particle creation or annihilation (i.e. when $\beta=0$), the above reduces to the well known form given in [9, 14, 15].

Using the equation of state (1) for the van der Waals gas,

[TABLE]

one finds $(\partial p/\partial T)_{n}=n(1+An).$ Substituting this and the equation of state into the thermodynamic identity (15) yields:

[TABLE]

This differential equation can be easily integrated:

[TABLE]

In the case of an ideal monoatomic gas with three translational degrees of freedom, the mass density is, approximately, $\rho=n[m_{0}+(3/2)T]$. The expression (22) for $\rho$ should agree with that for an ideal gas when ideal gas limit is applied for the van der Waals gas, that is, when $A$ and $B$ are both set to zero. This gives $\phi(T)=m_{0}+(3/2)T$. Namely, the energy density $\rho$, the number density $n$, and the temperature $T$ of the van der Waals gas are related via

[TABLE]

Thus, $(\partial\rho/\partial T)_{n}=(3/2)n$ and the dynamical system for the case of a van der Waals gas becomes:

[TABLE]

where $p(n,T)=nT(1+An)-Bn^{2}$ and $\rho(n,T)=n[m_{0}+(3/2)T]-Bn^{2}$.

There is a symmetry: dividing (26) by (24) gives:

[TABLE]

and this is independent of $H$.

Equation (27) can be easily integrated to get the temperature evolution law in terms of the particle number density:

[TABLE]

where the positive constant $\tau$ represent a temperature scale and is a third parameter of the model (in addition to $\alpha$ and $\beta$). Note that the temperature is independent of $\alpha$ and $\beta$.

Equation (27) and its solution are the same as the ones encountered in the case of absence of matter creation or annihilation [9].

The temperature can be excluded so that the system can be reduced to a two-component one:

[TABLE]

3 Analysis

There is a global first integral given by:

[TABLE]

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 defined by $(d/dt)K(\vec{x})=\mu(\vec{x})K(\vec{x})$. It is as an invariant, but only on a restricted subset, given by its zero level set [16]. As no trajectory can cross a hyper-surface defined by a second integral, the second integrals “fragment” the phase space into regions with separate dynamics (yet governed by the same dynamical system). For the two-component dynamical system, the ordinate $n=0$ is one such invariant manif old because $(d/dt)n=[3H(\beta H^{\alpha-1}-1)]n$. Similarly, the curve defined by $3H^{2}-\rho=3H^{2}-n[m_{0}+(3/2)T]+Bn^{2}=0$ is another second integral because $(d/dt)(3H^{2}-\rho)=-3H(3H^{2}-\rho).$ It will be called a separatrix — see Figure 1.

There is a value $\tau_{0}$ of $\tau_{0}$ for which the separatrix $3H^{2}-n[m_{0}+(3/2)T]+Bn^{2}=3H^{2}-n[m_{0}+(3/2)\,\tau\,n^{2/3}\,e^{2An/3}]+Bn^{2}=0$ is tangent to the $n$-axis at point, say $n_{0}$ (see Figure 1). Both $\tau_{0}$ and $n_{0}$ can be determined as follows. When $\tau=\tau_{0}$, the separatrix has a minimum at $n_{0}$ and that minimum is [math]. Thus, $(3/2)\,\tau_{0}\,n_{0}^{2/3}\,e^{2An_{0}/3}=Bn_{0}-m$ and $(d/dn)\left[n[m_{0}+(3/2)\,\tau\,n^{2/3}\,e^{2An/3}]-Bn^{2}\right]_{n=n_{0},\tau=\tau_{0}}=0$ with solutions $n_{0}=[2m_{0}A+B+(4m_{0}^{2}A^{2}+20m_{0}AB+B^{2})^{1/2}\,]/(4AB)$ and $\tau_{0}=(2/3)(Bn_{0}-m_{0})n_{0}^{-2/3}\,e^{-2An_{0}/3}$.

The energy density $\rho[n,T(n)]=n[m_{0}+(3/2)\,\tau\,n^{2/3}\,e^{2An/3}]-Bn^{2}>0$ may become negative over a certain range of $n$, depending on the choice of initial conditions, namely, depending on $\tau$. Such trajectories would temporarily violate the weak energy condition and, as this is admissible in phantom cosmology models [17], the validity of the model will not be restricted by this.

The stability matrix $L$ for the two-component dynamical system (29)–(30) is given by:

[TABLE]

There are three types of critical points for the dynamical system. Firstly, one has the critical points $(n^{*},H^{*}=0)$, where $n^{*}$ are the solutions of the equation $p(n^{*})=0$, that is $\tau(An^{*}+1)n^{*^{5/3}}e^{2An^{*}/3}-Bn^{*^{2}}=0$. This can be written as:

[TABLE]

with

[TABLE]

Depending on the parameter $\tau$ (i.e. on the choice of initial conditions), the number of intersection points of these two curves is one (the origin), two [the origin and a point $\widetilde{n^{*}}$ at which $T(n^{*})$ is tangent to $T^{*}(n^{*})$], or three — one of which is the origin and the other two are $\nu_{1,2}^{*}$ which tend to $\widetilde{n^{*}}$ as $\tau\to\widetilde{\tau}$ from below — see Figure 2a.

To determine the value of $\widetilde{\tau},\,$ for which $T(n^{*})=\widetilde{\tau}n^{*^{2/3}}e^{2An^{*}/3}$ is tangent to $T^{*}(n^{*})=Bn^{*}/(An^{*}+1)$, and to also determine the point $\widetilde{n^{*}}$ from the $n^{*}$-axis where these two curves are tangent to each other, consider the following. At point $\widetilde{n^{*}}$, the two curves intersect, i.e. $\widetilde{\tau}\,\widetilde{n}^{*^{2/3}}e^{2A\widetilde{n}^{*}/3}=B\widetilde{n}^{*}/(A\widetilde{n}*+1)$, and, also, the tangents to the two curves coincide, i.e. $[(d/dn^{*})T(n^{*})]_{(n^{*}=\widetilde{n}^{*},\tau=\widetilde{\tau})}=[(d/dn^{*})T^{*}(n^{*})]_{(n^{*}=\widetilde{n}^{*},\tau=\widetilde{\tau})}$. From these two simultaneous equations, one immediately determines that $\widetilde{n}^{*}=(\sqrt{3/2}-1)/A$ and that $\widetilde{\tau}=B\,\widetilde{n}^{*^{1/3}}e^{-2A\widetilde{n}^{*}/3}\,(1+A\widetilde{n}^{*})^{-1}=\sqrt{2/3}(\sqrt{3/2}-1)^{1/3}e^{2/3-\sqrt{2/3}}A^{-1/3}B$. (For the numerical example considered, one has $\widetilde{n}^{*}=22.47$ and $\widetilde{\tau}=19.84$.)

Focusing firstly on the case of $\alpha=2$, the components of the stability matrix at the critical point $(n^{*},H^{*}=0)$ are: $L^{*}_{11}=0,\,\,\,L^{*}_{12}=-3n^{*}$,

[TABLE]

The eigenvalues at this point are:

[TABLE]

Note that the point at which $L_{21}^{*}$ becomes zero, that is, the point at which the smaller eigenvalue $\lambda^{*}_{2}$ changes sign, is exactly equal to the determined earlier $\widetilde{n}^{*}=(\sqrt{3/2}-1)/A$ — the point at which $T(n^{*})$ is tangent to $T^{*}(n^{*})$ when $\tau=\widetilde{\tau}$. With the decrease of $\tau$ in $T(n^{*})$, the point at which the graphs of $T(n^{*})$ and $T^{*}(n^{*})$ are tangent bifurcates into two intersection points: $\nu_{1,2}^{*}$ (see Figure 2a). Thus, for critical points to the left of $\widetilde{n}^{*}$, where $L_{21}^{*}$ is positive, the eigenvalues are both positive or with positive real parts, while for critical points to the right of $\widetilde{n}^{*}$, where $L_{21}^{*}$ is negative, the eigenvalues are real with $\lambda_{1}^{*}$ being positive, while $\lambda_{2}^{*}$ — negative (see Figures 2b and 2c). In view of this, given that the eigenvalue $\lambda_{1}^{*}$ is always positive over the range of $n^{*}$ where it is real or it always has positive real part over the range of $n^{*}$ where it is complex, critical points $(n^{*},H^{*}=0)$ are never stable.

The eigenvalues $\lambda_{1,2}^{*}$ will be real numbers when the determinant $L_{22}^{*^{2}}-12L_{21}^{*}n^{*}$ is non-negative. This happens when $\beta>\sqrt{8B}/m_{0}=0.09$. When $\beta<\sqrt{8B}/m_{0}$, the eigenvalues will be complex numbers when $n^{*}$ is in the interval $0<n^{*}<N_{0}^{*}$, where $N_{0}^{*}$ (which is smaller than $\widetilde{n}^{*}$) is the only positive root of $L_{22}^{*^{2}}-12L_{21}^{*}n^{*}=0$:

[TABLE]

For example, for $\beta=0.02$, one has $N_{0}^{*}=19.70$, while for $\beta=0.05$, the value of $N_{0}^{*}$ is $9.18$.

Given that to the left of $\widetilde{n}^{*}$ one has $L_{21}^{*}>0$, the eigenvalues will have positive real parts (Figure 2b). Such critical points are unstable and the trajectories near them are unwinding spirals (Figures 4 and 5).

When $N_{0}^{*}<n^{*}<\widetilde{n}^{*}=(\sqrt{3/2}-1)/A=22.47$ and $\beta<\sqrt{8B}/m_{0}=0.09$, the eigenvalues are both real and positive (Figure 2b). The critical points are unstable nodes (Figures 4 and 5). When $n^{*}>(\sqrt{3/2}-1)/A$, the eigenvalues are both real — one positive and one negative (Figure 2b) and one has saddles (Figures 4 and 5).

When $\beta>\sqrt{8B}/m_{0}=0.09$, the eigenvalues $\lambda_{1,2}^{*}$ are both real and positive for $0<n^{*}<\widetilde{n}^{*}=(\sqrt{3/2}-1)/A=22.47$ (Figure 2c). These critical points are unstable nodes (Figure 6). And, finally, for $n^{*}>\widetilde{n}^{*}=(\sqrt{3/2}-1)/A$, the eigenvalues are both real with $\lambda_{1}^{*}$ being positive and $\lambda_{2}^{*}$ – negative (Figure 2c). Such critical points are saddles (Figure 6).

The difference between the cases of $\alpha=2$ and $\alpha>2$ is in the 22-component ($\partial f_{2}/\partial H$) of the stability matrix $L$. At the critical point $(n^{*},H^{*}=0)$, it is not zero when $\alpha=2$ and zero when $\alpha>2$. Consider next the $\alpha=4$ dynamical system and denote the stability matrix by $L^{(4)}$ in this case. One has $L^{(4)^{*}}_{22}=0$ and the eigenvalues at the critical points $(n^{*},H^{*}=0)$ are given by

[TABLE]

The eigenvalues are purely imaginary, $\lambda_{1,2}^{(4)*}=\pm i\omega$, when $A^{2}n^{*^{2}}+2An^{*}-1/2<0$. That is, for $n^{*}$ from zero to $(\sqrt{3/2}-1)/A$ — exactly the point $\widetilde{n}^{*}$ at which $T(n^{*})=T^{*}(n^{*})$ when $\tau=\widetilde{\tau}=\sqrt{2/3}(\sqrt{3/2}-1)^{1/3}e^{2/3-\sqrt{2/3}}A^{-1/3}B$.

For values of $n^{*}$ above $\widetilde{n}^{*}=(\sqrt{3/2}-1)/A$, the eigenvalues are purely real: $\lambda_{1,2}^{(4)*}=\pm q$.

For $\tau>\widetilde{\tau}$, the curves $T(n^{*})$ and $T^{*}(n^{*})$ intersect only at the origin, thus critical points $(n^{*},H^{*}=0)$ do not exists (see Figure 2a).

For $\tau<\widetilde{\tau}$, the curves $T(n^{*})$ and $T^{*}(n^{*})$ intersect, except at the origin, at points $\nu_{1,2}^{*}$ (see Figure 2a again) and the intersection points $\nu_{1,2}^{*}$ are on either side of $\widetilde{n}^{*}$. Thus, at $n^{*}=\nu_{1}^{*}$, the eigenvalues $\lambda_{1,2}^{(4)*}$ are purely imaginary while, at $n^{*}=\nu_{2}^{*}$, they are purely real (with opposite signs) and the corresponding critical points are saddles.

The behaviour of the trajectories near the critical points $(n^{*},H^{*}=0)$ for which the eigenvalues are purely imaginary, namely, for $n^{*}<(\sqrt{3/2}-1)/A$, are studied with the help of centre-manifold theory [18] in the Appendix. One finds that all critical points with purely imaginary eigenvalues are unstable — the trajectories near them are unwinding spirals [18] — see Figures 7a and 7c.

The origin is also a critical point. The analysis of its behaviour is done by expanding the dynamical equations near the origin and retaining only the leading terms. For any $\alpha\geq 2$, one has:

[TABLE]

Consider again the separatrix $3H^{2}-\rho=0$, i.e. the second integral given by $3H^{2}-n\left[m_{0}+(3/2)\tau n^{2/3}e^{2An/3}\right]+Bn^{2}=0$. Along the separatrix and near the origin, one has $3H^{2}=m_{0}n\,\,+$ smaller terms. Then, the equations of the dynamical system in terms of powers of $H$ not higher than 3, reduce to $\dot{n}=-3nH$ and $\dot{H}=-(3/2)H^{2}$. The solutions are:

[TABLE]

where $\sigma=$ sgn $\!\!(H_{0})$.

In view of the continuity, the behaviour of the trajectories near the separatrix will be the same as the behaviour along the separatrix. For the trajectories in the upper half-plane, one will therefore have $n(t)\simeq 1/t^{2}$, while for those in the lower half-plane, $n(t)$ will increase with time. Similarly, $H$ will decay to zero ($H\simeq 1/t$) for trajectories in the upper half-plane or $H$ will decrease with time for trajectories in the lower half-plane.

The origin will attract trajectories from the upper half-plane and repel those from the lower half-plane.

There are other critical points for the $\alpha\geq 2$ dynamical system $\dot{n}=-3nH(1-\beta H^{\alpha-1}),\,\,$ $\dot{H}=-(3/2)H^{2}-(1/2)(1-\beta H^{\alpha-1})p[n,T(n)]+(1/2)\beta H^{\alpha-1}\rho[n,T(n)]$.

Clearly, if $1-\beta H^{\alpha-1}=0$, then $\dot{n}=0$ immediately and for the points $(n^{**},H^{**})$ of the separatrix $3H^{**^{2}}=\rho(n^{**})$, for which $H^{**}=\beta^{\frac{1}{1-\alpha}}$, one will also have $\dot{H}=0$, provided that $n^{**}$ are the solutions of $m_{0}n^{**}+(3/2)n^{**}T(n^{**})-Bn^{**^{2}}-3\beta^{\frac{2}{1-\alpha}}=0$ which can be written as:

[TABLE]

with

[TABLE]

Thus, such $(n^{**},H^{**}=\beta^{\frac{1}{1-\alpha}})$ are critical points for the $\alpha\geq 2$ dynamical system, in addition to the critical points $(n^{*},H^{*}=0)$ and the origin. For these critical points one has:

[TABLE]

and this is greater than zero for all $n^{**}$.

Since the critical points $(n^{**},H^{**})$ are on the separatrix, one should solve the equation for the trajectory reaching or moving away from such critical point firstly while on the separatrix itself. Substituting $\rho=3H^{2}$ into the dynamical equation for $H$ yields:

[TABLE]

and then, expanding about $H^{**}$, gives:

[TABLE]

where $\kappa=(1/2)\beta(\alpha-1)H^{**^{\alpha-2}}[3H^{**^{2}}+p^{**}]\!=\!(1/2)\beta(\alpha-1)H^{**^{\alpha-2}}[(5+2An^{**})\beta^{\frac{2}{1-\alpha}}$ $+(2n^{**}/3)(1+An^{**})(Bn^{**}-m_{0})-Bn^{**^{2}}]=$ const.

The solution along the separatrix near the critical point $(n^{**},H^{**})$ is therefore:

[TABLE]

The sign of $\kappa$ is important. When $\kappa>0$, in order to get $H\to H^{**}$, it is necessary to have $t\to-\infty,$ i.e. the separatrix in this case is an unstable curve of a saddle or the critical point is an unstable node. When $\kappa<0$, one has $H\to H^{**}$ as $t\to\infty,$ i.e. the separatrix in this case is a stable curve of a saddle or the critical point is a stable node. In view of the continuity, trajectories close to the separatrix will exhibit similar behaviour.

The function $T^{**}(n^{**})$ has a minimum at $\beta^{\frac{1}{1-\alpha}}\sqrt{3/B}$. When $\beta$ equals $\beta_{0}\!=\!(12B/m_{0}^{2})^{\frac{\alpha-1}{2}}\!$, this minimum will occur at $n_{0}^{**}$ from the $n^{**}$-axis: $n_{0}^{**}=\beta^{\frac{1}{1-\alpha}}\sqrt{3/B}=m_{0}/(2B)$. For values of $\beta<\beta_{0}$, the graph of $T^{**}(n^{**})$ is entirely above the $n^{**}$-axis, while for $\beta>\beta_{0}$, the function $T^{**}(n^{**})$ has zeros given by $\nu_{1,2}^{**}=[m_{0}/(2B)]\,[1\pm(1-12Bm_{0}^{-2}\beta^{\frac{2}{1-\alpha}})^{\frac{1}{2}}]$ — see Figure 3a. When $\alpha=2$, for the numerical example considered one has $\beta_{0}=0.1095$, while for $\alpha=4,$ the corresponding value is $\beta_{0}=0.0013$.

Depending on the parameters $\beta$ and $\tau$, the number of intersection points of the curves $T(n^{**})$ and $T^{**}(n^{**})$ is one, two, or three — see Figures 3a and 3b. At some value $\hat{\tau}$ of $\tau$, for any given $\beta$, the curves $T(n^{**})$ and $T^{**}(n^{**})$ are tangent to each other at point, say $\hat{n}^{**}$. At this point, the tangents to the two curves coincide, thus one has the following two simultaneous equations: $T(\hat{n}^{**})=T^{**}(\hat{n}^{**})$ and $(d/dn^{**})[T(n^{**})]_{(n^{**}=\hat{n}^{**},\tau=\hat{\tau})}=(d/dn^{**})[T^{**}(n^{**})]_{(n^{**}=\hat{n}^{**},\tau=\hat{\tau})}.$ The solution of this system is $\hat{n}^{**}$, which satisfies

[TABLE]

and $\hat{\tau}$ given by

[TABLE]

When $\tau<\hat{\tau}$, that is, when points $\hat{\nu}_{1,2}^{**}$ exist, one has $\hat{\nu}_{1}^{**}$ to the left of $\hat{n}^{**}$ and $\hat{\nu}_{2}^{**}$ to the right of $\hat{n}^{**}$.

The components of the stability matrix $L$ at the critical points $(n^{**},H^{**}=\beta^{\frac{1}{1-\alpha}})$ are: $L^{**}_{11}=0,\,\,L^{**}_{12}=3(\alpha-1)n^{**}$,

[TABLE]

The eigenvalues are always real:

[TABLE]

Given that $\lambda_{1}^{**}<0$, the critical points $(n^{**},H^{**}=\beta^{\frac{1}{1-\alpha}})$ will be stable if $\lambda_{2}^{**}$ is negative, that is, if $L^{**}_{21}<0$ or if

[TABLE]

Otherwise, the critical points $(n^{**},H^{**}=\beta^{\frac{1}{1-\alpha}})$ will be saddles.

Four curves $T^{**}(n^{**})$ with different $\beta$ are shown on Figure 3c, together with the curve $Q(n^{**})$ which starts at point $(0,-2m_{0}/5)$, crosses the $n^{**}$-axis at $n_{0}^{**}=m_{0}/(2B)$ and has a horizontal asymptote at $2B/A$. When $\beta>\beta_{0}=(12B/m_{0}^{2})^{(\alpha-1)/2}$, the curve $T^{**}(n^{**})$, marked with (i) on Figure 3c, intersects the $n^{**}$-axis at points $\nu_{1,2}^{**}$. The $n^{**}$-coordinates of the intersection point of $T^{**}(n^{**})$ with the curve $Q(n^{**})$ are $\sigma_{1,2}^{**}$. As, while negative, $T^{**}(n^{**})$ cannot intersect the strictly positive $T(n^{**})$, no critical points $(n^{**},H^{**}=\beta^{1/(1-\alpha)})$ can exist for $T^{**}(n^{**})<0$. Thus, stable critical points for $\beta>\beta_{0}$ exist in the interval $\nu_{2}^{**}<n^{**}<\sigma_{2}^{**}$ — where the non-negative $T^{**}(n^{**})$ is smaller than $Q(n^{**})$. When $\beta=\beta_{0}$, the curve $T^{**}(n^{**})$, marked with (ii) on Figure 3c, is tangent to the $n^{**}$-axis at point $\chi_{1}^{**}=n_{0}^{**}$ — the point at which $Q(n^{**})$ crosses the abscissa. This curve intersects the curve $Q(n^{**})$ further — at point $\chi_{2}^{**}$. Critical points for which $\chi_{1}^{**}<n^{**}<\chi_{2}^{**}$ are stable.

There is a value of $\beta$, say $\beta_{Q}$, for which, at certain $n_{Q}^{**}$ from the $n^{**}$-axis, the curve $T^{**}(n^{**})$ is tangent to the $\beta$-independent curve $Q(n^{**})$. That is, at $n_{Q}^{**}$, the two functions are equal, $T^{**}(n_{Q}^{**})=Q(n_{Q}^{**})$, and their first derivatives are also equal, $(d/dn^{**})[T(n^{**})]_{(n^{**}=n_{Q}^{**},\,\beta=\beta_{Q})}=(d/dn^{**})[Q(n^{**})]_{(n^{**}=n_{Q}^{**},\,\beta=\beta_{0})}.$ Thus, $n_{Q}^{**}$ is found, for any $\alpha\geq 2$, as the only positive root of the cubic equation

[TABLE]

For the numerical example considered, one gets $n_{Q}^{**}=45.4587$ and, hence, $\beta_{Q}=0.03426$ for $\alpha=2$ and $\beta_{Q}=0.00004$ for $\alpha=4$.

When $\beta_{Q}<\beta<\beta_{0}$, curve $T^{**}(n^{**})$, marked with (iii) on Figure 3c, never intersects the $n^{**}$-axis. It intersects the curve $Q(n^{**})$ at points with $n^{**}$ coordinates given by $\xi_{1,2}^{**}$. Critical points with $\xi_{1}^{**}<n^{**}<\xi_{2}^{**}$ are stable. Finally, when $\beta<\beta_{Q}$, curve $T^{**}(n^{**})$, marked with (iv) on Figure 3c, never intersects the $n^{**}$-axis or the curve $Q(n^{**})$. There are no stable critical points in this case.

For the dynamical system in the case of $\alpha=2$, three sub-cases are considered: $\beta=0.02$ (Figure 4), $\beta=0.05$ (Figure 5), and $\beta=0.1$ (Figure 6). With these, all qualitatively different possibilities are analyzed. The case of $\alpha=4$ is similar — see Figure 7 where some representative cases are shown. The two Tables at the end should also be considered as all possibilities for the model parameters are summarized there and references are given to the corresponding Figures.

Many of the trajectories exhibit inflationary regime (Figures 4 – 7). This happens in the upper half-plane ($H>0$) and while $H$ is increasing ($\dot{H}>0$), thus $\ddot{a}>0$. The un-physical trajectories that diverge to $(n\to\infty,H\to\infty)$ have eternal inflation, while all other trajectories with inflation, after exiting their inflationary regimes, either extinguish at the origin $(n\to 0,H\to 0)$ in infinite time (Big Freeze); or at a stable critical point in infinite time; or diverge to a Big Crunch: $(n\to\infty,H\to-\infty)$.

4 Conclusions

A cosmological model with two matter components — dust and gas with van der Waals equation of state has been examined. In addition, the model includes a particle production term, proportional to a constant power, $\alpha$, of the Hubble parameter $H$. Models with $\alpha=2$ and $\alpha=4$ are studied in detail. However, the presented analysis can easily be extended to an arbitrary integer $\alpha$ (the special case of $\alpha=1$ deserves a special attention and will be provided elsewhere).

The time-evolution of the model is given by a nonlinear dynamical system of three equations: for the particle number density $n$, the Hubble parameter $H$ and the temperature $T$. This system admits a global first integral, which explicitly gives $T$ as a function of $n$ and one of the van der Waals gas parameters. Hence, the system is reduced to a two-component one: in the two dimensional $n$–$H$ phase space. The system exhibits a complex behavior which is influenced by the presence of the several model parameters. This behaviour is examined in detail using the phase-plane analysis for all possible parameter choices. The two second integrals of the system are represented by curves which separate the phase space into domains which can not be crossed by the trajectories. The full classification of the critical points is presented in the two provided tables. It is shown that the critical points can not be reached in a finite time (the stable critical points can only be reached for $t\to\infty$, the unstable critical points can be reached only for $t\to-\infty.$). The critical points provide important information about the large-time behaviour of the system. This includes both the distant future ($t\to\infty$) or the distant past ($t\to-\infty$). For example, considering trajectories which end at the origin, i.e. $(n,H)\to(0,0)$ as $t\to\infty$, from (8) one has $\dot{n}=-3nH$ asymptotically when $\alpha\geq 2$ and taking into account (10), it follows that $d\rho_{d}/dn=\rho_{d}/n,$ or $\rho_{d}=Cn$ for some constant $C$. Then

[TABLE]

when $(n,H)\to(0,0).$ Therefore, the ratio between the two fractions approaches a constant.

In the case of high particle creation $n\to\infty$ and $H\to\infty$ (this can be viewed as a critical point at infinity), in the distant past or future, i.e. when $t\to\pm\infty$, the asymptotic equations are

[TABLE]

giving $H^{2}=(1/2)\,\tau\,n^{5/3}\,e^{2An/3}\,+\,$ lower-order terms. Substituting this asymptotic form of $H^{2}$ into the Friedmann equation (5) yields:

[TABLE]

In other words $\rho_{d}/\rho\to 0,$ which means that in this case the dust component becomes negligible and all trajectories are drawn in the neighbourhood of the separatrix $3H^{2}=\rho$ as $t\to\pm\infty$.

Finally, sets of initial values can be identified for which the corresponding trajectories exhibit inflationary behavior.

[FIGURE:]

[FIGURE:]