Degenerate Bogdanov-Takens bifurcations in a bulk viscous cosmology

Using the dynamical system theory we show that the Friedmann-Robertson-Walker (FRW) cosmological model with bulk viscous fluid in the presence of cosmological constant is equivalent to a degenerate two dimensional Bogdanov-Takens normal form. The equation of state parameter, $\omega$, the bulk viscosity coefficient, $\xi$, and the cosmological constant, $\Lambda$, define the necessary parameters for unfolding the degenerate Bogdanov-Takens system. The fixed points of the system are discussed together with the variation of their stability properties upon changing the relevant parameters $\omega$, $\Lambda$ and $\xi$. The variation of the stability properties are visualized by the appropriate bifurcation diagrams. Phase portrait for finite domain and global phase portrait are displayed and the issue of the structural stability are discussed. Typical issues such as late acceleration or inflation that can be induced by viscosity and could have relevance to observational cosmology are also discussed.


Introduction
Dynamical system tools applied to cosmology are valuable for enabling a qualitative understanding for the behavior of cosmological models. Through a careful suitable choice of the dynamical variables one can capture the global behavior of the models through finding fixed points without the need for finding explicit solutions. The phase portraits can reveal the general properties of orbits and also provide us with a wealth of information with minimal labor, for example see [4]. Not only qualitative understanding but also a quantitative one through using powerful analytical and numerical methods applied to the models under consideration. The dynamical system tools was first applied to anisotropic cosmological models as in [1,2,3] and for a recent review one can consult [5] and references therein.
Bogdanov-Takens bifurcations have been shown to occur in Bianchi IX cosmological models in the frame work of Gauss-Bonnet gravity [6]. A more recent study [7] has also demonstrated the occurrence of such a bifurcation in Friedmann-Roberston-Walker (FRW) cosmology in the presence of cosmological constant without considering viscosity. The latter study is of limited scope due to neglecting viscosity which is a real physical dissipative effect which is essential for getting certain desirable properties of Bogdanov-Takens system such as the finiteness of the number of fixed points. Up to the best of our knowledge, the works in [6,7] are the only two instances in cosmological studies where the Bogdanov-Takens bifurcations occurred. In fact, investigating and classifying all possible solutions and their stability properties in cosmological models enhances our understanding of the models. Needless to say, the identification of what kind of bifurcation we have for our cosmological models is important not only for spotting where we are in the vast landscape of dynamical systems but also for learning how to tune our models to have certain desired properties.
In the realm of cosmology, bulk viscosity provides the only dissipative mechanism consistent with isotropy and homogeneity. For simplicity, we consider a bulk viscosity model as described in the context of the Eckart formalism [8] rather than using the full causal theory of viscosity that was developed in [9,10]. Several authors have investigated the introduction of viscosity into cosmology for several reasons and motivations. For examples; in [11] the viscosity was introduced to resolve the big-bang singularity, while in [12] for finding a unified model for the dark component of universe (dark energy and dark matter) that could fit cosmological observational data like type Ia supernovae [13,14] and power spectrum [15,16]. Others as in [17,18,19] introduced viscosity as a source for deriving inflation in the early cosmology or for deriving late acceleration as in [20]. The possibility of using some sort of viscous fluid to get a unified cosmic history starting by inflation and ending by late acceleration dominated by dark energy have been investigated in [21]. Furthermore, in [22], it was shown that a bulk viscous model with constant coefficient of viscosity can give a viable coherent description of the different phases of the universe.
The bulk viscosity besides its clear physical origin as a dissipative effect, it might also entails the cosmological dynamical system with structural stability in the sense that the qualitative behaviour of the dynamical system doesn't change under small perturbation. The structural stability is a desirable property to be processed by any realistic system and thus worthy to be studied and tested through applying the proper criteria.
The paper is structured as follows: in Section 2, Friedmann equations for bulk viscous cosmology are presented and then expressed in terms of ρ fluid density and H Hubble parameter as our suggested dynamical variables. In Section 3, The basic theories and notations of dynamical system are presented and explained. The theories and techniques developed in Section 3 are applied in Sections 4, 5 and 6. Section 4 is devoted for investigating the case of perfect fluid with linear equation of state p = ωρ where p is the pressure. Section 5 is devoted to the case of perfect fluid as in Section 4 with the inclusion of a cosmological constant Λ. In Section 6, the bulk viscous fluid is introduced in the presence of cosmological constant and investigated. Thus, this last case amounts to having three parameter namely ω, Λ and ξ (viscosity coefficients). Finally Section 7 is devoted for discussion and conclusion.

Einstein Equations for Bulk Viscous Cosmology
A homogenous and isotropic cosmological model is described by Fredimann-Roberston-Walker (FRW) metric whose line element is given as, ds 2 = g µν dx µ dx ν = −c 2 dt 2 + R 2 0 a(t) 2 dr 2 1 − kr 2 + r 2 dθ 2 + r 2 sin 2 (θ) dφ 2 , where x µ is the four dimensional coordinate, x µ ≡ x 0 = c t, x 1 = r, x 2 = θ, x 3 = φ , a(t) is the scale factor, c is the speed of light and k = {0, ±1} which is the curvature index, while R 0 is a constant carrying the dimension of length. The metric tensor g µν can be easily read from Eq.(1) to be diagonal and given by, The scale factor a(t) can be determined by applying field equations of General Relativity (GR) which , in the presence of cosmological constant Λ, assumes the following form: where R µν and R are the Ricci tensor and scalar respectively. G is the universal Newton gravitational constant while c as before denotes the speed of light. As to the energy-momentum tensor T µν describing a bulk viscous fluid, it assumes the form where the viscous fluid has density ρ, pressure p, viscosity coefficient ξ and velocity U µ . Also, notice that H is the Hubble parameter defined as H ≡ a −1 da dt .
The resulting Einstein field equations stemming from Eq.(3), in the comoving frame, are; The above equations, Eqs. (5), can be written as a first order equations for H and ρ as, It is advantageous to rewrite the cosmological equations in dimensionless form by introducing dimensionless variables as, where ρ ch is a some chosen constant characteristic density and the characteristic Hubble parameter H ch is chosen such that H 2 ch = 8 π Gρ ch . Thus, the dimensionless form of Eqs.(6) would take the form, while the first equation in Eqs. (5) would assume the form, Assuming a barotropic equation of state,p = ωρ, then cosmological equations Eqs.(8) become, Notice that ω is is an equation of state parameter with physically motivated range range given by ω ∈ [−1, 1]. As examples for some typical values, we have ω = 0 (Dust), ω = −1 (Dark energy), ω = 1/3 (radiation), and ω = 1 (stiff fluid). The equations as given in Eq.(10) constitute the dynamical system representing the cosmological model with dynamical variablesρ andH that determine the state of the dynamical system. It is clear that these two dynamical variables are unbounded

Basic Theories and Notations for Dynamical System Approach
For convenience and notational simplicity we introduce the vector state x, vector parameter α and vector function f defined as follows, The system of equations given in Eqs. (10) can be written compactly as,

Fixed Point Analysis and Classification
For any generic planer system,ẋ = f (x, α) not necessarily the one given in Eq. (12), the existence of a fixed point is determined through f (x 0 , α) = 0 and then the system can be expanded around the fixed point as,ẋ where Df (x 0 , α) = ∂fi(x0,α) ∂xj is the Jacobian matrix. The stability of the planer system can be examined through the eigenvalues of the Jacobian matrix. Here and later, the eigenvalues of the Jacobian matrix are denoted by λ 1 and λ 2 (they are conjugated to each other in case of being complex) while their corresponding eigenvectors by e 1 and e 2 . The stability of the fixed point can be decided according to the following criteria: • Stable (Sink), if λ 1 and λ 2 are real negative and repulsive center in case of being complex with negative real part .
• Unstable (Source), if λ 1 and λ 2 are real positive and attractive center in case of being complex with positive real part.
• Saddle, if λ 1 and λ 2 are real and have opposite sign.
For cases where one of the two eigenvalues or both equal to zero, degenerate fixed points (or nonhyberbolic ones), the stability can't be decided without referring to the higher order terms which means the failure of the linear stability theory. Classification of non-hyperbolic fixed points can be found in [23]. In fact, theses non-hyperbolic fixed points are known to form the germs of bifurcation.
We follow Ref. [24]) in classifying the planer dynamical systems whose fixed point lies at (x, α) = (0, 0) with double zero eigenvalues λ 1,2 (0) = 0. The Jacobian of this system can be brought into the form J = 0 1 0 0 by introducing new variables (y 1 , y 2 ) related linearly to (x 1 , x 2 ). Then the entire system can be written and organized as a power series in terms of (y 1 , y 2 ) aṡ y 1 = y 2 + a 00 (α) + a 10 (α) y 1 + a 01 (α) y 2 + 1 2 a 20 (α) y 2 1 + a 11 (α) y 1 y 2 + 1 2 a 02 (α) y 2 2 + O y 3 , where the coefficients a ij (α) and b ij (α) are smooth functions of α and satisfying The nondegeneracy conditions for the system are the following, (BT.0) the Jacobian matrix ∂fi ∂xj (0, 0) = 0, In our specific case, one can introduce the linear transformation y 1 = x 1 , y 2 = − 1 6 x 2 and then Eq.(12) can be expressed in terms of y ′ s as, One can notice the absence of O y 3 terms and the coefficients a ij (α) and b ij (α) as defined in Eq. (14) assume the following forms, In order to check the nondegeneracy conditions one needs the Jacobian matrix ∂fi ∂xj corresponding to the system in Eq. (12) which is easily found to be, All nondegeneracy conditions are fulfilled except the condition (BT.2) where b 20 (0) = 0, thus the dynamical system described in Eq.(12) is a degenerate Bogdanov-Taken system.

Behavior at Infinity and Poincaré Sphere
To get an idea and visual representation for the solution behavior of the dynamical system one can draw a phase portrait composed of all possible solution curves in the (x 1 , x 2 ) plane. As a matter of fact, this visual representation is limited to a finite domain in the (x 1 , x 2 ) phase space. Thus, one should seek for an alternative visual representation in order to get an idea about the global solution behavior and specially at infinity. This global picture can be achieved by introducing the so-called Poincaré sphere [25,26] where one projects from the center of the unit sphere S 2 = (X, Y, Z) ∈ R 3 | X 2 + Y 2 + Z 2 = 1 onto the (x 1 , x 2 )-plane tangent to S 2 at either north or south pole as shown in Fig.(1). Projecting the upper hemisphere of S 2 onto the (x 1 , x 2 )-plane, then one can derive the following relations between (x 1 , x 2 ) and (X, Y, Z), These clearly define a one-to-one correspondence between points (X, Y, Z) on the upper hemisphere of The circle x 2 1 + x 2 2 = a 2 on the (x 1 , x 2 )-plane corresponds to points on the circle The circle at infinity of (x 1 , x 2 )-plane corresponds to the equator of S 2 . The whole orbits induced by the dynamics described by Eqs. (12) can be mapped onto the upper hemisphere of the Poincaré sphere which is difficult to draw. In contrast, the orthogonal projection of the upper hemisphere of the Poincaré sphere on the unit disk in the (X, Y ) plane is much easier to draw and still captures all of the information about the behavior at infinity. Such a kind of flow on the unit disk , X 2 + Y 2 < 1, when drawn is called a global (or compact) phase portrait. It is possible to obtain the dynamical system in terms of (X, Y ) that corresponds to the dynamical system given in Eqs. (12) and after simple algebra one can get, The determination of the fixed points, at infinity, is rather involved if one works in terms of the coordinates (X, Y, Z). Fortunately, there is a simpler approach where one can introduce plane polar coordinates (r, θ) where x 1 = r cos θ and x 2 = r sin θ and the dynamical system represented by Eqs.(12) takes the following form,ṙ = cos θ f 1 (r cos θ, r sin θ, α) + sin θ f 2 (r cos θ, r sin θ, α) , [cos θ f 2 (r cos θ, r sin θ, α) − sin θ f 1 (r cos θ, r sin θ, α)] .
Assuming f 1 and f 2 are multinomial in x 1 and x 2 and organized as, where the integer superscripts, in f ′ s, indicate the power of the associated multinomial and m is the maximum power in the expansion. Then as r → ∞ the evolution of θ is dominated by terms of maximum power f m 1,2 (x, y, α) * , contained in the expansion of f 1,2 (x, y, α), leading tȯ θ ≈ 1 r cos θ f m 2 (r cos θ, r sin θ, α) − sin θ f m 1 (r cos θ, r sin θ, α) , Furthermore, one can factor r from Eq.(23) since it doesn't affect the sign ofθ to get, θ ∼ G m+1 (θ) = cos θ f m 2 (cos θ, sin θ, α) − sin θ f m 1 (cos θ, sin θ, α) .
The function G m+1 (θ) having only total powers of (m + 1) in sin θ and cos θ and thus G m+1 (θ + π) = ± G m+1 (θ) for odd and even m respectively. The zeros of G m+1 (θ) determine the fixed points at infinity and now it is evident if θ j is a zero of G m+1 (θ) then so θ j + π. For more details about Poincaré Sphere and capturing the behavior at infinity one can consult [25,26].

Normal Forms and Simplifications
It is the time now to introduce the normal form technique which enables us to simplify the equations describing the dynamical system. In this Subsection we follow closely the notation found in [26]. In order to understand what we mean by a simplification here, it is important to separate the equations describing the dynamical system into linear and nonlinear part as, where J, determining the linear part of the system, is simplified into one of the Jordon canonical forms.
As to the nonlinear part, it is organized as, where F i (x) means terms of order x i . Starting with simplifying the second order term by introducing the nonlinear transformation, where h 2 (y) is of order y 2 , when applied to Eq.(25) leads to, Dealing only with second order terms amounts to, To eliminate the second order term, one need to impose To be more concrete we introduce H 2 , the space of homogenous two column polynomails of degree 2, and the map L Using the map L J the space H 2 can be nonuniquely decomposed, direct sum composition, as where G 2 represents the space complementary to L (2) J (H 2 ). Thus the simplification takes place by eliminating F 2 , if it is in the range of L (2) J , through choosing a suitable h 2 (y) leaving terms belonging to G 2 .
Applying the technique of the normal form to the case of interest where J and H 2 are respectively given as, and the parameter α is kept without normalization for the sake of clarity and simplicity. The resulting L J (H 2 ) = Span The construction of G 2 is a little bit more involved as we have to find the orthogonal complement of L J (H 2 ). The determining properties are; where the bracket · · · | · · · indicates the Euclidean inner product. The vanishing of V L J | X for any X ∈ H 2 leads to the vanishing of V L . The easier way to get V is to construct a 6 × 6 matrix representation for L (2) J where considering the vector space corresponding to H 2 as The resulting matrix representation of L J is found to be, and the resulting zero eigen-space for L (2)T J and hence G 2 are found to be spanned by It is clear that L J (H 2 ) and G 2 , as given respectively in Eq.(35) and Eq.(39), are orthogonal but this is not necessary in direct sum composition introduced in Eq.(32). One can combines J (H 2 ) with elements in G 2 , found in Eq.(39), to find additional two realization for G 2 . Last, the two-dimensional dynamical systems characterized by J, in Eq.(33), in their simplest possible form containing quadratic terms are, where a and b are two independent constants. The processes of simplification using normal forms can be continued to the terms of O y 3 and that is the maximum we need in our present work. All procedures followed previously for simplifying second order terms can be straight forwardly applied to third order terms. The dynamical system, after simplifying second order terms, isẏ where F r 2 (y) are the simplified O y 2 terms whileF 3 (y) are the O y 3 terms in their unsimplified forms. The simplification ofF 3 (y) terms is achieved by making the following transformation, where for the notational simplicity we use the same name for y for new and old variables describing the dynamical system. The resulting necessary condition to simplify O y 3 terms is, One can define analogous to L The composition of H 3 as a direct sum of L while G 3 which is orthogonal to L As we know that G 3 is not necessarily to be orthogonal to L with G 3 to get other two alternatives for G 3 which are namely,

Analysis of Universe Filled with Perfect Fluid
It is tempting to apply the dynamical system theory to the system of Eqs. (12) in its full generality, but it might be better to first study special cases in order to get deep insight into the dynamical system represented by these equations. The first simple case is to set cosmological constant and viscosity coefficient to zero α = ω,Λ = 0,ξ = 0 T . Thus, the resulting equations arė Unless ω = −1 nor ω = − 1 3 , the system has only one finite fixed point at x = x 0 = (0, 0). Then the Jacobian matrix evaluated at the fixed point turns out to be, This clearly shows that the Jacobian matrix has a double zero eigen values, λ 1 = λ 2 = 0, while the corresponding generalized eigenvectors are determined to be e 1 = (1, 0) T and e 2 = (0, 1) T . Such a kind of system, where there are two zero eigenvalues, is termed as a Bogdanov-Taken system. The stability of such a system can't be decided according to the linear stability theory. Now let us turn to the case for ω = −1 where we find an infinite number of fixed points along the curve, x 2 = 3 x 2 1 , and the resulting Jacobian is, The eigenvalues resulting from this Jacobian are λ 1 = −2 x 1 and λ 2 = 0 while their corresponding eigenvectors are respectively e 1 = (1, 0) T and e 2 = (1, 6 x 1 ) T . The direction e 1 is a stable when (x 1 > 0) and unstable for (x 1 < 0). The other direction e 2 is along the tangent of the parabola curve (x 2 = 3 x 2 1 ) where all points along the parabola are fixed points.
The last remaining special case for ω = − 1 3 , where we find an infinite number of fixed points, this time, along the x 2 axis and leading to the following Jacobian, The Jacobian matrix has a double zero eigenvalues, λ 1 = λ 2 = 0, while the corresponding generalized eigenvectors are determined to be e 1 = (1, 0) T and e 2 = (0, 1) T . Once again, the occurrence of the double zero eigenvalues makes the stability analysis not possible according to the linear stability theory.
The fixed points at infinity and as explained in Section 2.2 can be determined by the zeros of the function G m+1 (θ), defined in Eq. (24), which for Eqs.(48) amounts to For ω = − 2 3 , all points at the circle of infinity are fixed points otherwise there are finite number of fixed points corresponding to θ = 0, π 2 , π, 3 π 2 . Considering the flow only along the circle at infinity and provided that (2 + 3 ω) > 0, the points (θ = 0) and (θ = π) can be shown to be respectively stable and unstable while the points (θ = π 2 ) and (θ = 3 π 2 ) are found to behave as saddle but of non-hyperbolic type since dG 3 (θ) dθ is vanishing at θ = π 2 or 3 π 2 . Having (2 + 3 ω) < 0, all directions of flow are reversed on the circle at infinity leading to switching fixed point from stable to unstable and vice versa. The saddle points keep their type unchanged.
As to the normal forms, the system in Eqs.(48) when compared to the form in Eq. (25), J has the form of Eq.(33) with α = − 1 6 (1 + 3 ω), F (x) turns out to be It is evident that F (x) contains two pieces the first one belongs to G 2 , see Eq.(35), while the second one to L (2) J ,see Eq.(35). Thus the piece belonging to L (2) J can be shown to be canceled by the following transformation, The resulting equations in terms of y i s variables turn out to be, One can get another alternative normal form aṡ which can be achieved by the following transformation, As is clear the reduction to normal forms produces terms of O y 3 which ,in our case, belongs to In this case, one can also prove analytically that there is always a solution along the H axis representing an empty expanding universe (H > 0) or contracting one (H < 0) with negative curvature (Milne universe). This feature can be easily observed in the phase diagrams as depicted in Fig.(2) and Fig.(3). In fact the presence of that particular solution, i.e. Milne universe, serves as a phantom divide separating zone (ρ + p = 0), which can never be crossed. One should remember the fact that the solution curves in phase diagrams can't intersect each other except at a fixed point.
The above case of a perfect fluid contains a collection of interesting cosmologies that includes different types of bounce cosmologies. For example, in the cases presented in Fig.2(A,a), if we started with an expanding universe at some point in time i.e., H > 0 (where, ρ > 0) the Hubble rate will keep decreasing till it vanishes, then becomes negative describing a collapsing universe. This case has a maximum scale factor a max and a minimum density ρ, in addition, the whole evolution occurs in a finite time since it does not contains any fixed points. Furthermore, the values for H and ρ > 0 are bounded from below but not bounded from above. But the cases presented in Fig.2(C,c) and Fig.2(D,d) describe cosmological bounces with values of H and ρ > 0 which are bounded from below and from above. These last cases have one fixed point, therefore, the whole evolution time is infinite.   2 . x 1 and x 2 respectively denote the dimensionlessH andρ as defined in Eq. (7). X and Y are the coordinates on the Poincaré sphere as defined in Eq. (19). The dotted circles represent fixed points.

Analysis of Universe Filled with Perfect Fluid in the Presence of Cosmological Constant
The second simple case is to ignore viscosity in Eq. (12) and thus the dynamical system reduces to, The fixed points are determined to be three fixed points. The first two fixed points together with their Jacobians are, The reality of fixed points necessitates thatΛ ≥ 0 and hence the real eigenvalues for the Jacobian in Eq.(59) together with their corresponding eigenvectors are, where the sign (±) indicates to fixed points having x 1 = ± Λ 3 . The types of fixed points are controlled by ω as follows; the fixed point x 1 = + Λ 3 , x 2 = 0 is a stable (sink) one for ω > −1 and a saddle otherwise while the fixed point x 1 = − Λ 3 , x 2 = 0 is a unstable (source) one for ω > −1 and a saddle otherwise. Here a typical behaviour of saddle-node bifurcation is observed, where forΛ < 0 there is no fixed point but atΛ = 0 a single fixed point appears at the origin and then forΛ > 0 two fixed point appear along the H(or x 1 ) axis. The stability of the two appearing fixed points depends the value of ω as just discussed previously. This finding concerning the saddle-node bifurcation can be conveniently depicted in the following diagram, Fig.(4), consisting of two parts depending on the value of ω. The third fixed point together with its Jacobian matrix are, The reality of this fixed point is ensured for all real values ofΛ and ω while the reality is not guaranteed for the eigenvalues of the associated Jacobian matrix. For this case the eigenvalues together with their eigenvectors are, The fixed point is of a saddle type forΛ (1 + ω) > 0 while of a center type forΛ (1 + ω) < 0. This persistent fixed point a long the x 2 axis changes its type from saddle to center according the sign of Λ (1 + ω) and this is also a typical behavior bifurcation called degenerate Hopf bifurcation. The bifurcation behavior can be neatly and conveniently depicted in the following diagram consisting of three parts depending on the value of ω. 1+3 ω , determining the fixed points along the x 2 axis as a function ofΛ for a fixed value of ω in the range specified. The half-filled circle and arrowed circle represent a saddle and a center respectively.
It is worthy to stress that the flow depicted by bifurcation diagrams in Fig. (4) is restricted to the flow along the x 1 axis while the proper flow should be inferred from a kind of graphs as provided in Fig. (6) and Fig. (7) where the true flow is a two dimensional one. Needless to mention that the flow depicted in Fig. (5) should be viewed in the proper context of two dimensional flow in the x 1 − x 2 plane where a fixed point as a center along the x 2 can have a meaning. In fact, this kind of reduction is intended for simplification and more clarification otherwise one should work in a plane describing the parameter space for ω andΛ divided into regions according to the behavior of the emerging fixed points. One should not take this kind of reduction too latterly and keep in mind that the whole picture that these emerging fixed point whatever saddle, stable, unstable and center are coexisting together as shown in various figures like Fig. (6) and Fig. (7). This kind of reduction proves to be more useful and convenient when viscosity is included where the parameter space would be a three dimensional one leading to a difficulty in visualization. Another remark, in both bifurcations diagrams in Fig. (4) and Fig. (5), the nature of the fixed point whenΛ = 0, namely the origin except at ω = −1 where there an infinite number of fixed points, should be inferred from the graphs in Fig. (2) and Fig. (3).
A careful treatment is required for the special case where ω = −1 which leads to, There is a family of fixed points determined by the relation x 2 1(0) = 1 3 x 2(0) +Λ . The reality of these fixed points is ensured by requiring x 2(0) +Λ ≥ 0. The fixed points and their associated Jacobian matrices are, The eigenvalues for the Jacobian in Eq.(64) together with their corresponding eigenvectors turn out to be, where the sign (±) denotes fixed points having . The direction e 1 is a stable when (x 1 > 0) and unstable for (x 1 < 0). The other direction e 2 is along the tangent of the parabola curve The other special case for ω = − 1 3 also requires a careful treatment and here is the equations governing this case as obtained from Eqs.(58) after substituting ω = − 1 3 , There are only two fixed points that are given as x 1 = ± Λ 3 , x 2 = 0 as opposed to case, in the absence of cosmological constant, where there an infinite number of fixed points along the ρ axis. Thus, the issue of the presence of an infinite number of fixed points is cured for that case of ω = − 1 3 after including cosmological constant.
In order to get real fixed points one should imposeΛ ≥ 0. The fixed points and their associated Jacobian matrices are, As is clear the system has degenerate eigenvalues ∓ 2 1+3 ω . ForΛ (1 + ω) > 0. The stability of the two fixed points along the x 1 axis are, the right one is stable (sink) while the left one is unstable (source). In contrast, for Λ (1 + ω) < 0, the two fixed points along the x 1 axis are of saddle type. Now, the third fixed point along x 2 axis, it is a saddle forΛ (1 + ω) > 0 and a center otherwise. The rest of figures in (7)(A,a,C,c,D,d) confirm the analytical analysis revealing that whenΛ < 0 and ω = − 1 3 , there is no fixed points a long the x 1 axis but only a one along the x 2 axis being a saddle forΛ (1 + ω) > 0 and a center forΛ (1 + ω) < 0. and ω = − 3 2 ,Λ = 1 . x 1 and x 2 respectively denote the dimensionlessH andρ as defined in Eq.(7). X and Y are the coordinates on the Poincaré sphere as defined in Eq. (19).
The dotted circles represent fixed points. and ω = − 2 3 ,Λ = − 1 2 . x 1 and x 2 respectively denote the dimensionlessH andρ as defined in Eq.(7). X and Y are the coordinates on the Poincaré sphere as defined in Eq. (19).
The dotted circles represent fixed points.
Remarks concerning the fixed points at infinity and the normal forms are in order. First, we find the same fixed points as the case without including the cosmological constant and the fixed points are determined by the same function found in Eq.(52). This is can be easily understood since the introduction of the cosmological constant adds only zero order terms and thus doesn't affect the behavior at infinity compared to the other present higher order ones. All figures in Figs. (6) and Figs. (7) for the compact phase portraits confirms this finding concerning the fixed points at infinity. A particular emphasis for the case, ω = − 2 3 , where the circle at infinity in its totality are fixed points as clear from Fig. 7(c). Second, as to the normal form one can use the following transformation, then the system in Eq.(58) will reduces to, The normal form corresponding to the case where ω = − 1 3 needs careful treatment since the transformation in Eq.(68) is singular. Introducing the variables z 1 = x 1 − Λ 3 and z 2 = x 2 then Eq.(66) would transform into,ż In this new form described by Eq.(70), the Jacobian, J, is clearly proportional to the identity and thus L J (H 2 ) = H 2 which enables us to remove any quadratic terms. Removal of quadratic terms is not for free but at the expense of introducing higher order terms. As an example one can try the following transformation that has a validity not at the whole region of the coordinates but at small neighborhood around the origin whose size is depending onΛ, where F = Λ 3 − 2 Λ 3 y 1 . The above transformation when applied to Eq.(70) results in the following, The dots in Eq.(71) and Eq.(72) indicates the neglected higher order terms. It is important to stress that there two extreme cases for the Jacobian where it is zero or proportional to the identity. In both cases the simplification introduced through normal forms losses its appealing and the reason behind is detailed as follows; For J = 0 we have L J (H 2 ) = 0 implying that any F 2 (second order terms) can't be transformed away, while for J proportional to the identity we have L (2) J (H 2 ) = H 2 which means that we can remove any second order terms but at the expense of introducing other higher order terms as obtained in Eq.(72).
The case of a perfect fluid with cosmological constant contains new interesting features in addition to bounce cosmologies, which is the appearance of a pair of fixed points a long the H-axis. This pair admits new type of cosmological models in which the universe is interpolating between two fixed points one in the negative H region and another in the positive H region. As presented in Fig.6(A,a,C,c), the universe could start with a fixed point along the negative H-axis and end up with another fixed point a long the positive H-axis passing through a bounce, i.e., H = 0 point. Another new feature here is the existence of oscillating cosmological solutions as shown in Fig.7(C,c). In this intersecting case, all solution are either bounces or oscillating cosmologies with finite evolution time and a minimum density ρ in the case of bounce or minimum and maximum values for both H and ρ in the case of oscillating cosmologies.
As a final remark we assert the existence of Milne type solution in the presence ofΛ as can be verified by direct integration of Eq.(58) when x 2 = 0. That particular type solution prohibits the crossing from solutions with positive ρ to ones with negative. In another equivalent way, there is no phantom divide crossing as shown in the previous section when the cosmological constantΛ is absent. There could be fixed points along the Milne type solution but to cross through these fixed points requires infinite time to reach them.

Analysis of Universe Filled with Bulk Viscous Fluid in the Presence of Cosmological Constant
The cosmological equations in their full generality, in the presence of ω,Λ andξ, arė These equations can be transformed into one of the standard normal form given as, This form corresponds to the normal form for a degenerate Bogdanov-Taken bifurcation as classified in [24,27]. The form can be achieved by the following transformation, given here with its inverse, (75) After performing the previous transformation, the parameters α 1 , α 2 , α 3 , b, d and e are found to be According to the classification and the study carried out in [24,27], the parameters b and d together with their combination d 2 + 8b shouldn't be vanishing. Actually, b is vanishing for ω = −1 and d for ω = − 5 3 while d 2 + 8b for ω = − 1 3 . Although the form in Eq.(74) is relevant to recognize the classification of the system described in Eq.(73) as a degenerate Bogdanov-Taken bifurcation, it turns out more simpler and transparent to study the bifurcation of the system in its original form in Eq.(73).
The starting point for this bifurcation study is to find and classify the fixed points and then investigating their behavior under changing parameters (ω,Λ,ξ). The first possibility is where we have a fixed point along the x 2 axis which is given together with its Jacobian as, . (77) The eigenvalues for the Jacobian in Eq.(77) together with the corresponding eigenvectors are, This fixed point is always present provided ω = − 1 3 . For nonvanishingξ > 0 and where ∆ 1 < 0 the fixed point is a repelling center. When ∆ 1 = 0, the fixed point turns out to be unstable (source) and continues to be unstable (source) whenever ∆ 1 < 9ξ 2 . When ∆ 1 ≥ 9ξ 2 , the eigenvalue λ 2 vanish at ∆ 1 = 9ξ 2 and then start to be negative leading to a saddle fixed point. To simplify matter for depicting the behavior of the fixed point as the parameters change, we fix ω at a specific values and then the condition ∆ 1 = 0 turns out to define a parabola in the (Λ,ξ) plane given byΛ = −9ξ 2 4(1+ω) . This parabola together withξ axis divide the (Λ,ξ) plane into four distinct regions and each region has a characteristic behavior for the fixed point. The property of the fixed point is changing with parameter according to Hophf bifurcation. This typical kind of bifurcation is shown in a bifurcation diagram in Fig.(8). The phase space diagrams, compact and uncompact ones, are also displayed for representative cases as in Figures Figs. (9-12). For fixed value ofξ = 0.1, we choose the other parameters (Λ, ω) in such a way to have only, whenever possible, a fixed point along the x 2 axis with a clear appearance as done in Figures Fig. (9) and Fig. (10). As to the fixed points not appearing along x 2 , but along the flat curve solution, we anticipate the results which will be explained later in this section. The figures in Fig.(11) and Fig.(12) are devoted to the case of stiff matter, (ω = 1) and thus (1 + ω) > 0, with varyingΛ to produce all possible scenarios for the fixed point along the x 2 axis. The figures in Fig.(13) are devoted to the case of phantom matter, (ω = −2) and thus (1 + ω) < 0, with varyingΛ to produce all possible scenarios for the fixed point along the x 2 axis. It is worth mentioning that in these set of figures we include the case for (ω = −1) where fixed point doesn't occur except at the x 2 axis at x 2 = −Λ and having a saddle character.
To have more quantitative results, we present in Table (    Representative cases are ω = −2 whileΛ = 2.5, 0.1575, 0.011, 0 and −0.1. x 1 and x 2 respectively denote the dimensionless H andρ as defined in Eq.(7). X and Y are the coordinates on the Poincaré sphere as defined in Eq. (19). The dotted circles represent fixed points. The second possibility where we have two fixed points that are given as, where The reality of the fixed points are ensured when ∆ 2 ≥ 0 (or equivalentlyΛ ≥ − 3ξ 2 (1+ω) 2 ). The real fixed points (x 1+ , x 2+ ) and (x 1− , x 2− ), when realized, are always located on the parabola describing flat solution x 2 1 = 1 3 x 2 +Λ . The Jacobian at these fixed points, Eq.(80), are found to be where the sign (±) respectively denotes the fixed points (x 1+ , x 2+ ) and (x 1− , x 2− ). The resultant eigenvalues and their associated eigenvectors are The case with two fixed points, along the flat curve solution, is more involved than the case of a single point along the x 2 axis. In the parameter space where ∆ 2 < 0 there are no fixed points at all. When ∆ 2 = 0 an emergent single fixed point (non-hyperbolic one) appears whose coordinates, associated eigenvectors and eigenvalues are, after using eq.(80) and eq.(83), The eigenvector e 1 corresponding to the zero eigenvalue is in the same direction as that of the tangent of the flat curve solution at x 1 =ξ (1+ω) , x 2 = 6ξ 2 (1+ω) 2 whenever (1 + ω) > 0 and opposite otherwise. As to the direction given by e 2 , it represents a stable direction whenever (1 + ω) > 0 and an unstable for (1 + ω) < 0.
In the parameter space where ∆ 2 > 0, the single fixed point at ∆ 2 = 0 is shattered into two fixed points as described by Eq.(80) and Eq.(83). The fixed point designated by (+) is a stable (sink) fixed point when (1 + ω) > 0 and of a saddle type for (1 + ω) < 0. The other fixed point designated by (−) doesn't behave in a simple manner as one designated by (+). When ∆ 2 = 9ξ 2 that leads toΛ = 0 provided that ω = −1, the fixed point turns out to be at the origin (x 1− = 0, x 2− = 0) and the associated eigenvalues and eigenvectors are, The direction e 1 corresponds to a stable direction while e 2 has a zero eigenvalue which means that fixed point is a non-hyperbolic one. Apart from this value of ∆ 2 and as 0 < ∆ 2 < 9ξ 2 the fixed point, (−), is an unstable (source) for (1 + ω) < 0 while of saddle type for (1 + ω) > 0. The behavior is switched off for ∆ 2 > 9ξ 2 which means getting unstable (source) fixed point for (1 + ω) > 0 while a saddle type for (1 + ω) < 0. The corresponding bifurcation diagram can be simplified by considering a fixed value for ω and depicting the condition ∆ 2 = 0 as a parabola curve in the plane (Λ,ξ) given byΛ = − 3ξ 2 (1+ω) 2 . This parabola divide the plane (Λ,ξ) into three distinct regions and each region has a characteristic behavior for the fixed points. All these behaviors are displayed in the bifurcation diagram in Fig.(14) showing a similar behavior to that of saddle-node bifurcation. The phase space diagrams, compact and uncompact ones, are also displayed for representative cases as in Figures Figs. (15) and (15). The finding for these representative cases are summarized in Table (2) with the same notations used in Table (1   The fixed points, at infinity, is determined by the zeros of the function G m+1 (θ), defined in Eq.(24), which for Eqs. (12) amounts to G m+1 (θ) m=2 = G 3 (θ) = − cos 2 θ sin θ (2 + 3 ω) − 18ξ cos θ .
For nonvanishing value ofξ, we have four fixed points corresponding to θ = π 2 , 3 π 2 , ϕ, ϕ + π where ϕ = tan −1 18ξ 2+3 ω . The non vanishing value ofξ prevents the occurrence of an infinite number of fixed points at infinity when ω = − 2 3 and reducing them to just pair of fixed points at θ = π 2 , 3 π 2 . Considering the flow only along the circle at infinity, the two fixed points at (θ = π 2 ) and (θ = 3 π 2 ) are behaving as saddles but of nonhyperbolic type, for any values of the relevant parameters, since dG 3 (θ) dθ is vanishing at θ = π 2 or 3 π 2 . Again by considering the flow along the circle at infinity, the other two fixed points corresponding to θ = {ϕ, ϕ + π} can be shown to be of opposite type such that one is stable and the other is unstable depending on the sign of (2 + 3 ω) and which quadrant the angle ϕ belongs to.
As we have seen the presence of viscosity prevents the occurrence of an infinite number fixed points wherever they are; at the finite domain or the circle at infinity. Moreover, the occurrence of periodic orbits are prohibited by the presence of viscosity. In fact, the absence of these two kinds of behaviors is crucial since it is among the basic requirement for the dynamical system to have structural stability according to the criteria presented in [26,25]. In fact, Peixoto theorem for a flow defined on a compact two-dimensional as in [26,25], which in our case the flow induced on the Poincare sphere, can be used to decide the presence of structural stability or not in the considered cosmological models. According to Peixoto theorem [25,26], the hyperbolcity of the fixed points is a necessary conditions to attain structural stability which can't be satisfied in our case since we have always non-hyperbolic fixed points at infinity corresponding to (θ = π 2 ) and (θ = 3 π 2 ). In this case when bulk viscosity is included, the curve ρ + p = 0 (phantom divide curve) which turns out to be a straight line given byρ (1 + ω) − 6ξH = 0 is not a solution curve as can be checked explicitly. Thus, there could be a solution curve that might cross the phantom divide curve in a finite time. The crossing of phantom divide can be noticed, as for examples, from the phase portraits presented in Fig.9(A,a) and Fig.10(C,c).
The cosmological model incorporating bulk viscosity in one of its simplest form can still lead to some interesting consequences that could be relevant to the actual physical universe. We find that our parameters (ω,Λ,ξ) could be adjusted to have three fixed points one along the x 2 (ρ)-axis and the other two along the flat curve solution. The one along the x 2 axis is a repulsive center and it represents a static universe. It is implausible to consider our physical universe had started in the neighborhood of this repulsive center since it contradicts with the standard scenario of initial big-bang and early inflation. Thus we are left with the two fixed points along the flat curve solution. It is convenient to restrict the dynamical study to the spatially flat case [4] which can be derived from the first equation in the set of Eqs.(73) and found to be, The fixed points corresponding to this flow are, The flows as depicted in Fig. 17(A) single out a one starting from initial big-bang and ending at a de-Sitter universe represented by the fixed point x 1 + as the only possible candidate describing our physical universe. This scenario of starting with big-bang and ending up with late acceleration could be achieved in the presence of viscosity without the need for including cosmological constant. This finding is consistent with what have been found in earlier study [11]. Also, the big-bang can occur in both closed and open universe and still ending up with late acceleration (fixed point x 1+ ) as evident from the plots in Fig.10(C,  c). The other two remaining possibilities are not good candidates for describing our physically observed universe as explained as follows. The first one staring from x 1− and ending up with x 1+ (starting with small value for Hubble parameter and ending with a larger one) can't describe the actual universe since the opposite behavior is required. As to the second one starting with x 1− and going to x 1 = −∞ which means passing through contracting phase and ending with a big crunch and this is clearly doesn't match the behavior of the observed universe which, at present, is expanding with acceleration. The bulk viscosity coefficient can be taken as a varying function depending on x 2 and for simplicity we assume linear dependence asξ (x 2 ) = α x 2 . Inserting this form of varyingξ into the cosmological equations in Eq.(73) would give fixed point as,

Discussion and conclusion
In this work, we present a complete dynamical study for a bulk viscous cosmology with a single fluid in the presence of a cosmological constant. For the sake of illustration and clarification we don't study, in a single step, the bulk viscous cosmological model in its full generality containing the three parameters namely ω,Λ andξ, but our investigation is carried out in three different stages. The first stage, we consider only ω to be non-vanishing and thenΛ is included while finallyξ is introduced. In each of these stages, the fixed points, whether they are at the finite domain of the phase space or at infinity, are studied and classified. Also, the normal forms are obtained for each stage together with phase space portraits for meaningful representative cases. Suitable and convenient bifurcations diagrams are plotted for illustrating the changing behavior of fixed points as the relevant parameters vary.
The dynamical system corresponding to bulk viscous cosmological model is shown, in Section 2, to be a two dimensional unregenerate Bogdanov-Takens system following the classification carried out in [24]. This point concerning the classification is a novel result up to the best of our knowledge. Another issue besides the classification which is worthy to be discussed is the structural stability which means that the qualitative behavior of the system is unaffected by small perturbations. In two dimensional dynamical system, simple criteria can be established for testing structural stability utilizing Peixoto theorem for a flow defined on a compact two-dimensional space as in [26,25]. In our study the flow induced on the Poincaré sphere can be used to shed some light on the structural stability of the considered cosmological models. One of the basic criteria is to have a finite number of fixed points and periodic orbits which are hyperbolic. Here the finiteness of the number of fixed points, as shown in section 6, can be achieved by introducing a non vanishing viscosity while the hyperbolicity of fixed points in the finite domain of phase space can be attained by restricting the relevant parameters (ω,Λ,ξ), as an example,Λ > 0 and ω +1 > 0. Unfortunately, as can be inferred from Eq.(86), we have at infinity fixed points, (θ = π 2 , 3 π 2 ), that are always nonhyberbolic for any choice of the parameters. Thus the structural stability for the considered cosmological models can't be achieved even after introducing viscosity and for any chosen region in the parameter space (ω,Λ,ξ).
The late acceleration behaviour is a confirmed feature of the observed universe due to the observations of distant supernova type Ia [13,14] and cosmic microwave background anisotropy measurements [15,16]. The bulk viscous fluid can provide us with a source of this late acceleration, even in the absence of cosmological constantΛ, as discussed in Section 6 and illustrated in Fig.17(A). In this case, the cosmological model interpolates between big-bang and late acceleration. This induced late acceleration can be attributed to the effect of negative pressure associated with viscosity as is clear from the expression of T µν in Eq. (4).
The bulk viscous fluid with viscosity coefficient dependent on density asξ(x 2 ) = α x 2 and in conjunction with cosmological constant can provide us with a non singular cosmological model. The model interpolates between an inflation point, x 1 = 1+ω 6α , and a late acceleration point, x 1 = Λ 3 , as discussed in Section 6 and illustrated in Fig.7(B). The complete study of this model including viscosity coefficientξ dependent on x 2 along the lines presented for the one of constantξ would be a subject for the future work. To confront the introduced cosmological models with observational data like Type Ia supernova, one should include an additional fluid component that represents matter besides the dark energy component represented byΛ and a viscous fluid. The introduced parameterΛ,ξ and α might enhance the agreement with observational data but that needs a detailed study which would be a subject for a future work.