Dynamical system analysis of FLRW models with Modified Chaplygin gas

We investigate Friedmann–Lamaitre–Robertson–Walker (FLRW) models with modified Chaplygin gas and cosmological constant, using dynamical system methods. We assume \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=(\gamma -1)\mu -\dfrac{A}{\mu ^\alpha }$$\end{document}p=(γ-1)μ-Aμα as equation of state where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}μ is the matter-energy density, p is the pressure, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α is a parameter which can take on values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<\alpha \le 1$$\end{document}0<α≤1 as well as A and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document}γ are positive constants. We draw the state spaces and analyze the nature of the singularity at the beginning, as well as the fate of the universe in the far future. In particular, we address the question whether there is a solution which is stable for all the cases.

Scientific Reports | (2021) 11:2750 | https://doi.org/10.1038/s41598-020-80396-w www.nature.com/scientificreports/ and modified Chaplygin gas as well as modified Chaplygin gas will behave similar to CDM in the far future. Li et al. 24 studied the modified Chaplygin gas model that interacts with cold dark matter. Bhadra and Debnath 25 studied an interacting modified Chaplygin gas model, using dynamical system methods. Fang et al. 26 considered a potential parameter Ŵ( ) and made a three dimensional analysis of quintessence scalar field.
In this study, we use modified Chaplygin gas with p = (γ − 1)µ − A µ α as equation of state where γ is a parameter which is in the interval 0 ≤ γ ≤ 2 , µ is the matter-energy density, p is the pressure, α is a parameter which can take on values 0 < α ≤ 1 as well as A is a positive constant.It is clear that when γ = 1 , generalized Chaplygin gas is recovered. If α = 1 as well, it reduces to pure Chaplygin gas while A = 0 recovers the equation of state of perfect fluid 15,27,28 . We also include the nonnegative cosmological constant in our study. Therefore, we obtain a three dimensional state space where a mixture of perfect fluid, generalized Chaplygin gas and cosmological constant exist.
The outline of the paper is as following. In section "Methods", we show how the dynamical system is obtained. In "Results" section, using Friedmann equation and weak energy condition we obtain the constraints of the state spaces; find the equilibrium points where all the derivatives are zero; obtain the eigenvalues using the linearized Jacobian matrix as well as construct the three dimensional complete state space. We end with some conclusions in "Discussion" section.

Methods
FLRW models are spatially homogeneous and isotropic models and represented by the following metric in a (t, r, θ, ψ) co-moving coordinate system. This metric which is called Robertson-Walker (RW) metric is independent of field equations. f k (r) function, which depends on the parameter k, is defined as following. According to this; universe is elliptical if k = +1 , Euclidean if k = 0 and hyperbolic if k = −1 . t and l(t) are called cosmic time and distance scale factor, respectively.
Defining Hubble scalar by H =l/l , the field equations and equations of motion become as following 29 Here Eqs. (3)(4)(5) are the Friedmann equation, the Raychaudhuri equation and the energy conservation equation, respectively where 3 R is the three-curvature of the symmetry surfaces, q is the deceleration parameter, µ is the energy density, p is the pressure and is the cosmological constant as well as a dot denotes the derivative with respect to cosmic time.
We use Eq. 6 as equation of state.
In order to make the variables dimensionless, it is necessary to norm them to a suitable variable. This suitable variable could be determined with the help of Eq. (3).
Matter-energy density µ is nonnegative because of the weak energy condition and 3 R ≤ 0 if k = 0 or k = −1 . If is nonnegative (� ≥ 0) as well, the term H 2 becomes greater than the other terms so it becomes the dominant term.
Therefore, the variables could be made dimensionless by dividing them by the powers of H if H = 0 . According to this, for k = 0, −1 and ≥ 0 , the following dimensionless variables can be defined where , , K, A and P are matter density, density, curvature, Chaplygin gas density and pressure parameters, respectively which are dimensionless.
Considering these dimensionless variables as well as defining a dimensionless time variable as dτ = H(t)dt , we obtain following evolution and constraint equations as well as the definition equation of q. (1) (2) f k (r) = sin r k = +1 r k= 0 sinh r k = −1 If we introduce a time variable as dτ = D(t)dt , we obtain following evolution and constraint equations where ′ denotes the derivative with respect to τ.

Results
For k = 0, −1 case, since H ′ equation decouples, we obtain a three-dimensional reduced dynamical system for , and A which are given in Eqs. (9)(10)(11). Considering the weak energy condition ( + P ≥ 0 ) and ≥ 0 inequality as well as constraint equation 13, we get the following conditions. Thus, the state space is three-dimensional and compact. Jacobian matrix becomes as following for this case: For k = 0, +1 case, D ′ equation decouples and we get a three-dimensional reduced dynamical system for Q, and A which are given in Eqs. (16)(17)(18). Considering the weak energy condition ( + P ≥ 0 ) and ≥ 0 inequality as well as constraint equation 20, we get the following conditions.: This state space is three-dimensional and compact as well. Jacobian matrix becomes as following for this case: Hence, systems of equations for A = 0 submanifold become as in Eqs. (26)(27)(28)(29)(30) and Eqs. (31-35) for k = 0, −1 and k = 0, +1 cases, respectively. These give the same state spaces that are constructed by Goliath and Ellis for perfect fluid matter source 14 .
State space for k = 0, +1 includes both contracting and expanding cases while state space for k = 0, −1 depends on the sign of H. Combining one contracting and one expanding state spaces of k = 0, −1 case with the state space of k = 0, +1 case on k = 0 submanifold, we obtain the full dynamical system which is three dimensional and seen in Figs. 1 and 2. z axis corresponds to A for all of the state space. On k = 0, +1 region (middle) vertical and horizontal axes correspond to and Q, respectively. On k = 0, −1 region (triangular parts on both sides), M-F and M-dS lines correspond to and axes, respectively. The limits of the state space determined by Eqs. (21 and 23) for k = 0, −1 and k = 0, +1 cases, respectively.
As seen in Figs. 1 and 2, the dynamical system is bounded by A = 0 (bottom), = 0 (back) and � A = α+1 √ γ � (front) invariant submanifolds. k = 0 invariant submanifold lies between the intersection of k = 0, +1 and k = 0, −1 state spaces. The state space has a different structure depending on the value of γ . Since the value of α only changes the positions of Einstein static solutions but does not change the behaviour of the system, we only give α = 1 case.
The dynamical system has a number of solutions which are given in Tables 1 and 2 where the eigenvalues are determined by the jacobian matrices (22) and (24). Except from E(Einstein static solution curve) which exists for γ >  Table 3.
The state space for γ > 2 3 is seen in Figs. 1, 3 and 5. For k = 0, −1 and H > 0 , future attractor is the CD equilibrium line. The universes, where and A are both equal to zero, first evolve to Milne model from F(flat Friedmann solution) before expanding to a CD universe. The other universes on this part of the state space expand from F and evolve directly to a CD model. The time reverse region occurs for H < 0.
For k = 1 and H > 0 , models which are past asymptotic to F and start with sufficiently small (� � + � A ) begin contracting after some expansion and collapse to a big crunch at F while the future attractor for the models with large enough (� � + � A ) is a CD universe. For H < 0 , models that start with large enough (� � + � A ) cross to H > 0 region and expand back to their initial conditions while universes with sufficiently small (� � + � A ) are future asymptotic to F. There are universes which are past asymptotic as well as future asymptotic to Einstein static model as well. For k = 0, −1 and H > 0 , Milne universe is the source and the universes that are past asymptotic to it, first expand to F and then evolve to a CD universe, if A and are both zero. Otherwise, they expand directly to a CD model. When H < 0 , time reverse region occurs.
For k = +1 , universes are past asymptotic to contracting versions of the equilibriums as well as future asymptotic to the expanding versions.
, respectively. z axis corresponds to A for all of the state space. On k = +1 region (middle) vertical and horizontal axes correspond to and Q, respectively. On k = −1 region (triangular parts on both sides), M-F and M-dS lines correspond to and axes, respectively. The right and left parts of the state space represent expanding and contracting universes, respectively. Subscripts on equilibriums refer to the sign of H there. Lines go from red to black. (MATLAB ver. R2019b).

Discussion
For γ > 2 3 case, past attractor for contracting universes is pressureless flat models (CD solution line). From this solution line, if contraction starts with enough(� � + � A ) and positive curvature, universe evolve to an expanding CD model. Otherwise, it contracts to matter dominated flat model (Flat Friedmann universe) either directly or first evolving to Milne or Einstein static universes. From this matter dominated flat model, it may bounce back and evolve to contracting flat Friedmann model again or to a CD model depending on the initial condition of the disturbance. Thus, if the universe starts with small (� � + � A ) and positive curvature as well as the configuration of the universe does not change at the singularity, universe may bounce back and forth between contracting and expanding flat Friedmann universes. Otherwise, it bounces back and forth between pressureless flat universes (CD solution line) either directly or first evolving to other models.
For γ < 2 3 case, universe expand to a pressureless flat model from either flat Friedmann or Milne universes, if it starts with positive Hubble parameter. Since there are no direction on this solution line, universe may evolve to any model between CH and dS models as long as Hubble parameter is positive. If, it bounces back from here and begin contracting, there are three possible directions for it. If it disturbed to positive curvature, contraction slows down and the universe begin expanding without a singularity and evolve back to expanding CD model. If it starts with zero curvature, universe contracts to flat Friedmann model and from this singularity it may evolve to an expanding flat Friedmann universe or bounce back to an expanding CD model or contracts to Milne universe.  www.nature.com/scientificreports/ From Milne singularity, it expands to a CD universe. Thus, we have a cycle starting with as well as finishing at a pressureless flat model. For both cases, we have a dynamical universe and there is no final state for it as long as we don't put any restriction on H. We have examined the evolution of FLRW models with modified Chaplygin gas where the equation of state is p = (γ − 1)µ − A µ α . It is seen that value of α has negligible effect on the structure of the state space as well as no effect on the stability of the equilibriums. For γ < 2 3 , there are 2 equilibrium points which are F and M as well as CD equilibrium line where � A = α+1 √ γ � . When γ > 2 3 , there is an additional equilibrium curve which corresponds to Einstein static universe. Unless we put a restriction on the value of H, there is no stable solution and universe may evolve from one equilibrium point to another continuously but if we assume that Hubble constant is positive, it evolves to a pressureless flat model. We are planning to apply this method to the Bianchi models in future papers.