Van der Waals universe with adiabatic matter creation

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.

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,5] 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 [6] 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 [7]. 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, [8,9,10,11,12], 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 [13]. 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

$$\begin{aligned} p = n T [1 + n F(T)], \end{aligned}$$

(1)

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.¹

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

$$\begin{aligned} ds^2= & {} g_{\mu \nu } dx^\mu dx^\nu \nonumber \\= & {} dt^2 - a^2(t) [dr^2 + r^2 (d \theta ^2 + \sin ^2 \theta \, d \phi ^2)], \end{aligned}$$

(2)

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,3]. This process manifests itself, geometro-thermodynamically [4, 5], through the appearance of “creation pressure” $\Pi$ in the cumulative energy-momentum tensor $T_{\mu \nu }$ [6, 14, 15]:

$$\begin{aligned} T_{\mu \nu } = (\rho + \rho _d + p + \Pi ) \, u_\mu \, u_\nu - (p + \Pi ) \, g_{\mu \nu }, \end{aligned}$$

(3)

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:

$$\begin{aligned}&\frac{\ddot{a}}{a} = - \frac{1}{6} [ \rho _d + \rho + 3 (p + \Pi ) ], \end{aligned}$$

(4)

$$\begin{aligned}&H^2 = \frac{1}{3} (\rho _d + \rho ), \end{aligned}$$

(5)

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:

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

(6)

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

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

(7)

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

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

(8)

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

$$\begin{aligned}&{\dot{\rho }} + 3H(\rho + p + \Pi ) = 0, \end{aligned}$$

(9)

$$\begin{aligned}&\dot{\rho _d} + 3H \rho _d = 0, \end{aligned}$$

(10)

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 [3]:

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

(11)

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

On the issue of the equivalence of bulk viscosity and matter creation, Calvaõ et al. [4] and Lima et al. [5] argue that the matter creation process, as described by Prigogine [3], 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. [17] 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 [4, 5, 17], 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 [4, 5, 17] 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 [17] – and with this, $\Pi = - \alpha n \Gamma / (3H)$ becomes the same as (11).

Substituting (11) into (6) gives:

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

(12)

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 = n a^3$ with respect to time, using ${\dot{a}} = a H$ and (8) gives ${\dot{N}} = 3 \beta N H^\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

$$\begin{aligned} T ds = d \left( {\rho \over n} \right) + p \, d \left( {1\over n} \right) = -\left( {\rho + p \over n^2}\right) \,\, dn + {1 \over n} \,\, d\rho .\nonumber \\ \end{aligned}$$

(13)

is

$$\begin{aligned} n \left( \frac{\partial T}{\partial n}\right) _\rho + (\rho + p)\left( \frac{\partial T}{\partial \rho }\right) _n = T \left( \frac{\partial p}{\partial \rho }\right) _n. \end{aligned}$$

(14)

The latter can be written as the following thermodynamic identity:

$$\begin{aligned} \rho + p = T \left( \frac{\partial p}{\partial T}\right) _n + n \left( \frac{\partial \rho }{\partial n}\right) _T. \end{aligned}$$

(15)

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

$$\begin{aligned} {\dot{\rho }}(n, T) = \left( \frac{\partial \rho }{\partial n}\right) _T {\dot{n}} + \left( \frac{\partial \rho }{\partial T}\right) _n \,\, {\dot{T}}. \end{aligned}$$

(16)

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

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

(17)

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

$$\begin{aligned}&- 3 (\rho + p) H(1 - \beta H^{\alpha -1}) \nonumber \\&\quad = -3 n H(1 - \beta H^{\alpha -1}) \, \left( \frac{\partial \rho }{\partial n}\right) _T + \left( \frac{\partial \rho }{\partial T}\right) _n \,\, {\dot{T}}. \end{aligned}$$

(18)

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

$$\begin{aligned} {\dot{T}}= & {} - 3 H (1 - \beta H^{\alpha -1}) T \left( \frac{\partial p}{\partial \rho } \right) _n\nonumber \\= & {} - 3 H (1 - \beta H^{\alpha -1}) T \frac{\left( \frac{\partial p}{\partial T} \right) _n}{\left( \frac{\partial \rho }{\partial T} \right) _n} \end{aligned}$$

(19)

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 [13, 18, 19].

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

$$\begin{aligned} p(n,T) = n T (1 + An) - B n^2, \end{aligned}$$

(20)

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

$$\begin{aligned} \left[ \frac{\partial }{\partial n} \left( \frac{\rho }{n} \right) \right] _T = - B. \end{aligned}$$

(21)

This differential equation can be easily integrated:

$$\begin{aligned} \rho = n [ \phi (T) - B n]. \end{aligned}$$

(22)

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

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

(23)

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

$$\begin{aligned} {\dot{n}}= & {} - 3 n H (1 - \beta H^{\alpha - 1}), \end{aligned}$$

(24)

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

(25)

$$\begin{aligned} {\dot{T}}= & {} - 2 (1 + A n) H (1 - \beta H^{\alpha - 1}) T, \end{aligned}$$

(26)

where $p(n, T) = n T ( 1 + A n ) - B n^2$ and $\rho (n, T) = n [m_0 + (3/2) T] - B n^2$. There is a symmetry: dividing (26) by (24) gives:

$$\begin{aligned} \frac{dT}{dn} = \frac{2T(1+An)}{3n}> 0 \, \text{ as } n > 0\text{, } \end{aligned}$$

(27)

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:

$$\begin{aligned} T(n) = \tau \, n^{\frac{2}{3}} \, e^{\frac{2 A n}{3}}, \end{aligned}$$

(28)

where the positive constant $\tau$ (not to be confused with proper time) represents 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 [13].

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

$$\begin{aligned}&{\dot{n}} \equiv f_1(n, H) = 3 n H (\beta H^{\alpha -1} - 1), \end{aligned}$$

(29)

$$\begin{aligned} {\dot{H}}\equiv & {} f_2(n, H) = \, - \frac{3}{2} \, H^2+ \frac{1}{2} \, \tau \, n^{\frac{5}{3}} \, e^{\frac{2An}{3}}\nonumber \\&\times \, \left[ (\beta H^{\alpha -1} - 1) \left( An + \frac{5}{2}\right) + \frac{3}{2} \right] \nonumber \\&+ \, \frac{1}{2} \, \beta \, (m_0 - 2Bn)\, n \, H^{\alpha -1} + \frac{1}{2} \, B \, n^2. \end{aligned}$$

(30)

3 Analysis

There is a global first integral given by:

$$\begin{aligned} I(n,T) = T \, n^{-\frac{2}{3}} \, e^{-\frac{2 A n}{3}} = \tau \text{= } \text{ const } > 0. \end{aligned}$$

(31)

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 [20]. 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 manifold 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 Fig. 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 Fig. 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 0. 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 = [ 2 m_0 A + B + (4 m_0^2 A^2 + 20 m_0 A B + 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 [21], 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:

Fig. 1

The separatrix $3H^2 - n[m_0 + (3/2) \, \tau \, n^{2/3} \, e^{2An/3}] + Bn^2 = 0$ is an open curve when $\tau > \tau _0 = 14.78$ for a van der Waals gas with parameters $A = 0.01$ and $B = 10$ and $m_0$, the typical mass of a representative particle, taken as 100. When $\tau < \tau _0$, the separatrix has a loop at low n and an open part at high n. When $\tau > \tau _0$, the trajectories to the right of the separatrix are those for dust component $\rho _d < 0$, while those above or below it are with $\rho _d > 0$. On the separatrix itself, $\rho _d = 0$. When $\tau < \tau _0$, the trajectories to the right of the open curve and those inside the loop are with $\rho _d <0$ while all others have $\rho _d > 0$. The curve with $\tau = \tau _0$ is tangent to the abscissa at $n_0 = \bigl ( 2 m_0 A + B + \sqrt{4 m_0^2 A^2 + 20 m_0 A B + B^2} \, \bigr ) / (4AB) = 73.59$. The energy density $\rho [n, T(n)] = n[m_0 + (3/2) \, \tau \, n^{2/3} \, e^{2An/3}] - Bn^2$ is positive for all values of n if $\tau > \tau _0$

Fig. 2

The critical points of the type $(n^*, H^* = 0)$ have $n^*$ determined by the intersection points of the curve $T(n^*) = \tau n^{*^{2/3}} e^{2An^*/3}$ with the curve $T^*(n^*) = B n^*/(1 + A n^*)$. These are: only the origin, when $\tau > {\widetilde{\tau }} = \sqrt{2/3} (\sqrt{3/2}-1)^{1/3} e^{2/3-\sqrt{2/3}} A^{-1/3} B$; the origin and ${\widetilde{n}}^* = (\sqrt{3/2} - 1)/A$ when $\tau = {\widetilde{\tau }}$; and the origin and $\nu _{1,2}^*$ when $\tau < {\widetilde{\tau }}$ (with $\nu _{1,2}^* \rightarrow {\widetilde{n}}^*$ when $\tau \rightarrow {\widetilde{\tau }}$ from below) – a. The eigenvalues $\lambda _{1,2}^*$ at critical points to the left of ${\widetilde{n}}^* = (\sqrt{3/2} - 1)/A$ are both positive or with positive real parts (depending on $\beta$), while at critical points to the right of ${\widetilde{n}}^*$ the eigenvalues are real with $\lambda _1^*$ being positive, while $\lambda _2^*$ – negative (see b and c). Given that one of the eigenvalues is always positive or has positive real part, critical points $(n^*, H^* = 0)$ are never stable

$$\begin{aligned} L_{11}= & {} \frac{\partial f_1}{\partial n} = 3 H (\beta H^{\alpha - 1} - 1), \end{aligned}$$

(32)

$$\begin{aligned} L_{12}= & {} \frac{\partial f_1}{\partial H} = 3 n (\alpha \beta H^{\alpha - 1} - 1), \end{aligned}$$

(33)

$$\begin{aligned} L_{21}= & {} \frac{\partial f_2}{\partial n} = \frac{1}{3} \, \tau \, n^{\frac{2}{3}} \, e^{\frac{2An}{3}} \nonumber \\&\times \left[ \left( An {+} \frac{5}{2} \right) ^2 (\beta H^{\alpha {-} 1} {-} 1) {+} \frac{3}{2} \, \beta A n H^{\alpha - 1} + \frac{15}{4} \right] \nonumber \\&+ \,\, \frac{1}{2} \, \beta m_0 H^{\alpha -1} + (1 - 2 \beta H^{\alpha -1})Bn, \end{aligned}$$

(34)

$$\begin{aligned} L_{22}= & {} \frac{\partial f_2}{\partial H}= -3 H + \frac{1}{2} \, \beta \, (\alpha - 1) \nonumber \\&\times \left[ \left( An + \frac{5}{2} \right) \, \tau \, n^{\frac{5}{3}} \, e^{\frac{2An}{3}} + (m_0 - 2Bn)\,n \right] H^{\alpha -2}.\nonumber \\ \end{aligned}$$

(35)

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:

$$\begin{aligned} T(n^*) = T^*(n^*) \end{aligned}$$

(36)

with

$$\begin{aligned} T^*(n^*) = \frac{B n^*}{A n^* + 1}. \end{aligned}$$

(37)

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 \rightarrow {\widetilde{\tau }}$ from below – see Fig. 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^*$,

$$\begin{aligned} L_{21}^*= & {} - \frac{1}{3} \, \frac{Bn^*}{An^* + 1} \, \left( A^2 {n^{*^2}} + 2 A n^* - \frac{1}{2} \right) , \text{ and } \end{aligned}$$

(38)

$$\begin{aligned} L_{22}^*= & {} \frac{1}{2} \, \beta \, \frac{n^*}{An^* + 1} \, \left[ - A B {n^{{*}^{2}}} + \left( m_0 A + \frac{B}{2} \right) n^* + m_0 \right] .\nonumber \\ \end{aligned}$$

(39)

The eigenvalues at this point are:

$$\begin{aligned} \lambda _{1,2}^* = \frac{1}{2} L_{22}^* \,\, \pm \,\, \frac{1}{2} \sqrt{L_{22}^{*^2} - 12 L_{21}^* n^*}. \end{aligned}$$

(40)

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 Fig. 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 Fig. 2b, c). 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} - 12 L_{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} - 12 L_{21}^* n^* = 0$:

$$\begin{aligned}&m_0^2 \beta ^2 - 8B + \left[ (m_0 \beta ^2 + 24A)B + 2 m_0^2 \beta ^2 A \right] n^* \nonumber \\&\quad +\left[ \frac{1}{4} B^2 \beta ^2 - (m_0 \beta ^2 + 48A)AB + m_0^2 \beta ^2 A^2\right] n^{*^2} \nonumber \\&\quad + \left[ A B^2 \beta ^2 + (2 m_0 \beta ^2 - 16 A) A^2 B \right] n^{*^3} {+} A^2 B^2 \beta ^2 n^{*^4} = 0.\nonumber \\ \end{aligned}$$

(41)

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 (Fig. 2b). Such critical points are unstable and the trajectories near them are unwinding spirals (Figs. 4, 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 (Fig. 2b). The critical points are unstable nodes (Figs. 4, 5). When $n^* > (\sqrt{3/2} - 1)/A$, the eigenvalues are both real – one positive and one negative (Fig. 2b) and one has saddles (Figs. 4, 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$ (Fig. 2c). These critical points are unstable nodes (Fig. 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 (Fig. 2c). Such critical points are saddles (Fig. 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

$$\begin{aligned} \lambda _{1,2}^{(4) *} = \pm \sqrt{\frac{Bn*}{An^* + 1} \left( A^2 n^{*^2} + 2An^* - \frac{1}{2} \right) } \end{aligned}$$

(42)

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 Fig. 2a).

For $\tau < {\widetilde{\tau }}$, the curves $T(n^*)$ and $T^*(n^*)$ intersect, except at the origin, at points $\nu _{1,2}^*$ (see Fig. 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 [22] in the Appendix. One finds that all critical points with purely imaginary eigenvalues are unstable – the trajectories near them are unwinding spirals [22] – see Fig. 7a, c.

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 \ge 2$, one has:

$$\begin{aligned}&{\dot{n}} = - 3nH + 3 \beta n H^{\alpha }, \end{aligned}$$

(43)

$$\begin{aligned}&{\dot{H}} = -\frac{3}{2} H^2 - \frac{1}{2} \tau n^{\frac{5}{3}} + \frac{1}{2} \beta m_0 n H^{\alpha - 1} + \frac{1}{2} Bn^2 + \cdots .\nonumber \\ \end{aligned}$$

(44)

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:

$$\begin{aligned}&n(t) = \frac{n_0}{\left[ 1 + \frac{1}{2} \sigma \sqrt{3 m_0 n_0}(t-t_0)\right] ^2}, \end{aligned}$$

(45)

$$\begin{aligned}&H(t) = \frac{H_0}{1 + \frac{3}{2} H_0 (t - t_0)}, \end{aligned}$$

(46)

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.

Fig. 3

Determination of the critical points of the type $(n^{**}, H^{**} = \beta ^{\frac{1}{1-\alpha }})$ for the $\alpha \ge 2$ dynamical system. The loci $n^{**}$ of the critical points are the solutions of $T^{**}(n^{**}) = T(n^{**})$ – a and b. c Shows where stable critical points of the type $(n^{**}, H^{**} = \beta ^{\frac{1}{1-\alpha }})$ can be found

Fig. 4

The case of $\alpha = 2$ with $\beta = 0.02$. As $\beta < \sqrt{8B}/m_0 = 0.0894$, the eigenvalues $\lambda _{1,2}^*$ are both complex with positive real parts for $0< n^* < N_0^* = 19.70$. The trajectories near them are unwinding spirals (see a and b). For values of $n^*$ between $N_0^* = 19.70$ and $(\sqrt{3/2}-1)/A = 22.47$, the eigenvalues are both real and positive. The critical points are unstable nodes. Finally, when $n^* > (\sqrt{3/2} - 1)/A = 22.47$, the eigenvalues are both real – one positive and one negative and one has saddles. As $\beta = 0.02 < \beta _Q = 0.0343$, both eigenvalues $\lambda _{1,2}^{**}$ are real and with opposite signs for all $n^{**}$, thus the corresponding critical points are always saddles

Fig. 5

Parts a, b, and c – the case of $\alpha = 2$ with $\beta = 0.05$. As $\beta < \sqrt{8B}/m_0 = 0.0894$, the eigenvalues $\lambda _{1,2}^*$ are both complex with positive real parts for $0< n^* < N_0^* = 9.19$. The trajectories near them are unwinding spirals (see a and d). For values of $n^*$ between $N_0^* = 9.19$ and $(\sqrt{3/2}-1)/A = 22.47$, the eigenvalues are both real and positive. The critical points are unstable nodes. Finally, when $n^* > (\sqrt{3/2} - 1)/A = 22.47$, the eigenvalues are both real – one positive and one negative and one has saddles (see c). In relation to the eigenvalues $\lambda _{1,2}^{**}$, one has $n_1^{**} = 18.27$ and $n_2^{**} = 66.45$. Critical points with $0< n^{**} < n_1^{**} = 18.27$ are with real eigenvalues with opposite signs (saddles, see b, d, f), those with $n^*$ between $n_1^{**} = 18.27$ and $n_2^{**} = 66.45$ are with real and negative eigenvalues (stable nodes, see b), and critical points with $n^*$ above $n_2^{**} = 66.45$ are with real eigenvalues with opposite signs (saddles, see c). Parts d, e, and f – the case of $\alpha = 2$ with $\beta = 0.05$

Fig. 6

Parts a, b, and c – the case of $\alpha = 2$ with $\beta = 0.1$. As $\beta > \sqrt{8B}/m_0 = 0.0894$, the eigenvalues $\lambda _{1,2}^*$ are real for all $n^*$ – they are both positive for $0< n^* <(\sqrt{3/2}-1)/A = 22.47$ (with the corresponding critical points being unstable nodes, see a and d) and positive and negative for $n^* > (\sqrt{3/2} - 1)/A = 22.47$ (with the corresponding critical points being saddles, see b, e). In relation to the eigenvalues $\lambda _{1,2}^{**}$, one has $n_1^{**} = 5.83$ and $n_2^{**} = 72.02$. Critical points with $0< n^{**} < n_1^{**} = 5.83$ are with real eigenvalues with opposite signs (saddles, see a and d), those with $n^*$ between $n_1^{**} = 5.83$ and $n_2^{**} = 72.02$ are with real and negative eigenvalues (stable nodes, see b), and critical points with $n^*$ above $n_2^{**} = 72.02$ are with real eigenvalues with opposite signs (saddles, see c. Parts d, e, and f – the case of $\alpha = 2$ with $\beta = 0.1$)

Fig. 7

Parts a and b: the case of $\alpha = 4$ – some representative cases. Parts c and d: the case of $\alpha = 4$ – some representative cases

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 \ge 2$ dynamical system ${\dot{n}} = - 3 n H (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 $3 H^{**^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^{**}) - B n^{**^2} - 3 \beta ^{\frac{2}{1-\alpha }} = 0$ which can be written as:

$$\begin{aligned} T(n^{**}) = T^{**}(n^{**}) \end{aligned}$$

(47)

with

$$\begin{aligned} T^{**}(n^{**}) = \frac{2}{3} (B n^{**} - m_0) + \frac{2 \beta ^{\frac{2}{1-\alpha }}}{n^{**}}. \end{aligned}$$

(48)

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

$$\begin{aligned} \rho ^{**} \equiv \rho [ n^{**}, T^{**}(n^{**}) ] = 3 \beta ^{\frac{2}{1-\alpha }} \end{aligned}$$

(49)

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:

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

(50)

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

$$\begin{aligned} {\dot{H}}= & {} \frac{1}{2}\beta (\alpha -1)H^{**^{\alpha - 2}}[3H^{**^2}+ p^{**}](H - H^{**}) \nonumber \\= & {} \kappa (H - H^{**}), \end{aligned}$$

(51)

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

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

$$\begin{aligned} \ln \left| \frac{H-H^{**}}{H_0-H^{**}} \right| = \kappa (t-t_0). \end{aligned}$$

(52)

The sign of $\kappa$ is important. When $\kappa > 0$, in order to get $H \rightarrow H^{**}$, it is necessary to have $t \rightarrow -\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 \rightarrow H^{**}$ as $t \rightarrow \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 Fig. 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 Fig. 3a, b. 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

$$\begin{aligned}&2 A B {\hat{n}}^{**^3} - (2 m_0 A + B){\hat{n}}^{**^2} + (6 A \beta ^{\frac{2}{1-\alpha }} - 2m_0) {\hat{n}}^{**} \nonumber \\&\quad + 15 \beta ^{\frac{2}{1-\alpha }} = 0, \end{aligned}$$

(53)

and ${\hat{\tau }}$ given by

$$\begin{aligned} {\hat{\tau }} = \left[ \frac{2}{3} (B {\hat{n}}^{**} - m_0) + \frac{2 \beta ^{\frac{2}{1-\alpha }}}{{\hat{n}}^{**}} \right] {\hat{n}}^{**^{-\frac{2}{3}}} e^{-\frac{2A{\hat{n}}^{**}}{3}}. \end{aligned}$$

(54)

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^{**}$,

$$\begin{aligned} L^{**}_{21}= & {} \frac{1}{2} T^{**}(n^{**}) \left( A n^{**} + \frac{5}{2}\right) - B n^{**} + \frac{m_0}{2} \nonumber \\= & {} \frac{1}{3} \left[ A B n^{** ^2} - (m_0 A + \frac{B}{2}) n^{**} - m_0 \right] \nonumber \\&\quad + \frac{A}{\beta ^{\frac{2}{\alpha - 1}}} + \frac{5}{2 \beta ^{\frac{2}{\alpha - 1}} n^{**}}, \end{aligned}$$

(55)

$$\begin{aligned} L^{**}_{22}= & {} - 3\beta ^{\frac{1}{1 - \alpha }}+ \frac{1}{2} \, (\alpha - 1) \, \beta ^{\frac{1}{\alpha - 1}} \, n^{**} \nonumber \\&\times \, \left[ T^{**}(n^{**}) (A n^{**} + \frac{5}{2}) - 2 B n^{**} + m_0\right] \nonumber \\= & {} - 3\beta ^{\frac{1}{1 - \alpha }} \,\, + \,\, (\alpha - 1) \, \beta ^{\frac{1}{\alpha - 1}} \, n^{**} \, L^{**}_{21}. \end{aligned}$$

(56)

The eigenvalues are always real:

$$\begin{aligned} \lambda _1^{**}= & {} - 3\beta ^{\frac{1}{1 - \alpha }} < 0, \end{aligned}$$

(57)

$$\begin{aligned} \lambda _2^{**}= & {} (\alpha - 1) \, \beta ^{\frac{1}{\alpha - 1}} \, n^{**} \, L^{**}_{21} = \frac{1}{2}\, (\alpha - 1) \, \beta ^{\frac{1}{\alpha - 1}} \, n^{**} \nonumber \\&\times \left[ T^{**}(n^{**}) \left( A n^{**} + \frac{5}{2}\right) - 2 B n^{**} + m_0\right] . \end{aligned}$$

(58)

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

$$\begin{aligned} T^{**}(n^{**}) < Q(n^{**}) \equiv \frac{2 B n^{**} - m_0}{A n^{**} + \frac{5}{2}}. \end{aligned}$$

(59)

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

Four curves $T^{**}(n^{**})$ with different $\beta$ are shown on Fig. 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 Fig. 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 Fig. 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 \ge 2$, as the only positive root of the cubic equation

$$\begin{aligned}&4 A^2 B n_Q^{**^3} - 2 A (m_0 A - 7 B) n_Q^{** ^2} \nonumber \\&\quad - 5 (2 m_0 A + B) n_Q^{**} - 5 m_0 = 0. \end{aligned}$$

(60)

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 Fig. 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 Fig. 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$ (Fig. 4), $\beta = 0.05$ (Fig. 5), and $\beta = 0.1$ (Fig. 6). With these, all qualitatively different possibilities are analyzed. The case of $\alpha = 4$ is similar – see Fig. 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 (Figs. 4, 5, 6, 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 \rightarrow \infty , H \rightarrow \infty )$ have eternal inflation, while all other trajectories with inflation, after exiting their inflationary regimes, either extinguish at the origin $(n \rightarrow 0, H \rightarrow 0)$ in infinite time (Big Freeze); or at a stable critical point in infinite time; or diverge to a Big Crunch: $(n \rightarrow \infty , H \rightarrow -\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\rightarrow \infty$, the unstable critical points can be reached only for $t \rightarrow -\infty .$). The critical points provide important information about the large-time behaviour of the system. This includes both the distant future ($t \rightarrow \infty$) or the distant past ($t \rightarrow -\infty$). For example, considering trajectories which end at the origin, i.e. $(n,H) \rightarrow (0,0)$ as $t \rightarrow \infty$, from (8) one has ${\dot{n}}=-3nH$ asymptotically when $\alpha \ge 2$ and taking into account (10), it follows that $d \rho _d / dn = \rho _d / n,$ or $\rho _d = C n$ for some constant C. Then

$$\begin{aligned} \frac{\rho _d}{\rho }=\frac{Cn}{ n [m_0 + \frac{3}{2} \, \tau \, n^{\frac{2}{3}} \, e^{\frac{2 A n}{3}} - B n ]} \rightarrow \frac{C}{m_0} = \mathrm {const} \end{aligned}$$

(61)

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

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

$$\begin{aligned}&{\dot{n}} = 3 \, \beta \, n \, H^{\alpha }, \end{aligned}$$

(62)

$$\begin{aligned}&{\dot{H}} = \frac{1}{2} \, \beta A \, \tau \, n^{\frac{8}{3}} \, e^{\frac{2An}{3}} \, H^{\alpha -1}, \end{aligned}$$

(63)

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:

$$\begin{aligned} \frac{1}{3} \left( 1 + \frac{\rho _d}{\rho }\right)= & {} \frac{\frac{1}{2} \, \tau \, n ^{\frac{5}{3}} e^{\frac{2An}{3}} + \cdots }{n [m_0 + \frac{3}{2} \, \tau \, n ^{\frac{5}{3}} e^{\frac{2An}{3}} - B n]}\nonumber \\&\rightarrow \frac{1}{3} \quad \text{ when } (n, H) \rightarrow (\infty , \infty ). \end{aligned}$$

(64)

In other words $\rho _d /\rho \rightarrow 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 \rightarrow \pm \infty$.

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

Appendix

1.1 Application of centre-manifold theory to the critical points with purely imaginary eigenvalues

To study 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$, centre-manifold theory [22] is applied. Firstly, the dynamical system is expanded near such $(n^*, H^* =0)$:

(65)

$$\begin{aligned}&{\dot{H}} = {-} \frac{3}{2} H^2 {+} \frac{1}{2} \beta \rho ^* H^3 {+} \frac{Bn^* ( A^2 n^{*^2} {+} 2An^* - \frac{1}{2})}{An^*+1}(n {-} n^*) \nonumber \\&\quad - \frac{1}{9} \, \frac{B( A^3 n^{*^3} + 9 A^2 n^{*^2} + \frac{21}{2} An^* - 2)}{An^*+1}(n - n^*)^2 \nonumber \\&\quad - \frac{1}{162} \, \frac{B( 4 A^4 n^{*^4} + 52 A^3 n^{*^3} + 150 A^2 n^{*^2} + 70 An^* - 5)}{(An^*+1) n^*}\nonumber \\&\quad \times (n - n^*)^3 + \cdots .\nonumber \\ \end{aligned}$$

(66)

Introduce new dynamical variables via: $n - n^* = \theta x$ and $H = \mu y$. Introduce also $\omega = 3n^* \mu / \theta$. Taking $\theta = 1$ and $\mu = (1/3) \sqrt{- B ( A^2 n^{*^2} + 2An^* - \frac{1}{2}) / (A n^* + 1)}$ [note that $\mu$ is real to the left of ${\widetilde{n}}^* = (\sqrt{3/2}-1)/A$ – where the analysis applies]. This yields $\omega ^2 = - B n^* ( A^2 n^{*^2} + 2An^* - \frac{1}{2}) / (A n^* + 1) > 0$.

The dynamical system can then be written as:

$$\begin{aligned} {\dot{x}}= & {} - \omega y + f(x, y), \end{aligned}$$

(67)

$$\begin{aligned} {\dot{y}}= & {} \omega x + g(x, y), \end{aligned}$$

(68)

where

$$\begin{aligned}&f(x, y) = - 3 \mu x y (1 - \beta \mu ^3 y^3) + 3 \beta n^* \mu ^4 y^4 , \end{aligned}$$

(69)

$$\begin{aligned}&g(x, y) = - \frac{3}{2} \mu y^2 + \frac{1}{2} \mu ^2 \beta \rho ^* y^3\nonumber \\&\quad - \frac{1}{9 \mu } \, \frac{B( A^3 n^{*^3} + 9 A^2 n^{*^2} + \frac{21}{2} An^* - 2)}{An^*+1}x^2 \nonumber \\&\quad - \frac{1}{162 \mu } \, \frac{B( 4 A^4 n^{*^4} + 52 A^3 n^{*^3} {+} 150 A^2 n^{*^2} {+} 70 An^* {-} 5)}{(An^*+1) n^*}\nonumber \\&\quad \times \, x^3 + \ldots . \end{aligned}$$

(70)

Then, at the critical point $(x = 0, y = 0)$, the stability parameter a – see (3.4.10) and (3.4.11) in [22] – is:

$$\begin{aligned} a= & {} \frac{1}{16} \, (f_{xxx} + f_{xyy} + g_{xxy} + g_{yyy}) \nonumber \\&+ \, \frac{1}{16 \omega } \, [f_{xy} (f_{xx} + f_{yy}) \nonumber \\&-\, g_{xy} (g_{xx} + g_{yy}) - f_{xx} g_{xx} + f_{yy} g_{yy}] \nonumber \\= & {} \frac{3}{16} \mu ^2 \beta \rho ^* > 0, \end{aligned}$$

(71)

provided $\rho ^* > 0$.

The energy density $\rho ^*$ at the equilibrium points $(n^*, H^* = 0)$ is, as in the case of $\alpha = 2$, non-negative for $n^* < n_0^* = \bigl ( 2 m_0 A + B + \sqrt{4 m_0^2 A^2 + 20 m_0 A B + B^2} \, \bigr ) / (4AB)$. And, given that one always has $n_0^* > {\widetilde{n}}^* = (\sqrt{3/2}-1)/A$, then $\rho ^*$ is positive in the entire region $0< n^* < {\widetilde{n}}^* = (\sqrt{3/2}-1)/A$ – where the eigenvalues are purely imaginary. This, in turn, means that a is positive in that region and thus all critical points with purely imaginary eigenvalues are unstable – the trajectories near them are unwinding spirals [22] – see Fig. 7a, c.