Dynamical systems of cosmological models for different possibilities of G and ρΛ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\Lambda }$$\end{document}

The present paper deals with the dynamics of spatially flat Friedmann–Lemaître–Robertson–Walker (FLRW)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(FLRW)$$\end{document} cosmological model with a time varying cosmological constant Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document} where Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document} evolves with the cosmic time t through the Hubble parameter H, that is, Λ(H)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda (H)$$\end{document}. We use the expression of Λ(H)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda (H)$$\end{document} in the form of Taylor series with respect to H keeping only the even powers of H because of the general covariance of the effective action of quantum field theory in a curved background. Dynamical systems for three different cases based on the possibilities of gravitational constant G and the vacuum energy density ρΛ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\Lambda }$$\end{document} have been analysed. In Case I, both G and ρΛ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\Lambda }$$\end{document} are taken to be constant. We analyse stability of the system by using the notion of spectral radius, behavior of perturbation along each of the axes with respect to cosmic time and Poincaré sphere. In Case II, we have dynamical system analysis for G=constant\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G=\hbox {constant}$$\end{document} and ρΛ≠\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\Lambda } \ne $$\end{document} constant where we study stability by using the concept of spectral radius and perturbation function. In Case III, we take G≠\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G \ne $$\end{document} constant and ρΛ≠\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\Lambda } \ne $$\end{document} constant where we introduce a new set of variables to set up the corresponding dynamical system. We find out the fixed points of the system and analyse the stability from different directions: by analysing behaviour of the perturbation along each of the axes, Center Manifold Theory and stability at infinity using Poincaré sphere respectively. Phase plots and perturbation plots have been presented. We deeply study the cosmological scenario with respect to the fixed points obtained and analyse the late time behavior of the Universe. The effective equation of state parameter ωeff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{eff}$$\end{document}, total energy density Ωtt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{tt}$$\end{document} are also evaluated at the fixed points for each of the three cases and these values are in agreement with the observational values in Aghanim et al. (Astron Astrophys 641(A6): 2020, 2018). We have also presented the EoS parameter for dark energy sector ωde(zr)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{de}(z_{r})$$\end{document}, the Hubble parameter H(zr)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H(z_{r})$$\end{document} and the deceleration parameter q(zr)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(z_{r})$$\end{document} as functions of redshift zr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{r}$$\end{document} for all the three cases and their plots over redshift are also provided. We analyse the quintessence-like, phantom-like or the purely cosmological-constant type dark energy, etc behavior when the EoS approaches the fixed point value near -1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-1$$\end{document}. The present values of ωde(zr)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{de}(z_{r})$$\end{document}, H(zr)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H(z_{r})$$\end{document}, q(zr)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(z_{r})$$\end{document} and zrt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_{rt}$$\end{document} have been tabulated in Table 4 and they fall within the range of cosmological observations. The transition redshift value (zrt)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(z_{rt})$$\end{document} for each of the three cases have also been evaluated. In each of the cases the developed model agrees with the fact that the Universe is in the epoch of accelerated expansion for suitable values of free parameters chosen. The developed cosmological models associated with each of the three cases have a deep connection with the accelerated expansion phenomena of the evolving Universe.


Introduction
In the past two decades many researchers have put tremendous efforts to develop and improve the plethora of theoretical models that explain the accelerated expansion of our Universe. Astrophysical measurements that reveal such a phenomenon put into the quest to give convincing theoretical explanations from various possible directions [1][2][3][4][5][6][7][8][9][10][11]. The dark energy model is one such proposed model that attributes the expansion phenomenon to an energy component with negative pressure so called dark energy which dominates the universe at late time. The simplest type of dark energy is the cosmological constant [12]. In this context of accelerated expansion the theory of general relativity (G R) modified by a cosmological constant term , which is known as the famous CDM model is one of the most popular one [13]. But, despite its fine agreement with the observation data, there are two major issues that have driven our young minds to focus sharply on some modifications to the assumed CDM model, namely, "the cosmological constant problem" which deals with the discrepancy between theoretical and expected values of the cosmological constant [14][15][16]; and "the cosmic coincidence problem" [17]. To mend up these issues, running cosmological models have been developed.
In the year 2002, Shapiro and Solá made the first development regarding the scaling evolution of the cosmological constant [18]. Various possible cosmological and astrophysical implications have been presented on variable cosmological constant as Planck scale effect and on running G as well as at low energies from the side of physics [19,20]. Among the running cosmological constant models that have been proposed, it is worthy enough to mention the time dependent cosmological constant motivated by quantum field theory [21][22][23], (t) cosmology induced by a slowly varying Elko field [24], a running vacuum in the context of supergravity [25], etc. The first paper where an action functional is associated to running vacuum was presented in [26]. The models of the type (H ) have important implications both on the theoretical arena as well as phenomenologically to smooth out some of the existing tensions of the C DM model. The most important implication on the theoretical side is that they can help to solve the cosmological constant problem [27,28]. Some of the related papers which can be of interest on these matters which are worthy of mentioning are [29][30][31][32]. A detailed theoretical and phenomenological explanation about these models can be found in [33,34] also. In Newtonian gravity, without any requirement of further constraints to be satisfied we can explicitly write the time variation of G. But in GR there are other constraints to be satisfied. For instance if we assume that the ordinary energymomentum conservation law holds then there should not be any variation in the gravitational coupling with respect to the space time or otherwise the ordinary energy-momentum conservation law will be violated [35,36]. A direct application of these ideas on running vacuum energy to study the variation of the fundamental constants of nature is given in [37,38]. In the light of Dirac's idea [39][40][41] which propose that some of the fundamental constants cannot remain constant forever, it is essential to do some modifications in GR field equations [42][43][44] if we are to consider this running cosmological constant term. In this regard, studying the cosmic scenario with varying G needs modified field equations as well as modified conservation laws. We can mention Brans-Dicke theory where there are modifications of GR with a varying G without violating the ordinary energy-momentum conservation law [45][46][47]. There are many other models that employ varying G theories that give a better understanding of the Universe regarding its late time behavior and nature [47][48][49][50][51][52][53][54][55][56][57][58][59][60]. As there are no rigorous proves that indicate whether the cosmological constant is running or not [60], one can study the cosmological implications of different possible theoretical assumptions of term. Motivated by the quantum field theory [19,20,61] and some theoretical motivations about the varying form [21,22], Aleksander Stachowski and Marek Szydtowski have also studied the dynamics of cosmological models with various forms of (t) [62].
In this paper, we consider a running vacuum model which evolves in power series of H . Our aim is to set up dynamical systems out of the cosmological field equations by introducing new set of variables and study the stability of the systems in the light of cosmological implications of the system. Based on the possibilities of the gravitational constant G and the vacuum energy density ρ , we develop different dynamical system for three cases and analyse the stability through different approaches by finding respective fixed points. The cosmological scenario associated with each fixed point has been discussed in detail. We arrange the paper in the following ways. In Sect. 2, we give preliminaries that provides a brief introduction on dynamical systems approach to cosmology with some definitions and theorems which will be required to understand the subsequent analysis in the paper. In Sect. 3, we have three cases. In case I of Sect. 3 we show the setting up of cosmological equations and dynamical system analysis where both G and ρ are taken to be constants which is the case of standard CDM cosmology. Under Case I we have three subsections based on analysis using spectral radius, perturbation function and stability at infinity using Poincaré sphere. We present, in Case II, the model dynamics where G = constant and ρ = constant. Under case II, we have two subsections based on analysis through spectral radius and using perturbation along each of the axes with respect to increase in cosmic time. In Case III we have dynamical system analysis where G = constant and ρ = constant. Under Case III we present three subsections on the basis of analysing stability by the use of perturbation function, Center Manifold Theory and Poincaré sphere. In Sect. 4 we give conclusion of our study. Stability analysis for each of the cases at the respective fixed points is presented and their corresponding cosmological implications along with the evaluation of various cosmological parameters at the respective fixed points are also obtained.

Preliminaries
Dynamical system is a mathematical system that describes the time dependence of the position of a point in the space that surrounds it, termed as ambient space. Here, we are approaching towards the system through an autonomous system of ordinary differential equations, (ASODE). ASODE is a system of ordinary differential equations which does not depend explicitly on time. Surendra and Chingtham [63] have also used this approach to study cosmological models in the presence of a scalar field using different forms of potential. From [63] we can also notice that in three dimensional dynamical system we can analyse stability by analysing the nature of perturbation along each of the axes. A dynamical system is generally written in the form of the following [64]: is an element of state space X ⊆ R n and the function f : X → X . The overhead dot denotes the derivative with respect to cosmic time, t.
. . , n denotes the first partial derivative of f i with respect to the i th component x i of the element The concept of Linear stability theory is one of the simplest ways used to understand the dynamics of a system near a fixed point. In Linear stability theory, the function f is assumed to be sufficiently regular so that we can linearise the system around its fixed point. The eigenvalues of the Jacobian matrix at a fixed point play an important role in studying the stability of the fixed point.
For hyperbolic fixed points if all the eigenvalues of J x o have positive real parts, then x o acts as a repeller and it is unstable as all the trajectories closed enough to it are repelled from it. x o is stable when all the eigenvalues of J x o have negative real parts. Here x o is called as attractor and it attracts all nearby trajectories towards it. If at least two eigenvalues have real parts with opposite sign then, x o behaves as a saddle fixed point which attracts trajectories in some directions and repels along other directions.
If at least one of the eigenvalues of the Jacobian matrix at a fixed point x o have zero real part then we can not do stability analysis by using eigenvalues of the Jacobian matrix. Such a fixed point is referred to as non-hyperbolic fixed point. To analyse stability of such fixed points we need a better approach other than the linear stability analysis like Center manifold theory, perturbation function, Lyapunov stability. Center manifold theory is the most popular method which reduces the dimensionality of the system and determines the stability of the critical points of the parent system according as the stability of the reduced system. Wiggins [64] and Carr [67] have discussed the Center manifold theory in detail.
The eigenvalues of the Jacobian matrix J with order n × n given in Definition 2.5 will have n eigenvalues. The eigenvectors of J associated to the eigenvalues with negative real part spans a vector space called stable space, J s and the eigenvectors associated with positive real part spans a vector space called the unstable space, J u . Similarly J c represents the vector space spanned by the eigenvectors associated with zero real part. Here, the superscript s, u, c denote the dimensions of the respective vector spaces. Also the spaces J s , J u and J c are the subspaces of R n . The space R n can be written as the direct sum of these three subspaces, that is, R n = J s ⊕ J u ⊕ J c . These results have been detailed in Carr's book [67], Elaydi's book [68] and Zhang's book [69]. If at least one eigenvalue of J at a fixed point x o has positive real part then x o will be unstable whether it is hyperbolic or not. But if x o is non-hyperbolic and no eigenvalues has positive real part, then we can use Center manifold theory to determine stability of the fixed point.
Let us consider a two dimensional dynamical system. Using a suitable coordinate transformation we can rewrite any system of the form (1) as follows: where A is a c × c matrix having eigenvalues with zero real parts, B is an s × s matrix having eigenvalues with negative real parts and (x,y) ∈ J c × J s . The functions f and g satisfy the following: for a sufficiently regular function h(x) on J s and δ however small it may be. The proofs of the existence of the center manifold for the system (2) is also provided in [67] and he has given the dynamics of the system (2) restricted to the center manifold as follows: for sufficiently small v ∈ R c .
Theorem 2.1 [70] Consider a flow defined by a dynamical system on R 2 where P 1 and P 2 are polynomial functions of x and y. Let P 1m and P 2m denote the m th degree term in P 1 and P 2 respectively. Then, the critical points at infinity for the m th degree polynomial system (9) occur at the points (X, Y, 0) on the equator of the Poincaré sphere where or equivalently at the polar angle θ j and θ j + π satisfying This equation has at most m + 1 pairs of roots θ j and θ j + π unless G m+1 (θ ) is identically zero. If G m+1 (θ ) is not identically zero, then the flow on the equator of the Poincaré sphere is counter-clockwise at points corresponding to polar angles θ where G m+1 (θ ) > 0 and it is clockwise at points corresponding to polar angles θ where G m+1 (θ ) < 0.
Theorem 2.2 [70] The flow defined by (9) in a neighborhood of any critical point of (9) on the equator of S 2 , except the points (0, ±1, 0), is topologically equivalent to the flow defined by the following system the signs being determined by the flow on the equator of S 2 as determined in Theorem 2.1.

Theorem 2.3 [70]
Let us consider a flow in R 3 defined bẏ where P 1 , P 2 and P 3 are polynomial functions of x, y, z of maximum degree m. The critical points at infinity for the m th degree polynomial system (10) occur at the points (X, Y, Z , 0) on the equator of the Poincaré sphere S 3 where X 2 + Y 2 + Z 2 = 1 and where P 1m , P 2m and P 3m denote the m th degree terms in P 1 , P 2 and P 3 respectively. Theorem 2.4 [70] The flow defined by the system (10) in a neighborhood of (±1, 0, 0, 0) ∈ S 3 is topologically equivalent to the flow defined by the system:

Dynamical system analysis for different possibilities of G and ρ
In this section we present the dynamical system analysis when G =constant and ρ =constant. This is a standard model and we present it as case I of our analysis.

Case I: Dynamical system analysis when G = constant and ρ = constant
The Einstein field equations in the presence of cosmological constant are given by where T μν is the ordinary energy-momentum tensor, T μν ≡ T μν + g μν ρ is the modified energy-momentum tensor and ρ = 8π G is the vacuum energy density in the presence of .
We assume that the universe is filled with a perfect fluid with velocity four-vector field V μ . With this consideration we have T μν = −p m g μν + (ρ m + p m )U μ U ν , where ρ m is the density of matter-radiation and p m = (γ − 1)ρ m is the corresponding pressure. In the similar way, the modified energy-momentum tensor can be expressed as where p tt = p m + p , ρ tt = ρ m + ρ and p = ω de ρ is the associated pressure in the presence of . With this substitution in the above expression we have By assuming a spatially flat Friedmann-Lemaître-Robertson-Walker(FLRW) metric along with the above modified energy-momentum tensor [71][72][73][74], we have the following gravitational field equations: 8π where the overhead dot denotes the derivative with respect to the cosmic time t.
With the help of FLRW metric and the Bianchi identities by respecting the Cosmological Principle embodied in the FLRW metric we have the following generalized local conservation law: If we put p m = (γ − 1)ρ m and p = ω de ρ , where we assume that ω de → −1 in the above equation we have the following balanced conservation equation: Since ρ is taken to be constant the right hand side of the above equation vanishes to give the following equation: Motivated by the work of Aleksander Stachowski et al. [62], let us consider that the cosmological constant evolves with time through the hubble parameter H with (H ) given in the form of Taylor series with respect to H .
where 0 = (H )| 0 and α 2n s, n = 1, 2, . . . are the coefficients in the Taylor series expansion of (H ) given by (19) denotes the dominant term when H is near the current value H o . The reason why only even powers of H appear in the running cosmological constant is due to the general covariance of the effective action of quantum field theories and there are no covariant structures which are odd in metric derivatives [75]. This was also explained in [18] and in more detail in [76]. The contribution of only the even powers of Hubbble parameter to the time varying (t) has also been presented in detail in [21]. In the year 2020, the particle production and the corresponding entropy increase in running vacuum model was studied and explanations about the cosmic history from the very early Universe upto the current Universe and further into the final de Sitter era was also given in [77]. The implication from these even powers of H for the early Universe and for the current Universe is clearly explained in [77]. At early stage of the Universe the presence of H 4 in the vacuum energy density triggers inflation while H 4 term becomes irrelevant in the late epochs of the current Universe [77]. The inflationary models triggered by the power H 4 were first proposed phenomenologically in [78] and later H n situation is also studied in [77,79]. In the subsequent studies, they were first supported from string theory [75,80]. Recent analysis has been made on the gravitational and chiral anomalies in the running vacuum Universe and the study of the matter-antimatter asymmetry is also given in [80]. Also, the study of the stringy-running-vacuum-model inflation from primordial gravitational waves and stiff axion matter to dynamical dark energy is done in [75]. In addition, a detailed analysis of the vacuum energy density of quantum field theory in F L RW spacetime is also presented in [28]. Also they have proven from first principles that a series of even powers of H emerges, hence predicting a generic mechanism of inflation fully based on quantum field theory in curved space time [28]. Thus it is no longer an ad hoc assumption but, it is a result of a calculation in quantum field theory in curved backgrounds. Now, using (19) in (15), we obtain the following relations: To set up the dynamical system we consider a new variable x = ( H 8π G ) 2 and substitute y = ρ m for our convenience in the subsequent analysis. With this substitution we can express (20) in terms of x and y as follows: where Since we want to do qualitative analysis on the current Universe and its behaviour in the late time, we consider only a few powers of H beyond the term C o so as to ensure a better CDM limit. All the other terms involving higher powers of H beyond H 2 term are neglected as their contribution is completely negligible at present [37,77,78].
Using (21) and the newly introduced variables in the above field equations, we obtain the following set of ordinary differential equations which will represent the required dynamical system: where = ln a denotes the logarithmic time with respect to the scale factor a. The overhead dash denotes the derivative with respect to while the overhead dot denotes the derivative with respect to cosmic time t.
that is, To analyse stability, firstly we need to find the fixed points of the system. For this we equate x = 0, y = 0, that is, This implies This implies either y = 0 or γ = 0. We can also have y → 0 in evaluating the fixed point. We need to observe both the possibilities and their implications to the evolving cosmological scenario. When y = 0 in the expression of x above we get . Again when γ = 0 then from (18) we see that ρ m = constant. Let us suppose that ρ m = ξ , that is, y = ξ . Then the second fixed point we have obtained for the case of γ = 0 is . When we consider y → 0 we will obtain a special case of non-hyperbolic fixed points called a normally hyperbolic fixed point which is actually a set of non-isolated fixed points. For normally hyperbolic fixed points stability is decided by the sign of real part of the remaining eigenvalue even if one of the eigenvalues of the Jacobian matrix vanishes. So when we choose y → 0 then we can write the fixed point as . Now let us evaluate the Jacobian matrices J F 1 , J F 2 and J F 3 at the respective fixed points to study the stability of the system. Let The Jacobian matrix at the respective fixed points are given by The above matrices are upper triangular matrices. We all know that the eigenvalues of the Jacobian matrices are given by the diagonal entries. So, the eigenvalues of The fixed points F 1 and F 3 are hyperbolic for γ = 0 as none of the eigenvalues vanishes. When γ = 0, As all the eigenvalues of J F 1 and J F 3 possess negative values for γ = 0, α 2 < 3, F 1 and F 3 are stable fixed points. If y → 0 is considered, though E V J 3 2 = 0 F 3 is still stable as the remaining eigenvalue (α 2 − 3) is negative for α 2 < 3. The fixed points F 1 and F 3 are stable and behaves as an attractor for α 2 < 3. Figures 1 and 2 shows the phase plot of F 1 for γ = 4 3 and γ = 2 respectively with α 2 = 2 < 3 where all the nearby trajectories are attracted towards it. When α 2 > 3, the eigenvalues of J F 1 possess opposite signs which shows that F 1 behaves as a saddle fixed point. Figure 3 shows the phase plot of the system for α 2 = 4 > 3 where trajectories in some directions are attracted towards F 1 while trajectories along some other directions are repelled away from it. For the fixed point F 2 we see that J F 2 is non-hyperbolic as one of the eigenvalues, namely, E V J 2 2 = 0. For non-hyperbolic fixed point F 2 we can not analyse stability using the above linear stability theory. Since it is a two dimensional dynamical system we can use the notion of perturbation function and spectral radius of the Jacobian matrix for the non-hyperbolic fixed point F 2 to analyse the stability. In the subsequent paragraph we will show the stability analysis using these methods.

A. Stability analysis for F 2 using the concept of Spectral radius:
Let's rewrite the Jacobian matrix at the fixed point F 2 as follows: The spectral radius of a matrix is the maximum of the absolute values of all the eigenvalues of the matrix. The stability of a fixed point (x, y) of a dynamical system can be determined by the value of spectral radius of its Jacobian matrix evaluated at the fixed point. The notion of spectral radius in discussing stability of a fixed point has been given in detail in [68].
The spectral radius of the above Jacobian matrix is given by By theorem [ [68], page 221], F 2 will be locally asymptotically stable if σ J F 2 < 1. And we can not determine stability when σ J F 2 = 1 and hence Center manifold theory is very useful in analysing stability. With reference to [ [68], page 200], spectral radius will be less than unity if and only if From the above arguments, F 2 is locally asymptotically stable for 3 < α 2 < 4 or 2 < α 2 < 3. It can be noted that we have assume α 2 = 3 here so that we can study our system with fixed points in finite phase plane.

B. Stability analysis for F 2 using the concept of Perturbation function:
To analyse stability in a simpler way we find perturbation function along each axis as a function of logarithmic time . It is noted that while studying perturbation along x-axis we assume y = 0 as we are analysing only along x-axis. We can make the interval where α 2 lies finer by analysing the stability from this side of perturbation function. Now to find the perturbation function we perturb the system by a small amount, that is, where η x and η y represent small perturbations along x and y axes respectively. With these perturbed system, (22) and (23) takes the following form: Solving the above differential equation we obtain η x as a function of logarithmic time, as follows: Similarly, When α 2 < 3, as tends to infinity the perturbation along x-axis, η x evolves to a constant value which is ξ (α 2 −3)8π G . In the above expression of η y if we consider → ∞, we get ∞ ∞ form. So we can apply L Hospital's rule of finding limit in the expression of η y to obtain its limiting value as −ξ for any value of γ . We can also directly put γ = 0 in (23) to get η y = 0 and obtain η y =constant. But by doing so we won't be able to show the nature of η y in terms of and further with (26) we can achieve the constant value towards which η y evolves in a finer way. As perturbation along both the axes evolve to a constant value when α 2 < 3, we conclude that F 2 is stable for α 2 < 3 and it is locally asymptotically stable for 2 < α < 3. If = {α 2 : η x → 0 or a constant ∧ η y → 0 or a constant}, then F 2 is stable for any α 2 ∈ where = (−∞, 3). The perturbation plots shown in Fig. 4 shows the variation of perturbation function along y axis with respect to for F 2 . From Fig. 4 we see that when γ = 0, η y becomes a constant function, but if γ = 0 then as → ∞, η y takes ∞ ∞ form. So by applying L Hospital's rule as → ∞, η y tends to −ξ which is a constant value. Figure 5 shows that the perturbation along x-axis tends to a constant value, namely, In the plot shown in Fig. 5 we take ξ = 1, 8π G = 1 and α 2 = 2.5 < 3 to show that η y tends to In terms of the variables x and y we obtain the value of effective equation of state ω e f f and total energy density tt as follows : where ρ tt = 24π Gx and 1+z r we express the EoS for dark energy sector ω de , the Hubble parameter H and the deceleration parameter q in terms of redshift z r as follows: , where ω de is the EoS for dark energy sector. This implies, Stable for α 2 < 3 and γ = 0, late time attractor; saddle point for α 2 > 3, unstable Stable for α 2 < 3, locally asymptotically stable for 2 < α 2 < 3 where E 1m , E 1 , b 1 , B 1 are arbitrary constants. Also we have the transition redshift, z rt as follows: At F 1 the value of EoS parameter ω de is calculated as −1 which assures the presence of negative pressure in the existing cosmological scenario with the numerical value of tt as tt ≈ 1. Thus the presence of this late time attractor contributes to our model with an accelerated expansion phase of the Universe with ω de = −1 and tt = 0.99 ≈ 1 which is in agreement with the observational data in [81]. Also when we evaluate the above cosmological parameters at the fixed point F 2 , we obtained 3 . The relative energy density at F 2 is found to be tt = 1. These results have been tabulated in Table 1. To analyse the phantom-like or quintessence like behavior, lets us discuss ω de by expressing it in terms of redshift z r .
The expression given by (30) is very useful to study the behaviour of the EoS at late time when the cosmic time t tends to infinity. We can analyse the phantom-like or quintessencelike behaviour when the redshift parameter z r approaches the fixed point value near −1. Figure 6 shows the plot for ω de with respect to z r . In this plot we see that the Universe is in the phantom-like phase at the present epoch (assuming the present value of scale factor to be unity) with the present value of EoS obtained as ω de −1.029 < −1. As z → −1 at late time it is observed that the curve approaches the phantom region. As we know that the deceleration parameter determines the accelerating or decelerating behaviour while the Hubble parameter decides the rate of expansion of the Uni- verse, it is also important to analyse the behavior of these parameters in redshift. From the plots shown in Figs. 7 and 8 respectively, we see that the Hubble parameter increases with the increase in time while the deceleration parameter decreases monotonically with deceleration-acceleration transition taking place at z rt 0.71. This value of z rt represents the the value of the transition point from early decelerated regime q(z r ) > 0(corresponding to (z r > z rt ), into current accelerated one q(z r < 0) (corresponding to (z r < z rt )) and it remains close to the C DM [61]. Also the present values of the above parameters are tabulated in Table 4 and these results are in agreement with those reported in [61,[82][83][84][85].

C. Stability at infinity and Poincaré sphere:
The detail explanation of Poincaré sphere and behavior at infinity is given in [70]. By using stereographic projection we can study the behavior of trajectories far from origin by considering the so-called Poincaré sphere where we project from the center of the unit sphere S 2 = {(X, Y, Z ) ∈ R 3 : X 2 + Y 2 + Z 2 = 1} onto the (x, y)-plane tangent to S 2 at the north pole [70] by using the transformation of coordinates given by The equations defining (X, Y, Z ) in terms of (x, y, z) are given by The critical points at infinity are mapped on the equator of the Poincaré sphere. We consider the following flow in R 2 : Let f (x, y) = 1 8π G (C o + β 1 x − (γ − 1)y), g(x, y) = −3γ y. The degree of this polynomial system is one and let f 1 and g 1 denotes the homogeneous polynomials in f and g of first degree, that is, f 1 = 1 8π G (β 1 x − (γ − 1)y), g 1 = −3γ y. In terms of the polar coordinates r , θ with x = r cos θ , y = r sin θ , we can express the above equations as Order of r in (37) as r → ∞ isī = 1 and that of (38) is j = 0. Let us denotek =ī −j = 1 − 0 = 1. And using (37), we have Then using Theorem 2.1 we find G 2 (θ ) which is also equal to the highest power term in r of the expression [86].
Solving θ for which G 2 (θ ) = 0 we get θ = nπ , where n = 0, ±1, ±2, . . . So we can conclude that G 2 (θ ) is not identically equal to zero but it becomes zero in those directions where θ takes the value nπ . Since G 2 (θ ) has at most 2 pairs of roots θ and θ + π , the equator of the Poincaré sphere has finite number of fixed points located at θ such that G 2 (θ ) = 0, that is, at θ = 0, π, π, 2π or equivalently θ = 0, π. At γ = 0, 4 3 and 2, G 2 (θ ) takes the following form: The flow on the equator of the Poincaré sphere is counterclockwise at points corresponding to polar angles {θ : The flow is clockwise at points corresponding to polar angles For γ = 4 3 the flow on the equator of the Poincaré sphere is counterclockwise at points corresponding to polar angles {θ : θ > tan −1 [(1 + α 2 )24π G])} where G 2 (θ ) > 0 and the flow is clockwise at points corresponding to polar angles {θ : θ < tan −1 [(1 + α 2 )24π G]} where G 2 (θ ) < 0. For γ = 2 the flow is counterclockwise at points corresponding to polar angles {θ : θ > tan −1 ((3+α 2 )8π G)} where G 2 (θ ) > 0 and the flow is clockwise at those points corresponding to polar angles {θ : θ < tan −1 ((3+α 2 )8π G)} where G 2 (θ ) < 0. By Theorem 2.1, the critical points at infinity for the system occur at the points (X, Y, 0) on the equator of the Poincaré sphere where X 2 + Y 2 = 1 and 8π G and g 1 (x, y) = −3γ y. Using (34), the above equation becomes Solving for X and Y from the above equations, we find that fixed point occurs at (±1, 0, 0). Also we see from the expression in (41) that for γ = 0 the flow on the equator of S 2 is clockwise for XY > 0 and counterclockwise for XY < 0. For γ = 4 3 , the flow on the equator of S 2 is clockwise for XY > 0 and −(1 + α 2 )XY > Y 2 24π G ; and the flow is counterclockwise for XY < 0. For γ = 2, the flow on the equator of S 2 is clockwise for XY > 0 and −(3 + α 2 )XY > Y 2 8π G ; and the flow is counterclockwise for XY < 0. Using Theorem 2.2, the behavior in the neighbourhood of the critical point (1, 0, 0) is topologically equivalent to the behavior of the following system, where the overhead dash denotes the derivative with respect to logarithmic time .
Putting the expressions of f, g in (42) and (43) we get The Jacobian matrix of the above system is This is a diagonal matrix. So the eigenvalues are given by the main diagonal entries, that is, m 1 = (α 2 − 3) + 3γ and m 2 = α 2 − 3. Both m 1 and m 2 are negative for α 2 < 3 and γ < |( α 2 −3 3 )|. So the critical point (1, 0, 0) behaves as a stable attractor which represents the late time accelerated expansion phase of the Universe for α 2 < 3, γ < |( α 2 −3 3 )|. As we can see when α 2 < 3 and γ = 4 3 , the critical point (1, 0, 0) behaves as a saddle point with m 1 > 0 and m 2 < 0 which is unstable representing the matter dominated phase of the evolving Universe. For α 2 > 3, both m 1 and m 2 are positive for any γ and the critical point (1, 0, 0) behaves as an unstable repeller representing the inflationary epoch of the evolving Universe. Figures 9 and 10 shows the phase plot of stable attractor as well as the unstable repeller respectively Fig. 9 The phase plot of stable attractor (0, 0) for analysing stability at infinity for case I when γ = 0, α 2 < 3 taking C o = 8π G = 1 for γ = 0. Figures 11 and 12 show the phase plot for unstable saddle point and repeller for respectively for γ = 4 3 . Since the degree of f (x, y) and g(x, y) is odd, the behavior at the antipodal point (−1, 0, 0) is exactly the same as the behavior at (1, 0, 0).
What is interesting in this model is that in infinite spacetime when we study fixed points at infinity the presence of the critical point (1, 0, 0) indicates the possibility of inflationary phase and it is represented also in the phase plots shown in Figs. 10 and 11 where all the nearby trajectories are repelled away from it.
Case II: Dynamical system analysis forĠ = 0 and ρ = constant Let us consider the following Einstein field equations: where G μν = R μν − 1 2 g μν R denotes the Einstein tensor. With general Bianchi identity ∇ μ G μν = 0, the above field equation gives the following relation: This implies that the local conservation law takes the following form which we named it mixed local conservation law: If we assume thatĠ = 0 and ρ =constant, then the above relation leads to the following equation which indicates a non-conservation of matter as G does not remain constant here: Fig. 10 The phase plot of unstable repeller (0, 0) for analysing stability at infinity for case I when γ = 0, α 2 > 3 taking C o = 8π G = 1 Fig. 11 The phase plot of unstable saddle point (0, 0) for analysing stability at infinity for case I when γ = 4 3 , Since we are inclined to qualitative study of the dynamics of the Universe, we set up a dynamical system for case-II by introducing a new variable: x = 8π G 3H 2 and putting y = ρ m . With these new variables the field equations can be rewritten as Fig. 12 The phase plot of unstable repeller (0, 0) for analysing stability at infinity for case I when γ = 4 3 Again using the Taylor series form of (H ) in the field equation 8π Gρ m + = 3H 2 , we get Now the dynamical system is represented by the following system of ordinary differential equations: Using the expression ofĠ,Ḣ and o 3H 2 we have found above, we get and In order to find the fixed points we equate x = 0 and y = 0. If x = 0, then either y = 0 or γ = 0 as x = 0 otherwise if x = 0, then (48) will be violated. Again if γ = 0 is considered then we get y = b where b is a real constant and x = a where a, b ∈ R satisfies a(b + ρ ) = 1. So the first fixed point we have obtained here is P = (a, b) where a(b + ρ ) = 1; a, b ∈ R. Now consider y = 0 when γ = 0 then x = 1 ρ , that is, Q = ( 1 ρ , 0) is the second fixed point. In studying the stability of the fixed points, Jacobian matrix of the system plays a leading role. The Jacobian matrix J 2 of the system is as follows: At the fixed points P, Q, J 2 takes the following form respectively: Since P is obtained when γ = 0, J P becomes a null matrix and hence the eigenvalues of J P are m 1 = 0, m 2 = 0. The eigenvalues of J Q are m 3 = 0, m 4 = −3γ . We see that at least one of the eigenvalues vanish at both the fixed points and hence both P and Q are non-hyperbolic. So we need to use the concept of perturbation function as it is easy to analyse the behaviour of the system from the nature of perturbation function expressed in terms of . As tends to ∞, if the perturbation along each of the axes grows then the fixed point is unstable whereas if the perturbation along each of the axes decays to zero or evolves to a constant value, then the fixed point is stable. We shall not employ Center manifold theory for two dimensional problems as it is simpler to use the method of perturbation function, but for higher dimensional problems as using Center manifold theory is one of the prominent tools to study stability of a system, we have also shown in the later part, namely, Case III of this section how the dynamics of the center manifold determines the dynamics of the entire system.

A. Stability analysis using the concept of Spectral radius of the Jacobian matrix at the respective fixed points:
The spectral radius of J P and J Q are given by Since σ P < 1, all the eigenvalues of J P lie inside a unit disc. So P is stable. When γ > 0, σ Q < 1 if γ < 1 3 and σ Q = 1 if γ = 1 3 . So, Q is stable for 0 ≤ γ < 1 3 and we can't say whether Q is stable or not if γ = 1 3 . In addition when γ = 1 3 one eigenvalue of J Q ,namely, −3γ , has absolute value equal to one, the other eigenvalue, that is, zero has absolute value less than one. In this case a bifurcation may occur where a small change in the parameter values of the system leads to a sudden qualitative change in terms of topological behavior of the system. We need to further our study from the concept of perturbations along each of the axes and study the behaviour of perturbations when → ∞.

B. Stability analysis using the concept of Perturbation function:
Let x = x P + η x and y = y P + η y , where x P , y P are the values of x, y at P and η x , η y are small perturbations along x-axis and y-axis respectively. Putting the perturbed value of x and y in the dynamical system equations (51) and (52) leads to the following relations: where c 1 is an arbitrary constant. Similarly, at fixed point Q we get where c 2 is an arbitrary constant.
As increases and tends to ∞, η y for P evolves to a constant value for all γ ∈ [0, 2] and η y for Q also converges to zero for all γ ∈ [0, 2]. Since the perturbation along each axis does not grow with the increase in , P is stable for all γ ∈ [0, 2], in particular for γ = 0. When γ = 0 η y → −b as → ∞ but if we directly put γ = 0 in the expression of η y above, η y becomes a constant function, η y = c 1 − b. Figure 13 shows the variation of perturbation along y-axis , η y with respect to as γ → 0 + for the fixed point P. From Fig. 13 we see that as γ → 0 from the right the curves gradually tends to η y = c 1 − b. Figure 14 shows that η y decreases exponentially as increases and ultimately decays to zero as tends to ∞ for Q for any positive value of γ . So it is obvious that η y → 0 as → ∞ for γ < 1 3 as determined from the concept of spectral radius. So Q is also no doubt stable for all 0 < γ < 1 3 . The relative energy density

We have calculated the value of effective equation of state
parameter ω e f f as follows : As a function of redshift parameter z r , we can express ω de , Hubble parameter H (z r ) and deceleration parameter q(z r ) respectively as follows:  This implies, where E i s are arbitrary constants. The expression for transition redshift is obtained as At the fixed point P, we get ω de = −1 − 2γ b ρ . When γ = 0 ω de at P becomes −1 and at Q ω de = −1 with tt = 1 for both P and Q which is in agreement with the observational data in [81]. Since ω de is found to be negative unity, the presence of the stable fixed point P indicates the presence of negative pressure in the developed cosmological model which contributes to our model with an accelerated expansion phase of the Universe. We tabulate the results in Table 2.
The plot for ω de versus z r is shown in Fig. 15. It is observed that the present epoch is dominated by phantom dark energy as ω de < −1 at present. Then as z r approaches the fixed point value −1, ω de tends to a fixed point value ω de −1 for any value of γ . As ω de −1, this indicates the complete dominance of the Universe in late time with purely cosmological-constant type dark energy which drives the Universe to expand with acceleration and continue expanding as cosmic time tends to infinity. This phenomena is again supported from Figs. 16 and 17 where we observe that the Hubble parameter increases with the increase in time while the deceleration parameter decreases with time and ultimately takes negative value. Since q(z r ) = 0 indicates the position of the transition point from early decelerated expansion to current accelerated expansion in the Universe, using (58) we obtain the value of transition redshift z rt 0.72 and this is supported from Fig. 17 also. The present values of these parameters are tabulated in Table 4 and these results are in agreement with the cosmological observations.
Case III: Dynamical system analysis forĠ = 0 anḋ ρ = 0 In this case both G and ρ are no longer constants, that is,Ġ = 0 andρ = 0. There are many possibilities  here. We consider the simplest one by assuming the standard local covariant conservation of matter-radiation, that is, ρ m + 3H (ρ m + p m ) = 0 [37], the relation in (46) leads to the following: We introduce the following new variables x, y, z to set up the corresponding dynamical system such that x = 8π G 3H 2 , y = ρ m , z = ρ . We take derivative of the newly introduced variables with respect to logarithmic time, and obtain the following relations: x =ẋ dt d , Using (59) in the above equation and the necessary substitutions we get z =ρ dt d , Putting the above expression of z in (61), we get the expression of y as follows: Finally putting the value of y above in (60), we get the expression of x as follows: We equate x = 0, y = 0, z = 0 using (64), (63) and (62) to obtain the fixed points. As y → 0, z → 0, then since x, y, z holds the relation 1 y+z = x, x must tend to infinity. If we view from the sequential approach of real analysis, any real sequence of the form 1 n converges to zero as n → ∞ but never equals to zero. For every > 0 there exist a positive integer m such that | 1 n −0| < for all n ≥ m, that is, in every neighbourhood of zero there contains infinite members of the sequence 1 n . Similarly when n → 0, 1 n → ∞. So as y → 0, z → 0 x must tends to infinity. To ensure that the fixed points obtained are physically feasible with the developed system, α 2 must be equal to 3 and with this consideration we can analyse our fixed points in the finite phase plane. Let us consider x = 0, y = 0, z = 0 at α 2 = 3, then as y → 0.0009, z → 0, x must also tends to a number, l = 1 (0.0009+0) = 1111. Let this fixed point be denoted by S = (x → l, y → 0.0009, z → 0). Stability of the above fixed points is determined by the eigenvalues of the Jacobian matrix J 3 of the above dynamical system which is obtained as follows: The above matrix is a 3 × 3 matrix. The eigenvalues of J 3 at the fixed point determines the stability of the fixed point. At S when γ = 0, J 3 takes the following form: The above matrix is a 3 × 3 matrix with eigenvalues 0, −16.74, −(α 2 − 3) = 0. Since some of the eigenvalues becomes zero, S is a non-hyperbolic fixed point. We analyse stability through perturbation function and Center manifold theory as it is a three dimensional problem with the fixed point as non-hyperbolic one and using these methods are more suitable.
A. Stability analysis for S using the concept of Perturbation function : We perturb the system by a small amount putting x = x F + η x , y = y F + η y , z = z F + η z where x F , y F , z F represent the values of x, y, z at the fixed point to be analysed for stability and η x , η y , η z denote the perturbations along x, y, z axes respectively. With these perturbed values in the dynamical system equations (64), (63) and (62) and necessary substitutions, we obtain the following perturbations as a function of logarithmic time : where C i , i ∈ κ are arbitrary constants and κ is the index set. Let = {α 2 : η x → 0 or c, η y → 0 or c, η z → 0 or c as → ∞, where c ∈ R is any real constant }. If we consider only the expression of η x obtained as a function of regardless of restricting the value of α 2 , then we can see that when → ∞, η x → C 1 − l for α 2 = 3, η x → −l for α 2 > 3, η y → C 2 for any positive value of α 2 . Similarly it is seen that η z exponentially increases for α 2 > 1.67. So we fail to obtain such value of α 2 where all of these η x , η y , η z decay or evolve to a constant value as tends to infinity. So is an empty set. Only when all of these η x , η y and η z decay to zero or tends to a constant value when → ∞, we can conclude that the fixed point is stable otherwise unstable if at least one of them go on increasing as → ∞. For S to be stable should not be an empty set. Figures 18, 19 and 20 show the perturbation plots for S at γ = 0. From Fig. 18, as α 2 → 3 − , the slope of the curve gradually decreases and as α 2 becomes exactly equal to 3, the slope of the curve equals zero and then as α 2 becomes just greater than 3, η x becomes an exponentially decreasing function of . So when α 2 > 3 as → ∞, η x exponentially decreases and evolves to a constant value, namely, −l. Figure 19 shows that η y → 0 as → ∞ for γ = 0 and any value of α 2 . But from Fig. 20 it is clear that when α 2 ≥ 3, η z exponentially increases as increases and continue to grow as → ∞. So S is unstable for any value of α 2 . Hence, S is unstable for α 2 = 3 also. In this case III, we have already presumed α 2 to be equal to 3 in order to ensure that the fixed point S obtained above is in the finite phase plane and remain physically feasible with respect to the dynamical system we have set up. So using  Fig. 19 Variation of η y with respect to for S at γ = 0 the above arguments we conclude that S is unstable from the side of perturbation function. We will also show the use of Center manifold theory in determining the stability of the fixed point S. Using Center manifold theory is one of the most rigorous mathematical technique to determine stability for non-hyperbolic fixed points as the nature of orbits on a center manifold reflects the nature of the system in the neighbourhood of the fixed point. To use Center manifold theory we need to transform the dynamical system equations into the standard form to study Center manifold theory. We know that S(x → l, y → 0.0009, z → 0) is a non-hyperbolic fixed point. Now using a suitable coordinate transformation we need to transformed the system in the required standard form to apply Center manifold theory for it will not change the nature of the fixed point. We present how to analyse stability using the Center manifold theory in the following section.

B. Stability analysis for S using Center manifold theory:
Firstly, we need to transform the dynamical system equations into the form required to use Center manifold theory. For this we need to shift the fixed point to origin (0, 0, 0) by doing suitable coordinate transformation as follows: In terms of this new coordinates our dynamical system equations (64), (63)and (62) with α 2 = 3 can be written as follows: Eur. Phys. J. C (2022) 82 :863 Fig. 23 The variation of η z with respect to at S at γ = 2 The Jacobian matrix of the above system at origin is The above Jacobian matrix has zero determinant which means at least one of the eigenvalues has become zero. To find the eigenvalues say m i we solve the characteristic equation det (J X =0,Y =0,Z =0 −m I ) = 0 and obtain m 1 = 0, m 2 = −0.05l, m 3 = −16.81. The minimal polynomial that annihilates J X =0,Y =0,Z =0 is given by m(m + 0.05l)(m + 16.81). As the linear factors occur exactly once in the minimal polynomial, J (X =0,Y =0,Z =0) is diagonalisable. To diagonalise J (X =0,Y =0,Z =0) to obtain the required form to use Center manifold theory, we need to find the stable subspace E s generated by the eigenbasis associated with the negative eigenvalues, the center subspace E c generated by the eigenbasis associated with the zero eigenvalue of above Jacobian matrix. The eigenspace associated with zero eigenvalue can be found out by solving for x 1 , x 2 , x 3 in the following matrix equation: where I 3×3 and O 3×3 represents the identity matrix and null matrix respectively. Solving the above equations we get the eigenbasis as Similarly we find the eigenbasis associated with the eigenvalues −0.05l and − 16.81 so that we can write stable subspace (E s ) as follows: Both E c and E s are the subspaces of R × R × R. Let us define a matrix P whose column vectors are formed by the above eigenbases as follows: P is a non-singular matrix with det (P) = −646.8l. So P is invertible matrix with P −1 as P −1 = 1 det (P) Ad j (P), where Ad j (P) denotes the adjoint of P. Therefore We again define a new co-ordinate transformation as: that is, In terms of the new coordinates U , V , W ; X , Y and Z can be expressed as follows: The definition of center manifold allows us to take h 1 and h 2 in Taylor's series form as V = h 1 (U ) = a 1 U 2 + a 2 U 3 and W = h 2 (U ) = b 1 U 2 + b 2 U 3 so that h 1 (0) = h 1 (0) = 0 and Dh 1 (0) = Dh 2 (0) = 0, where D = d dU . We then obtain the required standard form to apply Center manifold theory as follows: Now computing the above equations we obtain the following relations: The dynamics of the center manifold is given by: The tangency condition requires that By equating the coefficients of U 2 and U 3 in the tangency conditions (71) and (72), we can find the constants a 1 , a 2 and b 1 , b 2 where we unconsider all the powers of U higher than U 3 . Equating the coefficients of U 2 and U 3 in the tangency condition of V , we get a 1 = a 2 = 0 and from the tangency conditions of W comparing the coefficient of U 2 , we get Since l is a very large number, b 1 ∼ 0.8l and comparing the coefficient of U 3 we get Putting the values of a 1 , a 2 , b 1 , b 2 in the dynamics of center manifold we get where j 1 = (−54177l + 12186l 2 ) and j 2 = (14522l 2 − 32742l + 9748)(2l 2 − 2l + 0.6).
Since the first term of U is in even power of U , we deduce instability. If suppose j 1 = 0 then we will consider the next term which is in the odd power of U . Here if j 2 is negative then, it is stable otherwise if it is positive then we again achieve instability. But in our case j 1 never equals zero. So from the side of Center manifold theory we conclude that the fixed point S is unstable. Now when we take γ = 4 3 then, (64), (63) and (62)  Now when γ = 4 3 we have the Jacobian matrix at S as follows: 3 at α 2 = 3, S becomes non hyperbolic fixed point. We need to analyse stability through perturbation function and Center manifold theory. However stability analysis using Center manifold theory is similar to the above shown. So we will only analyse through perturbation function. From (65), (66) and (67), we see that for α 2 = 3 η x tends to a constant, namely, (C 1 − l) as → ∞ but η y exponentially increases as → ∞. η z is also an exponentially increasing function of and hence it fails to decay or evolve to a constant value as → ∞. Figure 21 shows the exponential increasing nature of η y and η z at γ = 4 3 , α 2 = 3. Figure 22 shows the perturbation plot for η x as tends to infinity. So S is unstable at α 2 = 3 and γ = 4 3 . As the perturbation along each of the axes fail to decay or evolve to a constant value we conclude that S is also unstable for γ = 4 3 . For γ = 2 also we can see from (67) that the perturbation along z axis, namely, η z is an exponentially increasing function of θ . So S is unstable for any value of α 2 for γ = 2 and this is shown in Fig. 23 also.
The Jacobian matrix of the above system at the fixed point (0, 0, 0) is a null matrix which has all the eigenvalues as zero. So it is a non-hyperbolic fixed point. We will analyse the stability by finding perturbation functions along each of the axes as a function of logarithmic time by perturbing the system (80) by a small amount. If the system comes back to the fixed point following the perturbation then the system is stable otherwise if the perturbation grows to make the system moves away from the fixed point then the system is unstable. Nandan Roy and Narayan Banerjee [87] has also used the concept of perturbation function to analyse stability for nonhyperbolic fixed points for three dimensional systems where linear stability fails. Now firstly consider the expression of (80) corresponding to +y, +z and +w respectively. Then we perturb our system (80) by taking y = η y , z = η z and w = η w . dη y d = (−9 + 56.24γ − 6.4α 2 )η 2 y − 3η 3 y , Integrating both sides of the above differential equation, we get where k = (−9 + 56.24γ − 6.4α 2 ). The domain of definition D of the above function at γ = 0 is The domain of definition D of the above function at γ = 4 3 and γ = 2 respectively are as follows: With the above domain and the choice of +y on the left side of (80), we cannot analyse our system for → ∞ as becomes bounded above and unbounded below as η y tends to 0, that is, when → −∞, η y → 0. Since we want to analyse the late time behaviour of the Universe as logarithmic time → ∞ we only consider the expressions of (80) corresponding to −y, −z and −w on the left sides of (80) as follows: With this consideration we get the expression of as a function of η y as follows: When → ∞, f (η y ) → ∞ which implies η y → 0. So as → ∞ the perturbation along y-axis decays to zero. For analysing the perturbation along z and w axes we consider the expression for +z and +w from (80) and find out the Fig. 24 The variation of with respect to η y for analysing stability at infinity for case III Fig. 25 The variation of η z with respect to for analysing stability at infinity for case III expressions of η z and η w as follows: where c 1 and c 2 are arbitrary constants of integration. As tends to infinity both η z and η w tends to zero. Figures 24, 25 and 26 show the projection of perturbations along y, z and w axes respectively for system (80). Since all of η y , η z and η w decay to zero as tends to infinity, we conclude that the fixed point (±1, 0, 0, 0) is a stable critical point. pagination The expressions of total energy density tt and effective equation of state ω e f f in terms of the variables x, y, z are as follows: Fig. 26 The variation of η w with respect to for analysing stability at infinity for case III Fig. 27 The plot for ω de with respect to z taking α 2 = 3.05 > 3 and γ = 4 3 for E 3 = 1 2.8 < 1 e , e = 2.7 for cosmological model associated with case III The EoS parameter for dark energy sector ω de , the Hubble parameter H (z r ) and the deceleration parameter, q(z r ) Fig. 29 The graphical nature of q(z r ) over z r for the cosmological model associated with case III. The free parameters E 5 = 1.45, γ = 4 3 , α 2 = 4 > 3 and c 1 = 0.0085, E 4 = 0.0005 and E 6 = 3.04 expressed as functions of redshift z r are as follows: This implies, where E i s are arbitrary constants. The transition redshift can be obtained by using the following relations: The graphical natures of cosmological parameters ω de , H (z r ) and q(z r ) are shown in Figs. 27, 28 and 29 respectively. From the plot, we find that the curve for ω de evolves from the phantom phase and remains in the phantom region at present where ω de < −1 and then it approaches the purely cosmological constant type dark energy dominated epoch where ω de = −1 and finally in the late time the Universe is dominated by quintessence where ω de > −1. At present time when z r = 0 assuming the present value of scale factor to be unity, we find the present value of the EoS parameter as ω de −1.061 < −1 which clearly depicts the phantom-like behaviour at present. When α 2 = 3, the value of ω de = −1 at present time as well as at late time, this indicates that the Universe is associated with purely cosmological constant type dark energy and hence the model behaves as a cosmological constant model. Also the curve for H (z r ) monotonically increases as time increases while the curve for q(z r ) decreases as time increases with the transition redshift z rt found to be z rt 0.70. This means as z r passes this point z rt the Universe begins to undergo accelerated expansion with deceleration parameter value being negative. The present value of deceleration parameter q(z ro ) is obtained as -0.35 and the present values of these parameters, namely, ω de (z ro ), H (z ro ), q(z ro ) are tabulated in Table 4 and they all are in agreement with the cosmological observations. Hence the present model associated with case III supports the fact that the Universe is currently in accelerated expansion epoch and will continue to expand in the late time.

Conclusion
In this work we have presented three accelerating cosmological models through dynamical systems approach in the presence of a time varying cosmological constant term under three different cases based on the possibilities of G and ρ . We use dynamical systems approach and find the fixed points to analyse stability as well as late time behavior of the evolving Universe. The cosmological parameters, namely, ω de (z r ), H (z r ), q(z r ), z rt are expressed in terms of redshift z r and their graphical behaviour with the change in redshift have been analysed. The cosmological models associated with case I, case II and case III show phantom behaviour at the present time with ω de < −1. It is found that both analytical and geometrical findings agree with the fact that the Universe is in the accelerated expansion phase with decelerationacceleration transition taking place at values of z r 0.72 which are also tabulated in Table 4 for all of the three cases and they are quite consistent with the spatially flat C DM model, the standard model of cosmology where the cosmological constant is the dark energy. Along with the dynamical system analysis in the finite phase plane what is more interesting in this work is the stability analysis of the systems for critical points at infinity by using the concept of Poincaré sphere where we use stereographic projection to study the behavior of trajectories far from origin.   q(z r 0 ) −0.520 −0.651 −0.350 −1.08 ± 0.29 [84] q(z rt ) 0.71 0.72 0.70 0.732 ± 0.0174 [61,85] In case I where both G and ρ are taken to be constants, we observe that the time varying cosmological constant model can be represented by a two dimensional dynamical system. It has three fixed points F 1 , F 2 , F 3 in the finite phase plane out of which F 1 behaves as a late time attractor as shown in Fig. 1 when α 2 < 3 with ω de = −1 and tt = 1. The presence of this stable attractor in the cosmological model associated with case I assures the presence of negative pressure representing the accelerated expansion phase of the Universe.Also more interestingly, when we analyse stability for the critical points at infinity using poincarè sphere, the critical point (1, 0, 0) behaves as an unstable repeller when α 2 > 3 for any value of γ . This repelling behaviour is also shown in the phase plots of Figs. 10 and 12. The presence of this repeller adds up to this model with the inflationary phase of the evolving Universe. The EoS parameter for dark energy sector remains in the phantom region at present with ω de (z ro ) −1.029 < −1 as seen from Fig. 6 while H (z r ) increases monotonically with the increase in time with H (z ro ) 73.5 (shown in Fig. 7). As seen from Fig. 8, the deceleration parameter q(z r ) decreases monotonically and vanishes at z rt 0.71 and remains negative after this transition point. It is at this point where the transition from early decelerated expansion to current accelerated expansion in the Universe takes place and this shows that the model associated with this case remains close to C DM. The present values of these parameters are tabulated in Table 4 and they all are in agreement with the cosmological observations. Thus, the behaviours of these cosmological parameters depict the accelerated expansion epoch of the model Universe behaving as phantom-like at present time.
In case II when G no longer remains constant, we obtain two non-hyperbolic fixed points P and Q both of which are stable attractors for 0 ≤ γ < 1 3 with the numerical values of ω de and tt obtained as − 1 and 1 respectively. From the expression of ω de as a function of redshift z r along with the graphical interpretation shown in Fig. 15, it is observed that as the redshift parameter z r approaches − 1, ω de tends to a fixed point value ∼ −1. So this particular cosmological model associated with case II shows the dominance of a purely cosmological-constant type dark energy [80] in the late time. Presently, this model shows a phantom-like phase with present value of EoS found to be ω de (z ro ) −1.06 < −1. The Hubble parameter also shows a monotonically increasing nature with the increase in cosmic time with the present value of Hubble parameter noted to be H (z ro ) 75. In this model the transition from deceleration regime to acceleration regime takes place at z rt 0.72. After this transition point the value of q(z r ) remains negative which depicts that the model Universe will continue to expand with acceleration in the late time also. However in case III, when both G and ρ are taken to be non-constants, we find that we can extend the system to a three dimensional dynamical system problem. Here we obtain a fixed point S in the finite phase plane which is physically feasible with the developed dynamical system with the value of α 2 presumed to be equal to 3. From the side of perturbation function approach as well as Center manifold theory, we get S to be unstable fixed point for any value of γ . But the main thing that distinguishes this approach is the presence of the stable fixed point (±1, 0, 0, 0) at infinity lying on the northpole of the Poincaré sphere S 3 with the use of projective geometry being carried over to higher dimension in R 3 . This fixed point behaves as a stable attractor in infinite spacetime. The model associated with this case evolves from the phantom-phase and gradually approaches the quintessence region in late time (shown in Fig. 27). At present it behaves as phantom-like with ω de −1.061 < −1. This model also supports the fact that the Universe is in acceleration era as we can see from Figs. 28 and 29 where H (z r ) tends to infinity as z r tends to −1 while q ( z r ) transits from positive to negative at z r 0.70. This point represents the transition redshift z rt where the transition from early decelerated expansion to present accelerated expansion of the Universe takes place. Throughout the entire work we find that the developed cosmological models associated with case I, case II and case III strongly support the accelerated expansion phenomena thereby depicting that the Universe is in the phase of expansion with acceleration and will continue to expand in late time also.