Two-fluid solutions of particle-creation cosmologies

TL;DR

This paper investigates cosmological models with two fluids, one with particle creation and one conserved, deriving a key differential equation for dynamics and exploring various particle creation rates.

Contribution

It introduces a general framework for two-fluid cosmologies with particle creation, deriving a nonlinear differential equation governing their evolution.

Findings

Derived a single nonlinear differential equation for two-fluid cosmologies with particle creation.

Explored specific particle creation rates including constant, inverse Hubble, and exponential forms.

Presented and analyzed singular algebraic solutions and their stability.

Abstract

Cosmological evolution driven incorporating continuous particle creation by the time-varying gravitational field is investigated. We consider a spatially flat, homogeneous and isotropic universe with two matter fluids in the context of general relativity. One fluid is endowed with gravitationally induced `adiabatic' particle creation, while the second fluid simply satisfies the conservation of energy. We show that the dynamics of the two fluids is entirely controlled by a single nonlinear differential equation involving the particle creation rate, $\Gamma (t)$. We consider a very general particle creation rate, $\Gamma (t)$, that reduces to several special cases of cosmological interest, including $\Gamma =$ constant, $% \Gamma \propto 1/H^{n}$ ($n\in \mathbb{N}$), $\Gamma \propto \exp (1/H)$. Finally, we present singular algebraic solutions of the gravitational field equations for the two-fluid particle creation models and discuss their stability.

Figures

Tables

Table 1. Table 1: The table presents different particle creation rates corresponding to different effective fluids.

Mimicking Fluid	Total equation of state ($w_{\mathrm{tot}}$)	Particle creation rate ($\Gamma$)
Radiation	$w_{\mathrm{tot}}=1/3$	$\Gamma =3H\left[1-\frac{4H^2-{\gamma }_1{\rho }_1}{\gamma \left(3H^2-{\rho }_1\right)}\right]$
Dust	$w_{\mathrm{tot}}=0$	$\Gamma =3H\left[1-\frac{3H^2-{\gamma }_1{\rho }_1}{\gamma \left(3H^2-{\rho }_1\right)}\right]$
Cosmological constant	$w_{\mathrm{tot}}=-1$	$\Gamma =3H\left[1+\frac{{\gamma }_1{\rho }_1}{\gamma \left(3H^2-{\rho }_1\right)}\right]$

Equations

$$d{s}^2=-d{t}^2+a^2(t)\left[d{r}^2+r^2\left(d{\theta }^2+{\sin }^2\theta d{\varphi }^2\right)\right],$$

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

$$2\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}$$

$$p_c=-\frac{\Gamma (t)}{3H}(p+\rho ),$$

$$\dot{n}+3Hn=n\Gamma \left(t\right),$$

$$nT\dot{s}=\dot{\rho }+3H\left(\left(\rho +p\right)+p_c\right),$$

$$\dot{\rho }+3H\left(1-\frac{\Gamma (t)}{3H}\right)(p+\rho )=0\iff \dot{\rho }+3\gamma H\rho =\Gamma (t)\gamma \rho ,$$

$$\rho =\frac{{\rho }_A}{a^{3\gamma }}\exp \left(\gamma \int \Gamma (t)dt\right)$$

$${\dot{\rho }}_1+3H(p_1+{\rho }_1)=0\iff {\rho }_1={\rho }_{1,0}{a}^{-3{\gamma }_1},$$

$$2\dot{a}\ddot{a}+\left(3\gamma -2\right)\frac{{\dot{a}}^3}{a}-{\rho }_{1,0}(\gamma -{\gamma }_1)a^{1-3{\gamma }_1}\dot{a}-\gamma \Gamma \left({\dot{a}}^2-\frac{{\rho }_{1,0}}{3}{a}^{2-3{\gamma }_1}\right)=0,$$

$$w_{\mathrm{tot}}=-1+\frac{1}{3H^{2}}\Bigg{[}\gamma_{1}\rho_{1}+\gamma\left(1-\frac{\Gamma}{3H}\right)(3H^{2}-\rho_{1})\Bigg{]}=-1+\frac{1}{\rho+\rho_{1}}\Bigg{[}\gamma_{1}\rho_{1}+\gamma\rho\left(1-\frac{\Gamma}{\sqrt{3(\rho+\rho_{1})}}\right)\Bigg{]},$$

$$\Gamma(a=a_{t})=\frac{3H}{\gamma(3H^{2}-\rho_{1})}\Bigl{[}(3\gamma-2)H^{2}+(\gamma_{1}-\gamma)\rho_{1}\Bigr{]}.$$

$$\Gamma \left(H\right)={\Gamma }_0+{\Gamma }_1{H}^{-1}+{\Gamma }_2{H}^{-2}+\sum _{i=3}^n{\Gamma }_i{H}^{-i},$$

$${\Gamma }_j=\frac{1}{j!}{A}^j{\Gamma }_0,j=\mathrm{1...}n\text{with }n\to \infty$$

$$\Gamma \left(H\right)={\Gamma }_0{e}^{A{H}^{-1}}$$

$$a_A\left(\tau \right)=a_{A0}{\tau }^{\frac{2}{3{\gamma }_1}}\text{, and }{a}_B\left(\tau \right)=a_{B0}{\tau }^{\frac{2}{3\gamma }}.$$

$${\rho }_{1,0}=\frac{4a_{A0}^{3{\gamma }_1}}{3{\left({\gamma }_1\right)}^2}\text{.}$$

$$a\left(\tau \right)=a_{\left(A,B\right)}\left(\tau \right)\left(1+\epsilon {\tau }^s\right),$$

$$s_{A1}=-1,s_{A2}=-\frac{2}{{\gamma }_1}\left(\gamma -{\gamma }_1\right).$$

$$s_{B1}=-1,s_{B2}=0.$$

$$s_{C1}=-1,s_{C2}=0.$$

$$s_{A1}=-1,s_{A2}=\frac{1}{2}.$$

$$\frac{a_A\left(\tau \right)}{a_{A0}}={\tau }^{\frac{1}{2}}+a_{A1}\tau +a_{A2}{\tau }^{\frac{3}{2}}+\sum _{i=3}^{\infty }{a}_{Ai}{\tau }^{\frac{1+i}{2}}$$

$$a_{A2}=-\frac{7}{8}{\left(a_{A1}\right)}^2,a_{A3}=\frac{5}{4}{\left(a_{A1}\right)}^3+\frac{3}{5}{a}_{A1}{\Gamma }_0,$$

$$a_{A4}=-\frac{273}{128}{\left(a_{A1}\right)}^4-\frac{13}{10}{\left(a_{A1}\right)}^2{\Gamma }_0,$$

$$a_{A5}=4{\left(a_{A1}\right)}^5+\frac{467}{140}{\left(a_{A1}\right)}^3{\Gamma }_0+\frac{3}{14}{a}_{A1}\left({\left({\Gamma }_0\right)}^2+2{\Gamma }_1\right),$$

$$a_{A6}=-\frac{8151}{1024}{\left(a_{A1}\right)}^6-\frac{2007}{224}{\left(a_{A1}\right)}^4{\Gamma }_0-\frac{3}{1400}{\left(a_{A1}\right)}^2\left(479{\left({\Gamma }_0\right)}^2+615{\Gamma }_1\right),$$

$$a_{A7}=\frac{33}{2}{\left(a_{A1}\right)}^7+\frac{a_{A1}{\Gamma }_0}{720}\left(17775{\left(a_{A1}\right)}^4+3174{\left(a_{A1}\right)}^2{\Gamma }_0+40{\left({\Gamma }_0\right)}^2\right)+\frac{22}{5}{\left(a_{A1}\right)}^3{\Gamma }_1+\frac{a_{A1}}{9}\left(3{\Gamma }_0{\Gamma }_1+4{\Gamma }_2\right).$$

$$\frac{a_B\left(\tau \right)}{a_{B0}}={\tau }^{\frac{1}{2}}+a_{B1}\tau +a_{B2}{\tau }^{\frac{3}{2}}+\sum _{i=3}^{\infty }{a}_{Bi}{\tau }^{\frac{1+i}{2}},$$

$$a_{B1}=\frac{2{\rho }_{1,0}}{{\left(a_{B0}\right)}^3},a_{B2}=-\frac{7}{8}{\left(a_{B1}\right)}^2+\frac{{\Gamma }_0}{6},a_{B3}=\frac{5}{4}{\left(a_{B1}\right)}^3-\frac{7}{15}{a}_{B1}{\Gamma }_0,$$

$$a_{B4}=-\frac{273}{128}{\left(a_{B1}\right)}^4+\frac{31}{240}{\left(a_{B1}\right)}^2{\Gamma }_0+\frac{5}{216}{\left({\Gamma }_0\right)}^2+\frac{1}{9}{\Gamma }_1,$$

Full text

I Introduction

Current astronomical observations show that the recent history of the universe is consistent with accelerated expansion that might be described either by some exotic dark energy fluid in the framework of Einstein’s gravitational theory or by a modification of the gravitational theory itself. This state of affairs has motivated the scientific community to look for alternative theories which can reproduce the effects of the dark energy or modified gravity models in a natural way. The theory of gravitational particle production has a long history. The production of quantum particles by the gravitational field was first investigated in the context of the early universe aniso as a device for damping initial anisotropies, following Schrödinger’s first look at quantum fields in an expanding isotropic universe Sch1939 , although such damping was likely to produce far more radiation entropy than is observed in the microwave background BMatz . Later, a classical counterpart to the quantum particle production was described by Grishchuk grish , and subsequently, cosmological particle production was examined extensively in the literature and used as a means of describing the generation of inhomogeneities during inflation. Many investigators Prigogine-inf ; Abramo:1996ip ; Gunzig:1997tk ; Zimdahl:1999tn ; SSL09 ; LJO10 ; LB10 ; JOBL11 ; LBC12 ; LGPB14 ; CPS14 ; FPP14 ; NP15 ; lima16 ; NP16 ; HP16 ; RCN16 ; PHPS16 (and references therein), have argued that the particle production process might be considered as another approach to explain different phases of the universe’s evolution. In particular, it has been found that such particle creation phenomena can describe early acceleration Prigogine-inf ; Abramo:1996ip ; Gunzig:1997tk ; Zimdahl:1999tn ; RCN16 and late acceleration in the expansion of the universe SSL09 ; LJO10 ; LB10 ; JOBL11 ; LBC12 ; LGPB14 ; CPS14 ; FPP14 ; NP15 ; lima16 ; NP16 ; HP16 ; PHPS16 . Particle creation might provide a way to unify the early- and late- accelerated expansions with intermediate radiation and matter dominated phases of the universe LBC12 ; PHPS16 . The development of a theory of particle production rests upon the choice of the creation rate $\Gamma(t)$, which is an unknown function of time and it can only be determined from the quantum field theory (QFT) in curved spacetime; however, QFT is not yet able to provide the exact functional form for $\Gamma(t)$. Therefore, it is convenient to study some different particle creation rates and build up a picture of the classes of cosmological evolution that arise from different choices for $\Gamma(t)$. This approach is phenomenological but it allows us to constrain the choice of the particle creation rates from their effects. It is possible to identify the particle creation rates with the dark energy equation of state, or the modified gravity effects, because the introduction of a particle creation rate is equivalent to a time-varying (effective) equation of state bar01 . Following this, several phenomenological choices for the form of $\Gamma(t)$ have already made by different authors Prigogine-inf ; Abramo:1996ip ; Gunzig:1997tk ; Zimdahl:1999tn ; RCN16 ; SSL09 ; LJO10 ; LB10 ; JOBL11 ; LBC12 ; LGPB14 ; CPS14 ; FPP14 ; NP15 ; lima16 ; NP16 ; HP16 ; PHPS16 . However, a theory for a general form of $\Gamma(t),$ recovering some well known choices as special cases, is an appealing goal.

It is useful to elaborate on the connection between the quantum particle production studied in references like aniso and the classic picture of particle production that we use here by introducing a non-adiabatic term in the conservation equation for one of our fluids. In effect, this is equivalent to endowing that fluid with a bulk viscosity bulk , which is the only form of dissipative stress allowed in the isotropic models. There are bulk viscous cosmological solutions that display an increasing density of particles for a period, starting from zero at an initial curvature singularity, before their density begins to fall adiabatically because of the domination of the expansion effects. These effects lead to entropy increase so long as the bulk viscous coefficient is positive bulk . Our model described in eq. (6) below is a cosmology with a classical bulk viscous stress that can induce particle creation in the form of particle density increases.

We will consider a spatially flat homogeneous and isotropic universe where the gravitational sector is described by general relativity and the total matter sector is divided into two fluids where one fluid (named as fluid I) is endowed with “adiabatic” particle production process and the second fluid (fluid II) is independently conserved. With this set up, we find that the dynamics of the universe is governed by a nonlinear differential equation which is dependent on the choice of the particle creation rate. We explore the dynamical evolution of this cosmological model using exact solutions of Einstein’s equations for a generalized choice of the particle creation rate, $\Gamma(t)$. By applying the method of singularity analysis to the nonlinear differential equations cotsakis1 ; cotsakis2 ; cotsakis3 ; paperleach ; anleach ; Paliathanasis:2016vsw , we find that their solution can be written in terms of the Laurent series around their movable singularity.

The work has been organized in the following way. In section II, we describe the field equations for a two-fluid system endowed with a particle creation mechanism in the framework of Einstein gravity. In section III we describe the analytic solutions for the cosmological model with a general series form for the particle creation rate. Finally, in section IV we summarise our main findings.

II Field equations for a two-fluid model with particle creation

We consider the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric,

[TABLE]

where $a(t)$ is the expansion scale factor and $t$ is comoving proper time.

We assume that the gravitational sector is described by the Einstein gravity and the matter sector is described by the two-fluid system, with energy densities and pressures: $(\rho,p)$ (fluid I) and $(\rho_{1},p_{1})$ (fluid II) where one fluid (fluid I) is endowed with “adiabatic” particle creation mechanism. With the term “adiabatic” we mean that the entropy per particle remains constant.

The gravitational equations for such two fluid system can be written as (in the units $8\pi G=1$)

[TABLE]

where an overhead dot represents the cosmic time differentiation; $p$, $\rho$ are respectively the thermodynamic pressure and the energy density for fluid I with $p=(\gamma-1)\rho$, ($0<\gamma\leq 2$ is the barotropic state parameter), while $p_{1}$, $\rho_{1}$ represent the same quantities for fluid II with $p_{1}=(\gamma_{1}-1)\rho_{1}$, ($0<\gamma_{1}\leq 2$ is the barotropic index of fluid II). The quantity $p_{c}$ is the creation pressure due to the production of particles and this is related to fluid I by

[TABLE]

where $H=\dot{a}/a$ is the Hubble rate of the FLRW universe and $\Gamma\geq 0$ is the rate of particle creation, which is an unknown quantity. The exact functional form for $\Gamma(t)$ can only be determined from the quantum field theory in curved spacetimes, otherwise it must be modelled by a general functional dependence on physical quantities. Since that subject is not fully developed yet, we assume specific functional forms for $\Gamma$ and try to explain the ensuing cosmological evolution. An important observation about the pressure term (3) is that if fluid I describes the cosmological constant, with $p=-\rho$, then $p_{c}$ vanishes.

We note that the expression (3) follows from the Gibbs equation together with the property that the entropy per particle is constant. In particular, if $n$ stands for the particle number density, then the nonconservation of fluid I particles follows

[TABLE]

while the Gibbs equation, $Tds=d\left(\frac{\rho}{n}\right)+pd\left(\frac{1}{n}\right),$ gives

[TABLE]

in which $p_{c}$ is defined  through equation (3) and ‘$s$’ denotes the entropy per particle which has been considered to be constant. Hence, for the fluid I one finally has

[TABLE]

from which one can solve the evolution for $\rho$ bar01 :

[TABLE]

where $\rho_{A}>0$ is an integration constant. As already noted in bar01 , the evolution equation for $\rho$ in (7) actually provides an equivalent description for a fluid with a dynamical equation of state. We can interpret the first of (6) as the conservation for a fluid with a bulk viscous coefficient equal to $\Gamma(t)\gamma\rho/9H^{2}$, bulk . For different creation rates, different fluids can be realized for a fixed $\gamma$. The second fluid, (fluid II), is assumed not to have any interaction with fluid I, and so obeys the usual conservation equation:

[TABLE]

where the integration constant $\rho_{1,0}$ denotes the present value of $\rho_{1}(t)$.

Using the field equations (1), (2), together with the evolution of the second fluid from (8), we can derive the master differential equation for this two-fluid system, as

[TABLE]

where, as already mentioned, $\gamma,\gamma_{1}>0$. We remark that in the limit $\gamma_{1}=\frac{2}{3}$, the dynamics reduces to the previous work with curvature of the universe added, see bar01 .

We continue with the definitions of the total equation of state $w_{\mathrm{tot}}$ and the deceleration parameter $q(a)$ providing a clear picture of the different phases of the universe. The total equation of state $w_{\mathrm{tot}}$ of the two-fluid system is defined by $w_{\mathrm{tot}}=P_{\mathrm{tot}}/\rho_{\mathrm{tot}}$, where $P_{\mathrm{tot}}=p+p_{c}+p_{1}$, and $\rho_{\mathrm{tot}}=\rho+\rho_{1}$. Using the field equations (1) and (2), this can be recast into the following simplified expression

[TABLE]

where $\rho_{1}$ can be found from (8).

We note from the total equation of state in (10), that we can derive different bounds on it. In other words, the total fluid could behave like a radiation fluid ($w_{\mathrm{tot}}=1/3$), or a dust fluid ($w_{\mathrm{tot}}=0$), or a pure cosmological constant ($w_{\mathrm{tot}}=-1$), depending on the particular different particle creation rate $\Gamma,$ independent of $\gamma$ and $\gamma_{1}$. In Table 1 we have explicitly listed different particle creation rates corresponding to distinct phases of the universe’s evolution.

Following the above equations, we can see that the rate of particle creation, $\Gamma$, plays a key role in determining different stages of the universe’s evolution. In addition, it is interesting to calculate the particle creation rate at the transition scale factor $a=a_{t}$ by solving the equation $q(a=a_{t})=0$, which gives,

[TABLE]

In the next section, we introduce a generalized model for $\Gamma$ and apply the singularity test paperleach ; anleach ; Paliathanasis:2016vsw to find the exact analytical solutions for the master equation (9).

III Singularities and analytic solutions

In this section we present the solutions for the master differential equation (9) by applying the singularity analysis for a general matter creation rate $\Gamma$. We begin our analysis by introducing the following generalized series for the matter creation rate:

[TABLE]

where the $\Gamma_{i}$’s ($i=0,1,...n$) are constants. From equation (11), we see that different choices of the constants, $\Gamma_{i}$, give different matter creation rates. In particular, we can recover some simplest matter creation rates, such as, $\Gamma=\Gamma_{0}$, $1/H$, $1/H^{2}$, as well as quadratic forms, $\Gamma=\Gamma_{0}+\Gamma_{1}/H$, or $\Gamma_{0}+\Gamma_{2}/H^{2}$, and various others. In addition, we see that in the early universe, where $H\rightarrow\infty$, it follows that $\Gamma\left(H\right)\simeq\Gamma_{0}$ from eqn. (11). This is the focus of our investigation because we are looking for singular solutions.

The series (11) can also describe other forms of particle creation function for specific values of the coefficients $\Gamma_{i}$; for instance, when

[TABLE]

then, the series (11) is the Taylor expansion of the exponential function

[TABLE]

around the point at which $H\rightarrow\infty$. Similarly, series (11) can describe other particle creation models for specific values of the coefficients, and our analysis is valid for all the particle creation models which can be described by the series (11), and hence, the present work offers a general study of these particle creation models. The solutions that we derive can also be seen as approximate or exact solutions for all the particle creation models close to the singular solution which corresponds to the dominant term in the series for $\Gamma(H)$.

In order to understand the evolution of the cosmological model given by (9) for the prescribed particle creation rate (11) in Figures 1 and 2, we present the numerical evolution of the deceleration parameter $q\left(a\right)$ for the initial conditions $a\left(t_{0}\right)=1$ and $\dot{a}\left(t_{0}\right)=70$, where the present value of the Hubble constant is $H_{0}=\frac{\dot{a}\left(t_{0}\right)}{a\left(t_{0}\right)}$.

The numerical solution in Fig. 1 is for $\gamma=1$ and $\gamma_{1}=\frac{4}{3}\,,$ and we select $\rho_{1,0}=3\Omega_{r0}H_{0}^{2}$, with $\Omega_{r0}\simeq 5\times 10^{-3}$.  On the other hand, Fig. 2 is for $\gamma=\frac{4}{3}$ and $\gamma_{1}=1$ with $\rho_{1,0}=3\Omega_{m0}H_{0}^{2}$ and $\Omega_{m0}\simeq 0.3$. In both figures the $\Gamma_{i}$ coefficients of the particle creation function $\Gamma\left(H\right)$ have been considered to be positive.

From Fig. 1 and Fig. 2, we observe that in the early universe, the universe is dominated by the matter source and there always exists a de Sitter point as a future attractor. These de Sitter points can be easily calculated in a similar way as in bar01 .

III.1 ARS algorithm

In order to derive the analytic solutions of the master equation (9) we work with the approach of Kowalevskaya Kowa , and determine if the second-order differential equation (9) possesses the movable singularities. The existence of a movable singularity means that the solution of equation (9) near the singularity is described by the power-law function $a\left(\tau\right)\simeq\tau^{p},~{}\tau=t-t_{0},$ where $p$ is a negative number and $t_{0}$ denotes the position of the singularity. The position of the singularity changes according to the initial condition of the problem which means that different initial conditions provide us with different locations for the singular point.

Ablowitz, Ramani and Segur (ARS) Abl1 ; Abl2 ; Abl3 , proposed an algorithm to test if a given differential equation is integrable when an algebraic solution exists given by a Laurent expansion.

The ARS algorithm is briefly described by the following three main steps paperleach :

(a) Determine the leading-order behaviour, at least in terms of the dominant exponent. The coefficient of the leading-order term may or may not be explicit.

(b) Determine the exponents at which the arbitrary constants of integration enter.

(c) Substitute an expansion up to the maximum resonance into the full equation to check for consistency.

For the singularity analysis the exponents of the leading-order term needs to be negative integers or a non-integral rational number. However, this is not so restrictive since we can always perform a change of coordinates in order to obtain a negative exponent for the leading-order term.

In order that the differential equation passes the Painlevé test, the resonances have to be rational numbers, so that the solution can be written as a Painlevé series. Alternatively, if at least one of the resonances is an irrational number, the differential equation passes the weak Painlevé test. Moreover, the value minus one $\left(-1\right)$ should always appear as one of the resonances. The existence of that resonance is important in order for the singularity to be movable. For a right Painlevé series (right Laurent expansion) the resonances must be nonnegative; for a left Painlevé series (left Laurent expansion) the resonances must be non-positive while for a full Laurent expansion the resonances have to be mixed. Obviously, the possible Laurent expansions for second-order differential equations are either left or right Painlevé series. For a review and various applications we refer the reader to Ref. buntis , while different applications of singularity analysis in cosmological studies can be found in Paliathanasis:2016vsw ; sin1 ; sin2 ; sin3 ; sin3b ; sin4 ; sin5 and references therein.

III.2 Leading-order behaviour

We continue by applying the first step of the ARS algorithm to determine the leading-order behaviour. We substitute, $a\left(t\right)=a_{0}\left(t-t_{0}\right)^{p}$, in the master equation (9), where we find that for $\gamma\neq\gamma_{1}$, there are two possible leading-order behaviours are described by the power-law functions111The leading-behaviour $a_{B}\left(t\right)$ exists only when $\gamma\geq\gamma_{1}.$

[TABLE]

The scale factors $a_{B}\left(\tau\right)$ describe solutions with perfect fluids. The function $a_{A}\left(\tau\right)$ corresponds to the matter solutions in which the perfect fluid with $\gamma_{1}$ dominates, while function a_{B}\left(\tau\right)\describes a solution in which the fluid term with $\gamma$ dominates.

Furthermore, the coefficient parameter $a_{B0}$ is found to be arbitrary, while $a_{A0}$ is related to the energy density $\rho_{1,0}$ by

[TABLE]

The coefficient $a_{B0}$ is the second-integration constant which controls the generic solution of the master equation (9). The position in the Laurent expansion of the second integration constant in the leading-order behaviour, $a_{A0}$, is determined by the values of the resonances. The later will be found below. Finally, note that in the case where $\gamma=\gamma_{1}$, there exists only one leading-order term, the $a_{B}\left(t\right),$ with an arbitrary value for the coefficient $a_{B0}$.

III.3 Resonances

The second step in the ARS algorithm is the determination of the resonances. We substitute the following expression,

[TABLE]

in the master equation (9) then we linearize around $\varepsilon=0.$ The coefficients of the leading-order behaviour provide a polynomial equation of order two, whose solutions are the resonances, $s$, of the leading-order behaviours $a_{A}\left(\tau\right)$ and $a_{B}\left(\tau\right)$.

For the dominant term, $a_{A}\left(t\right),$ the resonances are calculated to be

[TABLE]

This means that, for $\gamma>\gamma_{1}$, the solution will be given by a left Painlevé series, while when $\gamma<\gamma_{1}$ the solution is given by a right Painlevé series.

For the dominant term, $a_{B}\left(t\right)$, the resonances are calculated to be

[TABLE]

Resonance $s_{B2}$ indicates that the second integration constant is $a_{B0}$ – which we calculated above.

In the case where $\gamma=\gamma_{1}$, the two resonances are found to be

[TABLE]

From the values of these resonances we can extract more information than the position of the second integration constant. In particular, from the discussion in anleach we find that when the solution is given by a left Painlevé series, the leading-order term describes an attractor solution, while when the solution is given by a right Painlevé series the leading-order behaviour is described by a source point and so is an unstable solution. Furthermore, when the solution is given by a full Laurent expansion the leading-order behaviour is described by a saddle point.

III.4 Consistency tests and analytic solutions

We continue with the consistency test and we write the analytic solution in terms of Laurent expansions. However, in order to continue it is necessary to select values of the equation of state parameters of the two perfect fluids.

The consistency tests that we present are for combinations of dust and radiation:

(A) $\left(\gamma,\gamma_{1}\right)=\left(1,\frac{4}{3}\right)$, (B) $\left(\gamma,\gamma_{1}\right)=\left(\frac{4}{3},1\right)$, (C)$~{}\gamma=\gamma_{1}$ with $\gamma=\frac{4}{3},~{}$and (D) $\gamma=\gamma_{1}$ with $\gamma=1$.

The coefficient terms that we calculate are the first terms until the parameter $\Gamma_{2}$ of (11) appears.

III.4.1 Dust with particle creation and radiation

When $\left(\gamma,\gamma_{1}\right)=\left(1,\frac{4}{3}\right)$, the resonances which correspond to the leading-order term $a_{A}\left(\tau\right)=a_{A0}\tau^{\frac{1}{2}},$ take the values

[TABLE]

This means that the solution starts from the radiation era and is given by a right Painlevé series. Moreover, we can infer that the solution in the radiation era is a source, that is, it is an unstable solution.

The analytic solution is given by the right Painlevé series with step $\frac{1}{2}$,

[TABLE]

in which $\rho_{1,0}=\frac{4}{3}\left(a_{A0}\right)^{3}$.

The second integration constant of the solution is the coefficient term $a_{A1}$, while the rest of the coefficients, $a_{Ai}$, are functions of $a_{A1}$ and the parameters $\Gamma_{i}$, as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

From these expressions it is clear that the higher polynomial terms of the particle creation function (11) play a crucial role in the analytic solution as we evolve far from the radiation era.

III.4.2 Radiation with Particle creation and pressureless fluid

In order to test the consistency of the second solution with leading-order term $a_{B}\left(\tau\right)$, we select $\gamma=\frac{4}{3}$ and $\gamma_{1}=1$.

Hence, the analytic solution is given by the right Painlevé series with step $\frac{1}{2}$,

[TABLE]

where now $a_{B0}$ is the second integration constant. The first six coefficients are given by the following expressions

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Again, we observe that as we go far from the leading-order term, i.e. from the radiation era, the coefficients of the higher polynomial terms are involved in the solution.

III.4.3 Radiation-like fluid with particle creation and radiation

Now, in the case in which $\gamma=\gamma_{1}$ and $\gamma=\frac{4}{3}$, the generic analytic solution is given by the Laurent expansion

[TABLE]

where $a_{C0}$ is arbitrary and the first non-zero coefficients are

[TABLE]

[TABLE]

We note that the non-zero coefficients are the $a_{Ci}=a_{C\left(2k\right)}$, $~{}k\in\mathbb{N}.$

III.4.4 Dust with particle creation and pressureless fluid

In a similar way, when $\gamma=\gamma_{1}$ and $\gamma=1$, the general analytic solution is expressed by the Laurent expansion

[TABLE]

which starts from the matter-dominated era.

The coefficient $a_{D0}$ is the second integration constant of the solution, while the first non-zero coefficients are calculated to be,

[TABLE]

[TABLE]

We note that the non-zero coefficients are the $a_{Di}=a_{D\left(3k\right)}$, $~{}k\in\mathbb{N}.$

IV Conclusions

* *An explicit form of dark energy or the presence of modifications to general relativity gravity are two independent roads towards an explanation for the observed acceleration of the universe. Here we consider a third alternative through which we can describe the observational results in a convenient way. A theory of gravitational particle production has shown that different phases of the universe can be explained Prigogine-inf ; Abramo:1996ip ; Gunzig:1997tk ; Zimdahl:1999tn ; SSL09 ; LJO10 ; LB10 ; JOBL11 ; LBC12 ; LGPB14 ; CPS14 ; FPP14 ; NP15 ; lima16 ; NP16 ; HP16 ; RCN16 ; PHPS16 and so it was argued that the particle creation formalism might be considered as a viable alternative to the dark energy or modified gravitational theories. The key role in gravitational particle production is played by the rate of particle creation $\Gamma(t)$ which has an equivalent character to an effective dynamical equation of state, and so we can consider various phenomenological models for $\Gamma(t)$. But, except for some simple choices of the creation rate, the dynamics of the universe cannot be obtained in an analytic way. This motivated the present paper, where for the first time we present the exact solutions of the gravitational field equations for a system of two fluids model including particle creation.

In particular, assuming a spatially flat Friedmann-Lemaître-Robertson-Walker universe described by the Einstein gravity, we consider two perfect fluids with barotropic equations of state where one of them is endowed with particle creation, as if it possessed a bulk viscosity, while the other fluid obeys the standard conservation law. Interestingly, we found that the dynamics of such a two-fluid particle creation system can be concisely described by a single nonlinear differential equation (9) that involves the particle creation rate $\Gamma(t)$. We recall that the above two-fluid particle creation system might be considered to be an equivalent cosmological scenario to a single fluid associated with particle creation in the presence of curvature bar01 . However, since the particle creation rate $\Gamma(t)$ is any unknown function and its detailed functional form is unknown (it must be derived from QFT in future), we widen our investigations by allowing a general particle creation rate (11) that recovers some specific particle creation models as special cases, namely, $\Gamma=\Gamma_{0}=$ constant, $\Gamma\propto 1/H^{n}$ ($n\in\mathbb{N}$), $\Gamma=span\{\Gamma_{0},H^{-1},H^{-2},...,H^{-n}\}$ ($n\in\mathbb{N}$), as well as some other exceptional but interesting choices like $\Gamma\propto\exp(1/H)$. Following a singularity analysis applied to the nonlinear differential equation defining the theory, we see that the master equation in this work, i.e., equation (9) can pass the singularity test and hence, the solution for the gravitational equations can be written in terms of the Laurent series around the movable singularity. We note that we do not consider choices like $\Gamma\propto H^{2}$ in the general model of $\Gamma$ in (11), since for such models, the governing differential equation does not pass the singularity test bar01 .

Therefore we see that for the present two-fluid cosmological model with gravitationally induced “adiabatic” particle creation, we can obtain the singular algebraic solutions to the gravitational field equations for a class of particle creation models given in (9). Finally, we mention that the present work can be extended to more than two fluids model in presence of the particle creation, although the dynamics could be more complicated .

Acknowledgements.

The authors thank the referee for some specific comments on the article. SP acknowledges the research grant under Faculty Research and Professional Development Fund (FRPDF) Scheme of Presidency University, Kolkata, India. JDB was supported by the Science and Technology Facilities Council (STFC) of the United Kingdom. AP acknowledges the financial support  of FONDECYT grant no. 3160121.