Dynamical systems of cosmological models for different possibilities of $G$ and $\rho_{\Lambda}$

The present paper deals with the dynamics of spatially flat Friedmann-Lemaitre-Robertson-Walker (FLRW) cosmological model with a time varying cosmological constant $\Lambda$ where $\Lambda$ evolves with the cosmic time (t) through the Hubble parameter (H). We consider that the model dynamics has a reflection symmetry $H \rightarrow -H $ with $\Lambda(H)$ expressed in the form of Taylor series with respect to H. Dynamical systems for three different cases based on the possibilities of gravitational constant G and the vacuum energy density $\rho_{\Lambda}$ have been analysed. In Case I, both G and $\rho_{\Lambda}$ 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 axis with respect to cosmic time and Poincare sphere. In Case II, we have dynamical system analysis for G=constant and $\rho_{\Lambda} \neq $ constant where we study stability by using the concept of spectral radius and perturbation function. In Case III, we take $G \neq$ constant and $\rho_{\Lambda} \neq$ 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 axis, Center Manifold Theory and stability at infinity using Poincare 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. Our model agrees with the fact that the Universe is in the epoch of accelerated expansion. The EOS parameter $\omega_{eff}$, total energy density $\Omega_{tt}$ are also evaluated at the fixed points for each of the three cases and these values are in agreement with the observational values in [1].


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 [2,3,4,5,6,7,8,9,10,11,12]. 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 [13]. In this context of accelerated expansion the theory of general relativity (GR) modified by a cosmological constant term Λ, which is known as the famous Λ CDM model is one of the most popular one [14].
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 [15,16,17]; and "the cosmic coincidence problem [18]. To mend up these issues, running Λ cosmological models have been developed.
Shapiro et al. [19,20,21,22] made the first development regarding the scal-ing evolution of the cosmological constant. 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 [22,23,24], Λ(t) cosmology induced by a slowly varying Elko field [25], a running vacuum in the context of supergravity [26], etc. 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 energy-momentum conservation law holds then there should not be any variation in the gravitational coupling with respect to the space time or otherwise the ordinary energymomentum conservation law will be violated [27,28]. In the light of Dirac's idea [29,30,31] which propose that some of the fundamental constants cannot remain constant forever, it is essential to do some modifications in GR field equations [32,33] 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 [34,35,36].
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 [36,37,38,39,40,41,42,43,44,45,46,47,48,49]. As there are no rigorous proves that indicate whether the cosmological constant is running or not [49], one can study the cosmological implications of different possible theoretical assumptions of Λ term. Motivated by the quantum field theory [20,21,50] and some theoretical motivations [22,23] about the varying Λ form. Aleksander Stachowski, Marek Szydtowski [51] have also studied the dynamics of cosmological models with various forms of Λ(t).
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 analyze 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 section 1 we have given the introduction part, in section 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 section 3, we have three cases. In case I of section 3 we show the setting up of cosmological equations and dynamical system analysis where both G and ρ Λ are taken to be constant 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 axis 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 section 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. S. Surendra et al. [52] have also used this approach to study cosmological models in the presence of a scalar field using different forms of potential. From [52] we can also notice that in three dimensional dynamical system we can analyse stability by analysing the nature of perturbation along each of the axis. A dynamical system is generally written in the form of the following [53]: where x = (x 1 , x 2 , ......x n ) is an element of state space X ⊆ R n and the function The overhead dot denotes the derivative with respect to cosmic time, t. The for all n ∈ N. Otherwise, the fixed point x o will be called unstable; (ii) the fixed point x o is said to be attracting if there exists ζ > 0 such that (iii) the fixed point x o is said to be locally asymptotically stable if it is both stable and attracting. If in the previous item ζ = ∞, then x o is said to be globally asymptotically stable. Jacobian matrix of dynamical system at a fixed point: The Jacobian matrix of the dynamical system given in (1) at a fixed point x o is given by where δfi δxi , i = 1, 2, ..., n denotes the first partial derivative of f i with respect to the i th component x i of the element x = (x 1 , x 2 , ...x n ) ∈ X ⊆ R n .
Linear stability theory is one of the simplest method 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. 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. Centre 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 [53] and Carr [56] have discussed the centre 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 [56], Elaydi's book [57] and Zhang's book [58]. 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 coor-dinate transformation we can rewrite any system of the form (1) as follows: x = Ax + f (x,y), 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: Definition 2.6. [56] Centre Manifold: it can be locally represented as for a sufficiently regular function h(x) on J s and δ however small it may be.
The proofs of the existence of the centre manifold for the system (2) is also provided in [56] and he has given the dynamics of the system (2) restricted to the centre manifold as follows: for sufficiently small v ∈ R c .
[59] 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. [59] 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. [59]
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 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. [59] 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 fourvector 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 Λ = −ρ Λ 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 [60,61,62,63], we have the following gravitational field equations: 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 Λ = −ρ Λ and p m = (γ − 1)ρ m 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: In addition let us consider that there is reflection symmetry with respect to H, that is , H → −H. So, if the system has λ(t) as its solution then, λ(−t) is also a solution of the system. As a result only the terms containing even powers of H are present in the above power series (19). Shapiro and Solà [22] have also considered in detail the contribution of only the even powers of Hubbble parameter to the time varying Λ(t).
Using (19) in (14), we have where Λ 0 = Λ(H)| 0 and α n ′ s, n = 2i, i = 1, 2, ... are the coefficients in the Taylor series expansion of Λ(H) given by α n = 1 n! d n Λ(H) dH n | 0 , n = 2i, i = 1, 2, ... To set up the dynamical system we consider the following set of new variables: x = ( H 8πG ) 2 and y = ρ m . With this substitution we can expressed (20) in terms of the new set of variables as follows: where 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.
x ′ = 1 8πG Here 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 are neglected as their contribution is completely negligible at present [64] 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 . 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 F 2 = ( −Co−ξ (α2−3)8πG , y = ξ). 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 eigenvalue 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 F1 , J F2 and J F3 at the respective fixed points to study the stability of the system.
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 The fixed points F 1 and F 3 are hyperbolic for γ = 0 as none of the eigenvalues vanishes.

Fig. 1 and
When α 2 > 3, the eigenvalues of J F1 possess opposite signs which shows that F 1 behaves as a saddle fixed point. Fig. 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 F2 is non-hyperbolic as one of the eigenvalues, namely, EV J2 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 [57].
The spectral radius of the above Jacobian matrix is given by 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, x = −co−ξ (α2−3)8πG + η x and y = ξ + η y , 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 con-sider Θ → ∞, 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 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. Fig. 5 shows that the perturbation along x−axis tends to a constant value, namely, ξ (α2−3)8πG when α 2 < 3. 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 ω ef f and total energy density Ω tt as follows: where vacuum energy density, Ω Λ = the observational data in [1]. Also when we evaluated the above cosmological parameters at the fixed point F 2 , for any value of α 2 and ξ we obtained ω ef f = −1. The relative energy density at F 2 is found to be Ω tt = 1. The above results have been tabulated in TABLE I:   shows the phase plot for stable F 1 at γ = 2, α 2 < 3.  point. Fig.4 shows variation of η y with respect to Θ for F 2 .
0.5  C. Stability at infinity and Poincaré sphere: The detail explanation of Poincaré sphere and behavior at infinity is given in [59]. 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 the (x, y)−plane tangent to S 2 at the north pole [59] 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 : 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, 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 (30) as r → ∞ isī = 1 and that of (31) isj = 0. Let us denotē Then using Theorem 2.1 we find G 2 (θ) which is also equal to the highest power term in r of the Θ ′ expression [65].
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 {θ : where f 1 (x, y) = (α 2 − 3)x − (γ−1)y 8πG and g 1 (x, y) = −3γy. Using (27), 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 (33) that for γ = 0 the flow on the equator of S 2 is clockwise for XY > 0 and counterclockwise for XY < 0. For gamma = 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 gamma = 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, Putting the expressions of f, g in (34) and (35) we get The Jacobian matrix of the above system is 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). Fig. 8 and Fig. 9 show the phase plot for unstable saddle point and repeller respectively. Fig . 6 shows 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. Fig. 7 shows 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.  at infinity for case I when γ = 4 3 , α 2 < 3 taking C o = 8πG = 1. Fig. 9 shows the phase plot of unstable repeller (0, 0) for analysing stability at infinity for case I when γ = 4 3 , α 2 > 3 taking C o = 8πG = 1.
Case II-Dynamical system analysis forĠ = 0 and ρ Λ =constant Let's rewrite the General Relativity field equations (11) as follows: 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: But if we takeĠ = 0 as well asρ Λ = 0 assuming the standard local covariant conservation of matter-radiation (18), (38) leads to the following equation: Since we are inclined to qualitative study of the dynamics of the Universe, we set up a dynamical system for case-II by introducing new variables: With these new variables the field equations can be rewritten as 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 (41) 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, 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 alone 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 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 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 (44) and (45) 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. Fig. 10 shows the variation of perturbation along y−axis , η y with respect to Θ as γ → 0 + for the fixed point P . From Fig. 10 we see that as γ → 0 from the right the curves gradually tends to η y = c 1 − b. Fig. 11 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 γ = 4 3 also which is 1 3 as determined from the concept of spectral radius. So Q is also no doubt stable for all 0 < γ < 1 3 . We have calculated the value of effective equation of state parameter ω ef f = −1 − γxy and relative energy density Ω tt = Ω m + Ω Λ , where Ω m = xy, Ω Λ = Λo 3H 2 + α2 3 = 1 − xy. At both the fixed points P and Q, we get ω ef f = −1, Ω tt = 1 which is in agreement with the observational data in [1]. Since ω ef f 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 tabulated the results in     Fig. 11 shows the variation of η y with respect to Θ for Q at Case III-Dynamical system analysis forĠ = 0 andρ Λ = 0 In this case both G and ρ Λ are no longer constants, that is,Ġ = 0 anḋ ρ Λ = 0. The relation in (38) now becomeṡ We introduce the following new variables to set up the corresponding dynamical system: 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: Using (46) in the above equation and the necessary substitutions we get Putting the above expression of z ′ in (48), we get the expression of y ′ as follows: Finally putting the value of y ′ above in (47), we get the expression of x ′ as follows: The expression of total energy density Ω tt and effective equation of state ω ef f in terms of the variables x, y, z is as follows: where p tt = (γ − 1)y − z and ρ tt = y + z. We equate x ′ = 0, y ′ = 0, z ′ = 0 using (51), (50) and (49) 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 analyzed for stability and η x , η y , η z denote the perturbations along x, y, z axes respectively. With these perturbed values in the dynamical system equations (51), (50) and (49) and necessary substitutions, we obtain the following perturbations as a function of logarithmic time Θ: for any γ and α 2 = 3. (54) C 2 e (57−6.4α2)Θ , γ = 4 3 ; C 2 e (88.5−6.4α2)Θ , γ = 2.
where C i , i ∈ κ are arbitrary constants and κ is the index set.
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. Fig. 12, Fig. 13 and Fig. 14 show the perturbation plots for S at γ = 0.
From Fig. 12, 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. Fig. 13 shows that η y → 0 as Θ → ∞ for γ = 0 and any value of α 2 . But from Fig. 14 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 physically feasible with respect to the dynamical system we have set up. So using 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. Center manifold theory is one of the most powerful tools 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 (51), (50)and (49) with α 2 = 3 can be written as follows:  The Jacobian matrix of the above system at origin is 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 ) Adj(P ), where Adj(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 We then obtain the required standard form to apply central manifold theory 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 (60) and (61), 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).
Now when γ = 4 3 we have the Jacobian matrix at S as follows: theory. However stability analysis using Center manifold theory is similar to the above shown. So we will only analyze through perturbation function. From (54), (55) and (56), 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 Θ → ∞. Fig. 15 shows the exponential increasing nature of η y and η z at γ = 4 3 , α 2 = 3. Fig. 16 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 axis 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 (56) 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. 17 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 axis as a function of logarithmic time Θ by perturbing the system (69) 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 [66] has also used the concept of perturbation function to analyse stability for non-hyperbolic fixed points for three dimensional systems where linear stability fails. Now firstly consider the expression of (69) corresponding to +y, +z and +w respectively. Then we perturbed our system (69) by taking y = η y , z = η z and w = η w .
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: (−∞, k 3 ), α 2 > 16.15, k = 103.4 − 6.4α 2 < 0 With the above domain and the choice of +y on the left side of (69), 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 (69) corresponding to −y, −z and −w on the left sides of (69) 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 (69) and find out the expression 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. Fig. 18, Fig. 19 and Fig. 20 show the projection of perturbation along y, z and w axes respectively for system (69). Since all of η y , η z and η w decays to zero as Θ tends to infinity, we conclude that the fixed point (±1, 0, 0, 0) is a stable critical point.  Fig . 18 shows the variation of Θ with respect to η y for analysing stability at infinity for case III. Fig. 19 shows the variation of η z with respect to Θ for analysing stability at infinity for case III.

Conclusion
In this work we have presented  1) and (2) have supported these analytical results. With the notion of spectral radius we obtained a finer region of α 2 where F 1 is stable, that is, 2 < α 2 < 3. Fig. 3 shows that behaves as an unstable repeller representing the inflationary epoch of the evolving Universe. Fig. 6 and Fig. 7 show the phase plot of the stable attractor and the unstable repeller respectively. For γ = 4 3 , m 1 > 0 and m 2 < 0 when α 2 < 3 and the critical point (1, 0, 0) behaves as a saddle point which is unstable representing the matter dominated phase of the evolving Universe. When α 2 > 3, both m 1 and m 2 are positive and the critical point (1, 0, 0) behaves as an unstable repeller. Fig. 8 and Fig. 9 also support the above analytical results for γ = 4 3 . For γ = 2, the behavior is same as that of γ = 4 3 . Since the degree of the polynomial system 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). In case II of section 3, we present the case when ρ Λ = constant but G no longer remains constant. By introducing new variables, we represent the model with a two dimensional dynamical system where we obtain two non-hyperbolic fixed points P, Q. We present the stability analysis of these fixed points by using spectral radius as well as perturbation function where we have found that both are stable for γ ∈ [0, 1 3 ) with Ω tt = −1 and effective equation of state ω ef f = −1. Also for P , both η x and η y converge to a constant value as Θ tends to infinity. When γ = 0 η y → −b as Θ → ∞ but if we directly put γ = 0 in the expression of η y , it becomes a constant function, that is, η y = c 1 − b. Fig. 10 shows the variation of perturbation along y−axis , η y with respect to Θ as γ → 0 + for the fixed point P . From Fig. 10 we see that as γ → 0 from the right the curves gradually tends to η y = c 1 − b. For Q, η x evolves to a constant value and η y decays to zero as Θ gradually increases and tends to infinity as shown in Fig. 11 for γ < 1 3 . So both the fixed points are stable which gives the dark energy model which forms the strong base for the fact that the Universe is undergoing not just expansion but expansion with acceleration. When we take both G and ρ Λ to be non-constants, then we see from case III of section 3 that we can extend the system to a three dimensional problem. We have analysed the system when α 2 = 3 under three different values of γ, that is, γ = 0(dark energy model),γ = 4 3 (radiation dominated model), γ = 2(stiff fluid model) and study the system about its stability and corresponding cosmological implications. At γ = 0 the fixed point S is nonhyperbolic as some of the eigenvalues of the Jacobian matrix vanishes. Since S is non-hyperbolic, we do the stability analysis by studying how the perturbation along each of the three axis vary with the increase in Θ. As the set Φ = φ, S is unstable. Fig. 12, Fig. 13, Fig. 14 shows the perturbation plots for S. We have also used Center manifold theory to analyze stability by using a suitable coordinate transformation where we obtain the standard form to apply Center manifold theory. As the dynamics of the center manifold is unstable we deduce that S is unstable. From both approaches we find that S is unstable. For γ = 4 3 as well as γ = 2, S is non-hyperbolic and unstable. Fig. 15 and Fig. 16 show the perturbation plots of S for γ = 4 3 . The perturbation function along each of the axis fail to decay or evolve to a constant value as Θ → ∞ which shows that S is unstable. Fig. 17 shows that η z continues to increase exponentially as Θ increases which indicates that S is unstable for γ = 2 also. To analyse stability at infinity we use the concept of Poincaré sphere as any polynomial system in rectangular coordinates can be extended to the Poincaré sphere [65]. Here since the system is a three dimensional system, the ideas of projective geometry has been carried over to higher dimension to analyse stability for flows in R 3 [59].
The critical points at infinity occur at the points (±1, 0, 0, 0) on the equator the Poincaré sphere S 3 . Since the perturbation along each of the axis η y , η z and η w decays to zero as cosmic time Θ tends to infinity as shown in Figs. 18, 19 and 20, we conclude that the fixed point (±1, 0, 0, 0) is a stable attractor.
Throughout the entire work the developed cosmological models strongly support the fact that the Universe is in the phase of expansion with acceleration thereby depicting that our model has a deep connection with the accelerated expansion phenomena.

Declaration
The authors declare that there is no conflict of interest regarding the publication of this paper.