Bifurcations in Ratra-Peebles quintessence models and their extensions

We have used the dynamical system approach in order to investigate the dynamics of cosmological models of the flat Universe with a non-minimally coupled canonical and phantom scalar field and the Ratra-Peebles potential. Applying methods of the bifurcation theory we have found three cases for which the Universe undergoes a generic evolution emerging from either the deSitter or the static Universe state and finishing at the deSitter state, without the presence of the initial singularity. This generic class of solutions explains both the inflation and the late-time acceleration of the Universe. In this class inflation is an endogenous effect of dynamics itself.


Introduction
The main aim of cosmology is to study the structure and the evolution of the Universe at the large scale. In this context a crucial role plays the notion of a cosmological model obtained from general relativity. If we assume a spacetime with the maximal symmetry of the spacelike section of constant time then the evolution of the Universe can be analysed with tools of dynamical systems methods. In principle, there are two advantages of such a treatment of the dynamics. Firstly, evolutional paths are represented in a geometrical way in the phase space-the space of all states of the system under consideration. Secondly, it gives many possibilities of visualisation how solutions (represented in the phase space by trajectories) depend on the initial conditions. Are they exceptional or typical in the phase space? This question is strictly related to a problem of the stability of solutions.
Right-hand sides of the dynamical system depend not only on the state variables (which form the vector of state) but also on parameters. These theoretical parameters are present in any effective theory of the Universe. We believe that their interpretation will emerge if we find more fundamental theory giving us better insight into their nature. Of course, theoretical investigations should be complemented with a statistical estimation of model parameters from the astronomical and astrophysical data.
From the theoretical point of view, it is important to know how the structure of the phase space changes under variation of model parameters. For such an investigation the bifurcation theory is dedicated. Its methods allows to obtain some critical (bifurcation) values of model parameters for which a topological structure of the phase space (modulo homeomorphism preserving the direction of time along trajectories) is not equivalent in the topological sense. If the dynamical system, decribing the evolution of the physical system, undergoes the bifurcation at some value of its parameter, one can classify all possible evolutional paths for all theoretically admissible values of parameters. Hence, in principle, two qualitatively different types can be distinguished. The first case is a the model with the parameter at the bifurcation value which decribes some state, epoch, of the Universe evolution. In this sense the model is fine-tuned as some exact value of the model parameter is required to describe the Universe dynamics. The second case is a model with other, non-bifurcation values of its parameter. The dynamics of the model is qualitatively indifferent with respect to the choice of the value of the model parameter.
Dynamical systems methods are a powerful tool to analyse the dynamics of physical systems and the application of these methods to analyse the cosmological models is very important as the stability conditions allow to constrain or even rule out some models [4]. The dynamical analysis can be enriched further with using the biffurcation theory methods. Recently, the problem of bifurcations in the FRW cosmology with perfect fluid and the cosmological constant was investigated in the work of Kohli and Haslam [14]. They demonstrated that as the cosmological constant parameter is varied, expanding de Sitter and contracting de Sitter universes emerge via fold bifurcation. This type of bifurcation occurs in a neighbourhood of the Minkowski spacetime in the phase space. Also, bifurcations play important role in the idea of the Jungle Universe. Following this approach to study the dynamics of FRW models filled with different barotropic fluids, the jungle cosmological dynamical system assumes the form of a generalised Lotka-Volterra system where the competitive species are replaced with fluids [19]. The problem of bifurcation is also strictly related to the problem of structural stability which denotes changes of the phase space structure under small change of right-hand sides of the system [15]. It is also interesting to discover the phenomenon of period-doubling bifurcation in strongly anisotropic Bianchi I quantum cosmology [3].
Basic description of the current acceleration of the Universe is provided by the cosmological constant Λ, whose natural interpretation is the energy of the quantum vacuum. However, this energy could not be constant during the evolution of the Universe from early epochs to the current state. Hence, in order to include both inflation and late time acceleration in one model, it is necessary to introduce for example a tool making the vacuum energy dependent on time. Since the Einstein-Hilbert action must be covariant, one cannot just put Λ(t). The solution is to lead in a scalar field φ together with potential V (φ), driving the evolution of φ. Such a postulate has been suggested for the first time in [18], assuming the power-law form of potential, V (φ) ∝ φ n , where n = const. The scalar field approach for describing the late time acceleration of the Universe is known under a name of quintessence [23] (for the alternative approach to quintessence see [13]).
Apart from standard general relativity terms, it is considered the additional term ξφ 2 R in the Einstein-Hilbert action, describing the coupling between the gravity (represented by the Ricci scalar R) and the scalar field φ, where ξ is the coupling constant. For ξ = 0 this coupling is minimal, while for ξ = 0 it is non-minimal. The non-minimal coupling arises from quantum corrections to the scalar field theory [1,5,12] and is necessary for the renormalisation procedure of the scalar field in the curved space [7,10]. Moreover, a non-minimal coupling of the scalar field to the gravitation is widely used for modelling inflation [9,17]. The generic cosmological evolution in the scalar field cosmology with non-minimal coupling, which does not possess a singularity was investigated by Hrycyna [11]. These types of evolution correspond to trajectories starting from an unstable de Sitter state and going to a stable de Sitter state.
Observational constraints on scalar field cosmological models have been recently investigated in [2]. Using the Bayesian methods different models of the scalar field cosmology were analysed. The evidence for which as approximation the information criterion BIC was used, indicated that the Ratra-Peebles potential were favoured.
We present the analysis of the dynamical system, describing the evolution of the universe 1 with the dark energy as a scalar field with a potential. The application of bifurcation theory methods enables us to fully understand the dynamics of the system, which changes its qualitative properties under the variation of its parameters. The dynamical system under consideration includes both the canonical and the phantom scalar field models, which are distinguished by a value of the discrete parameter ε. We assume the power-law form of the potential with the exponent n as a model parameter and the constant non-minimal coupling between the gravity and the scalar field, described by the non-zero parameter ξ.
For the purpose of our investigation, we use elements of local bifurcation theory, whose basics we depict in Section 2. We discuss local bifurcations of codimension 1, as they appear in models under consideration. For getting more information about the bifurcation theory we refer readers to [6,8,16,20,22].
In Section 2 some short introduction to bifurcation theory is given. In Section 3 we derive the dynamical system from cosmological equations and the potential function, in terms of dimensionless phase space variables. In Section 4 we find equilibria of the dynamical system, analyse their stability properties, designate conditions they represent the de Sitter universe, find, on this basis, possible scenarios of the de Sitter-de Sitter evolution-which avoids the initial singularity and explains the inflation, as well as the late time acceleration-then prepare the full analysis of bifurcations of the local stability of equilibria for these scenarios. Next, Section 5 is devoted to the preparation the phase portraits of the system for non-singular scenarios on the Poincaré sphere. Finally, in Section 6, we plot the evolution of some physical parameters over the time, within non-singular evolutionary scenarios.

Elements of applied bifurcation theory
Let us start with some basics on the applied bifurcation theory. First, we give necessary, from the point of bifurcation theory's view, notions of dynamical system theory. Then, we discuss definitions related with the topological equivalence of dynamical systems and bifurcations. Finally, we present a basic description of most important types of local bifurcations. We limit the discussion to bifurcations of codimension one in two-dimensional dynamical systems since this is the only kind of bifurcations that appears in our analysis. Definitions 1-9 are quoted from Kuznetsov's book [16].

Dynamical systems
In deterministic processes, future and past states of a system can be revealed by knowing its present state and a law describing the evolution. The mathematical formalisation of this fact is the notion of a dynamical system.
A set X of all possible states of a dynamical system is called a state space (or phase space). Usually, one distinguishes finite-dimensional systems defined in X = R n from those defined on manifolds. The state x t ∈ X of a system changes with time t ∈ T , where T is a number set. In case T = R 1 a system is called a continuous-time dynamical system, while for T = Z it is a discrete-time dynamical system.
An evolution law of a dynamical system determines the state x t of the system at time t, provided the initial state x 0 is known. Description of the evolution is given by an evolution operator φ t , which is a map φ t : X −→ X, transforming an initial state x 0 ∈ X into some state x t ∈ X at time t, Deterministic nature of dynamical systems is reflected in following two properties of the evolution operator. First of all, where id is the identity map on X, id x = x for all x ∈ X. The second property is which implies x t+s = φ t (φ s x 0 ). Now we are ready to give the formal definition of the dynamical system, which comes from [16].
where T is a time set, X is a state space, and φ t : X → X is a family of evolution operators parametrised by t ∈ T and satisfying properties (1) and (2).
In case of continuous-time systems, a family φ t t∈T of evolution operators is called a flow. Moreover, in case of both continuous-and discrete-time systems, an orbit Γ is an ordered subset of the state space X, Orbits in the state space of a dynamical system compose a phase portrait, which is the geometrical representation of the dynamical system.
Very often one identifies a continuous-time dynamical system with differential equations describing implicitly a law of evolution. Suppose that the state space of a system is X = R n . Then, differential equations related to the system arė whereẋ = dx/dt, the state x = (x 1 , . . . , x n ) ∈ R n , α = (α 1 , . . . , α m ) ∈ R m is the parameter vector, and f : R n × R m (x, α) → (f 1 , . . . , f n ) ∈ R n is a function (a vector field) of the C k class, with k ≥ 1. We will frequently refer to dynamical systems without underlining its dependence on the parameter vector α. Thus, differential equations will have the forṁ Let us quote the following three important definitions, taken from [16], which are related to dynamical systems. First of them is an equilibrium, which is the special case of an orbit.
An equilibrium is thus a point in which the system remains forever. For differential equations (4) in an equilibrium x 0 we have Subsequent definitions concern the notion of another special kind of orbits-cycles.
Definition 3 A cycle is a periodic orbit, namely a non-equilibrium orbit L 0 , such that each point x 0 ∈ L 0 satisfies φ t+T 0 x 0 = φ t x 0 with some T 0 > 0, for all t ∈ T .
Definition 4 A cycle of a continuous-time dynamical system, in a neighbourhood of which there are no other cycles, is called a limit cycle.

Topological equivalence and bifurcations
Before introducing the notion of bifurcation, we should clarify what it means that phase portraits of two dynamical systems have the same qualitative features. For this purpose, let us quote from [16] the definition of topological equivalence of two dynamical systems.
Definition 5 A dynamical system {T, R n , φ t } is called topologically equivalent to a dynamical system {T, R n , ψ t } if there is a homeomorphism h : R n → R n mapping orbits of the first system onto orbits of the second system, preserving the direction of time.
Dynamical systems are very often studied in a vicinity of an equilibrium since it is possible to lead useful linear approximations there. Therefore, one may elaborate the topological classification of phase portraits near equilibrium points and hence the need for the following modification of Definition 5.
Definition 6 A dynamical system {T, R n , φ t } is called locally topologically equivalent near an equilibrium x 0 to a dynamical system {T, R n , ψ t } near an equilibrium y 0 if there exists a homeomorphism h : R n → R n that is (i) defined in a small neighbourhood U ⊂ R n of x 0 ; (ii) satisfies y 0 = h(x 0 ); (iii) maps orbits of the first system in U onto orbits of the second system in V = h(U ) ⊂ R n , preserving the direction of time.
Consider the dynamical system (3) dependent on parameters α ∈ R m . As parameters vary, the phase portrait of the system also varies. This results in two possibilities: either the system remains topologically equivalent to the original one, or its topology changes.

Definition 7
The appearance of topologically nonequivalent phase portrait under variation of parameters is called a bifurcation.
Thus, a bifurcation is a change of the topological type of the system. A value of parameters at which this change happens is called a bifurcation (critical) value.
If a bifurcation occurs when we fix any small neighbourhood of an equilibrium, then it is called a local bifurcation (or a bifurcation of an equilibrium). Otherwise, when a bifurcation cannot be determined by analysing only a small vicinity of an equilibrium, it is a global bifurcation.
It is often convenient to visualise how the occurrence of bifurcations depends on parameters of the system. For this purpose, let us introduce the definition of a bifurcation diagram taken from [16].
Definition 8 A bifurcation diagram of the dynamical system is a stratification of its parameter space induced by the topological equivalence, together with representative phase portraits for each stratum.
In the simplest case, the bifurcation diagram consists of a finite number of regions in the parameter space R m , inside which the phase portrait is topologically equivalent. These regions are separated by bifurcation boundaries, which are smooth submanifolds in R m (e.g. curves, surfaces).
Definition 9 A codimension of a bifurcation is the difference between the dimension of the parameter space and the dimension of the corresponding bifurcation boundary. The equivalent and more practical definition says that a codimension is the number of independent conditions determining the bifurcation.
Bifurcations are divided into many types that differ in their course. It is possible to lead the most simple differential equations that represent a given type of bifurcation. These equations are called normal forms. In case of local bifurcations, every dynamical system having a local bifurcation of a specified type is locally topologically equivalent near an equilibrium to the normal form of this type of bifurcation.

Local bifurcations of codimension one in continuous-time systems
Let us introduce most important types of local bifurcations in continuous-time systems of at most two variables. We limit the discussion to bifurcations of codimension one. We present four types of bifurcations. Three of them are bifurcations of the saddle-node family: an elementary saddle-node, a transcritical, and a pitchfork bifurcation, while the last is a Hopf bifurcation.
Bifurcations of the saddle-node family (also called fold bifurcations) arise as a collision of equilibria, when a real eigenvalue crosses zero in a response to a change of a parameter. In the case of an elementary saddle-node bifurcation, two equilibria-one stable, the other unstable-coalesce and then disappear. In the case of a transcritical bifurcation, two equilibria-one stable, the other unstable-coalesce and then separate again, exchanging their stability properties. Finally, in the case of a pitchfork bifurcation there are two possibilities: in first of them, called supercritical, three equilibria-one unstable, surrounded by two stablecoalesce into one stable equilibrium, while in the other, called subcritical, three equilibria-one stable, surrounded by two unstable-coalesce into one unstable equilibrium. Normal forms for saddle-node family of bifurcations are where α is the parameter with the critical value equal to zero. Bifurcation diagrams of the foregoing normal forms are shown in Figure 1.
A Hopf bifurcation (known also as an Andronov-Hopf bifurcation) arises, when real parts of two complex and conjugate eigenvalues crosses zero. It presents a collision of a stable or unstable equilibrium with a cycle, which shrinks to a point when the collision occurs. The normal form for this kind of bifurcation iṡ which in polar coordinates (ρ, ϕ) isρ where α is the parameter with the critical value equal to zero.

Dynamical equations
The scalar field approach describes the dark energy by the scalar field φ, which is affected by the potential U (φ) > 0. The total action is where is the Einstein-Hilbert action, and with ε = ±1 corresponding to the canonical and the phantom scalar field respectively, κ 2 = 8πG and c = = 1. The parameter ξ is called the coupling parameter; if ξ = 0, the scalar field is non-minimally coupled to the gravity. Otherwise, the scalar field is minimally coupled to the gravity. The variation of the total action (11) with respect to the metric tensor g µν yields the field equation where the energy-momentum tensor for the non-minimally coupled scalar field is given by while the variation with respect to the scalar field φ gives the Klein-Gordon equation where U ,φ := dU/dφ. We consider the spatially flat (k = 0) Friedmann-Lemaître-Robertson-Walker (FLRW) metric where a(t) is the scale factor. This produces the Friedmann equation in the form the acceleration equation and the Klein-Gordon equation where ρ φ is the energy density and p φ is the pressure of the scalar field, which are connected by the linear barotropic equation of state where w φ is the equation of state parameter. Let us note, that from the Friedmann equation (18) we see that the total energy density of the scalar field consists of the kinetic energy density term 1 2 εφ 2 , the potential energy density term U , and the term related to the gravity-scalar field interaction 3εξHφ(Hφ+2φ). Dividing this equation by the critical density 3H 2 κ −2 we obtain an energy conservation equation in the form where is the scalar field kinetic energy parameter and is the scalar field potential energy parameter. We assume the inverse power-law form of the potential function (known as the Ratra-Peebles potential [18,21]) where n is a dimensionless parameter, and M > 0 is a dimensional constant. The pair of parameter values ε = +1 and ξ = 0, together with the usage of the power-law form of the potential, refer to the standard Ratra-Peebles quintessence model, while other values constitute extensions of this model. Let us introduce following dimensionless real phase space variables Applying new variables, we obtain the following condition from the Friedmann equation (18) M n+4 which, in terms of energy parameters (23-24), reads where The acceleration equation (19) yieldsḢ where the equation of state parameter, in terms of new variables, is .
The left-hand side of equation (30) is equal to −q − 1, where q = −äaȧ −2 is the deceleration parameter; q < 0 for accelerated expansion (or decelerated contraction) of the universe, and q > 0 for decelerated expansion (or accelerated contraction) of the universe. Hence, we obtain which implies that for w φ < − 1 3 expansion of the universe is accelerated (contraction is decelerated), whereas for w φ > − 1 3 expansion is decelerated (contraction is accelerated). Finally, using equations (20), (27) and (30), we obtain the dynamical system, expressed by dimensionless variables (26) where f = df d ln a = H −1ḟ . Formulae for w φ (31) and for u in the system (33) diverge to infinity as In order to analyse the dynamics at this limit, we can multiply right-hand sides of the system (33) by the non-negative term 1 3 εv 2 − 2ξ(1 − 6ξ) 2 . This operation will produce a dynamical system which will be determined at the limit (34) and will have the same dynamical properties in other points as the previous system 2 Thus, we obtain the system 2 Such an operation is also related to the reparameterisation of the time. Consider the system dx/dt = f (x) and multiply its both sides by a function ξ(x) > 0. As a result of this operation, we obtain the new system dx/dτ = ξ(x)f (x), where dτ := dt/ξ(x), which has the same equilibria and is topologically equivalent to the original one.
which is determined also at the limit (34).
From initial assumptions we have U (φ) > 0 equivalent to φ > 0, which implies v > 0. It also implies that the left-hand side of equation (27) is greater than zero. It yields the following condition of physicality on variables (u, v) 4 Equilibria, non-singular evolutionary scenarios and bifurcation diagrams of local stability In this section, we will find equilibria of the dynamical system derived above, inspect stability of these points and conditions for representing the de Sitter state by them, distinguish, on this basis, evolutionary scenarios without singularity, and-finally-prepare bifurcation diagrams for these scenarios. As the phase space (u, v) is infinite, we will investigate firstly the finite region, and then go to the description at infinity. Finite equilibria of the system (35) together with existence conditions (resulting from the fact that u, v ∈ R) are shown in Table 1, while eigenvalues and values of the equation of state parameter w φ for equilibria are summarised in Table 2. Figures  2-6 present, in turn, bifurcation diagrams of the local stability of equilibrium points A, B, C and E, whose stability properties depend on values of parameters ξ, n and ε. Moreover, in these diagrams, there are put formulae for one-dimensional boundaries separating areas of the local topological equivalence.
We notice from Table 1 that, apart from the equilibrium line F , finite equilibrium points are located only on the axes of phase space. If for a finite point v = 0 but u = 0, then εΩ φ,kin → +∞. According to equation (28), provided the condition (36) is satisfied, we have Ω φ,pot → +∞. On the other hand, when u = 0 and v = 0, then Ω φ,kin = 0 and Ω φ,pot = 1 − 6εξv −2 , which is greater than zero when the condition (36) is satisfied.
The region (36) has its boundaries at v = 0 always at points A and B. The boundary of this region is the limit of the physicality of the system. If v = 0 we have H 2 → +∞ and |φ| → +∞ on this boundary. If the line F exists, the boundary has its extremum on this line at the point u = −6ξ. We denote this point as K.
In this paper, we will be investigating evolutionary scenarios of the universe, starting and for ξ = 0 *There exist only one-sided limits of w φ as v approaches the line F . These limits diverge to ±∞, depending on the values of u, ξ, n and ε.
finishing at the de Sitter state. The de Sitter universe is the solution to the Einstein field equations which assumes the dynamics of the universe to be dominated by the cosmological constant Λ, so the matter component (both baryonic and dark) is neglected. In this model, the pressure p and the energy density ρ satisfy so the equation of state parameter w is For the spatially flat universe (k = 0) the scale factor a, within the de Sitter model, depends on time as Solutions with the '−' sign in the exponent are frequently called anti-de Sitter states in contrast to the solutions with the '+' sign, called de Sitter states. The Hubble function H in this case is equal By looking at the de Sitter universe presence conditions in Table 2 we can distinguish, using bifurcation diagrams (Figures 2-6), sets of parameters for which an evolution from the initial (stable in the past) de Sitter to the final (stable) de Sitter state is possible. We can divide these sets into two groups: representing generic and non-generic de Sitter-de Sitter evolution. The generic evolution takes place, when there exists a family of solutions favouring given conditions, while the non-generic evolution takes place if only one particular solution corresponds to conditions. For example, when a family of orbits comes out of a stable in the past equilibrium (like an unstable node or a focus) representing the de Sitter state (w = −1),   and then finishes in a stable point also representing the de Sitter state, it is generic evolution.
On the other hand, if for example a starting or a final w = −1 point is a saddle, then there exists only one trajectory-a separatrix-which is able to reach this saddle (in the past or in the future, respectively), so the evolution is non-generic. Sets of parameters for which the evolution from the de Sitter state to the de Sitter state occurs have been shown in Table 3.   We will analyse generic evolutionary scenarios found in this work since they are more interesting than non-generic ones. Let us start with a discussion of bifurcations for these scenarios. Figures 7-9 present bifurcation diagrams of the local stability of equilibria on the v = 0 nullcline and u = 0 line under variation of the parameter n for fixed ξ = 3 16 and ε = ±1. This corresponds to generic evolutionary scenarios (a) and (b) in Table 3.
In Figure 7 we see the occurrence of two transcritical bifurcations on the v = 0 nullcline.
The first takes place for n = −6 (belongs to the general bifurcation condition n = −4 − 6ξ(6ξ − 1)/ξ) and concerns points B and C. For n < −6 the point B is a saddle and the point C is an unstable node, while after the 'collision' they exchange their stability properties. The second transcritical bifurcation happens for n = −2 (belongs to the general condition n = −4 + 6ξ(6ξ − 1)/ξ) and concerns points A and C. The situation looks similarly as in the previous case: for n < −2 the point A is a saddle and the point C is an unstable node, and these properties are exchanged between points after the bifurcation.
The local bifurcation diagram on the u = 0 line for ε = +1 is shown in Figure 8. The only transition visible in this diagram refers to the point E which undergoes a change from a stable focus via a stable degenerate node (the purple dot) for n = −4 4 7 (in general n = (6−64ξ)/(7ξ)) into a stable node. In fact, this transition is not a bifurcation, because-despite the fact that flows near nodes and foci are neither orbitally nor smoothly equivalent-they are topologically equivalent, which has been shown in Example 2.1 of [16].  the line F , while after the bifurcation the flow, as well below as above the line F , reverses.
Furthermore, combining Figures 7, 8 and 9, we can notice that a bifurcation takes place also for n = −4 in the point (u, v) = (0, 0). For ε = +1, as n approaches −4 from below, the saddle C on the negative part of the u = 0 nullcline and the stable node E on the positive part of the v = 0 axis are approaching the (0, 0) point. They both reach that point for n = −4, and then, when n exceeds that value, the point C becomes a stable node (intercepts the stability from the point E) and is moving towards increasing values of u, while the point E disappears (its coordinates become complex). It is so to speak a combination of the saddle-node and the transcritical bifurcation: the point E behaves as in the saddle-node bifurcation, while the point C behaves as in the transcritical bifurcation. In turn, for ε = −1 the bifurcation is reversed. For n < −4 there is only the saddle C on the negative part of v = 0 axis, which reaches the point (0, 0) for n = −4, and when n exceeds −4, it appears the saddle E, moving towards increasing values of v, while the point C becomes a stable node (gives its instability to the point E) moving towards increasing values of u. In that case, as well as in the previous one, the point E behaves as in the saddle-node bifurcation and the point C behaves as in the transcritical bifurcation. In general, these transitions take place for all values of ξ different than 0 and 1 6 . Let us remark that within the range of parameters of both the generic evolutionary scenarios (a) and (b) there does not occur any bifurcation, therefore it is sufficient to prepare only one phase portrait for each scenario, and this portrait is representative for the whole set of parameters of a given scenario. Also for the non-generic scenarios (e) and (f), we see no bifurcation.
Bifurcation diagrams on both u and v axes under variation of the parameter ξ are presented in Figures 10-11. The choice of the values of fixed parameters, i.e., n = −2 and ε = −1, corresponds to generic evolutionary scenario (c).
In Figure 10 we see the occurrence of two pitchfork bifurcations on the v = 0 nullcline under variation of ξ. The first, taking place for ξ = 0, is subcritical-stable point C is surrounded by unstable points A and B for ξ < 0, and at the critical value of ξ they coalesce into one unstable (saddle) point C (this bifurcation emerges in general at ξ = 0 and n > −4). The second pitchfork bifurcation at ξ = 1 6 appears when all points are unstable-two unstable nodes A and B, surrounding the saddle C emerge, as ξ exceeds this critical value, from the single saddle C (in general, this bifurcation takes place for ξ = 1 6 and n > −4). Finally, for ξ = 3 16 (belongs to the general bifurcation case n = −4 + 6ξ(6ξ − 1)/ξ), it occurs the transcritical bifurcation between points A and C from which, before the bifurcation, the first Table 4: Equilibrium points of the system (35) at infinity together with existence conditions and values of the equation of state parameter.
**The sign of the infinite limit of w φ for G and H depends on values of ξ, n and ε.
is an unstable node and the latter is a saddle. Bifurcation at this point can be seen also in Figure 7, however there the parameter n which was varied.
On the line u = 0, Figure 11 shows other bifurcations under variation of ξ for fixed n = −2 and ε = −1. The saddle D and the neutral line F , which exist for ξ < 0, coalesce in ξ = 0. Above this value, D and F disappear and it appears the stable point E; the line F emerges again at ξ = 1 6 . For ξ = 1 3 the stable point E and the neutral line F coalesce and, as ξ exceeds the bifurcation value, the point E becomes a saddle, while the line F remains neutral, but the flow below and above it reverses.
Combining Figures 10 and 11 we see, that for ξ = 0, n = −2 and ε = −1 all equilibrium points A-F coalesce in (u, v) = (0, 0). It happens for any value of n and ε = ±1 at ξ = 0 and refers to the case of the minimally coupled scalar field. Again, within the range of parameters of the generic scenario (c), we see no bifurcations, which means that one phase portrait will be representative for that case. In turn, the non-generic scenario (g) would require two portraits to be fully represented-the first for 0 ≤ ξ < 1 6 and the second for 1 6 ≤ ξ < 3 16 , while in case of the scenario (h) one portrait would be enough.

Phase portraits of the system
Now, we are ready for a discussion of phase portraits of the system for generic de Sitterde Sitter evolutionary scenarios. Figures 12-14 present phase portraits of system (35) for the evolutionary scenarios (a)-(c). Each of them is representative for its case, as we noticed during the discussion of bifurcation diagrams of equilibria. Portraits have been drawn on the plane (U, V ), which is an orthographic projection of the northern hemisphere W > 0 and the equator W = 0 of the unit Poincaré sphere on which the (u, v) infinite plane had been projected. According to the definition of the Poincaré sphere, coordinates U , V are expressed by u, v as follows [20] and the U 2 + V 2 = 1 circle is corresponding to the infinity of the (u, v) plane, i.e., points for which √ u 2 + v 2 → +∞. From an analysis of the flow of the system (35) on the (U, V ) plane we obtain the information of the equilibria of the system at infinity; it has been summarised in Table 4.
On the V = 0 line we have, according to equations (28) and (29), the kinetic energy parameter εΩ φ,kin → +∞. Moreover, at the point G, the density of energy related to the gravity-scalar field coupling has the limit 6εξ(1 + 2u)v −2 → −∞ · sign(εξ), while at the point H this limit is +∞ · sign(εξ). On the V = √ 1 − U 2 semicircle, if V = 0, the kinetic energy has a finite value Ω φ,kin = εU 2 V −2 , and the coupling part of the energy is equal to zero, thus if Ω φ,kin > 1, the point at infinity is non-physical, according to the condition (36). Hence, at points I, J we have Ω φ,kin = ε and Ω φ,pot = 1 − ε.
The phase portrait for the generic scenario (a), taking place for the canonical scalar field (ε = +1) with ξ = 3 16 and n > 0, is presented in Figure 12. Within the physical region (36), the point E, representing the de Sitter universe, is a global attractor to which trajectories from points A (w φ = −1) and B (w φ = 0) are approaching. De Sitter-de Sitter trajectories (from A to E-marked in purple) can be first either attracted by the saddle C (w φ < − 9 5 ) or both two saddles I and J in turn (w φ = 1 4 both) before reaching the final point. As we can see, in this second case, the value w φ = − 1 3 can be exceeded during evolution, which means the sign ofä can be changed to negative for a certain time.
The flow for the generic scenario (b), representing the phantom scalar field (ε = −1) for parameters values ξ = 3 16 and −2 < n < 0, has been shown in Figure 13, while in the left-hand side of Figure 14 enlargement of an area of the initial stage of an evolution has been outlined. In this case, an evolution from the de Sitter to the de Sitter state also starts at the point A and finishes at the point E. These orbits are transiting through the line F at the point K; when a trajectory is approaching this point from below (in the sense of values of the V coordinate), the equation of state parameter w φ → −∞, and, after exceeding, w φ decreases from the limit +∞-at this point then, the system undergoes the so-called w-singularity. After passing the point K, orbits are heading towards the stable focus E.
Finally, the flow for the generic scenario (c), for which ε = −1, 3 16 < ξ < 1 4 , and n = −2, is topologically equivalent to the flow for scenario (b), as we can see from bifurcation diagramscomparing Figure 9 with 11 and Figure 7 with 10 we notice that in the range of parameters for scenarios (b) and (c), mutual positions of equilibria are the same, except the pair of points A and C, which undergoes a swap of positions and stability properties while passing between those two scenarios; we cannot, however, say it is a transcritical bifurcation of codimension   2 between these points, despite we have to vary values of both ξ and n parameters, because it is enough to change the value of one of these parameters to obtain the same bifurcation. Considering that the non-singular evolution starts from the point A for scenario (b) and from the point C for scenario (c), only the differences between portraits of the (b) and (c) cases are letters marking initial points of this kind of evolution and the shape of the non-physical region described by the condition 36 (it always has boundary at the point A). The last indicates that trajectories from the initial de Sitter point to the point K being attracted on their way by the saddle, which lies on the V = 0 nullcline, are also physically acceptable in the case (c), in contrast to the case (b). The comparison of differences of phase portraits for the scenarios (b) and (c) has been presented in Figure 14.

Evolution of physical quantities
In this section, we will investigate the evolution of physical quantities, such as H,Ḣ + H 2 , Ω φ,kin , Ω φ,pot , over the cosmological time t (the proper time of fundamental observers).
Let us notice firstly that the model we use allows us to determine explicitly only H 2 , φ 2 and the quotientφ/H, which means that we are not able to calculate signs of H anḋ φ, because the model does not distinguish these signs 3 . Knowing that currently H > 0, we can determine regions where H goes to negative values by searching for H 2 = 0 points and checking a value ofḢ + H 2 there.
In order to present the time evolution of the scale factor a and its derivatives, we need to transform equations (27) and (30); thus we obtain respectively where the overdot symbol (˙) denotes a derivative with respect to t. The reparameterisation of time for dynamical system (35) gives the new time τ which is related with t by To reconstruct the time t we need to determine the sign of H and integrate the foregoing equation. In this work we assume arbitrarily τ = 0 such that w φ (τ ≥ 0) ≈ −1 at least to the third decimal place. Moreover, we put t(τ = 0) = 0. Having the time t, it is possible to calculate the scale factor directly from the definition of the Hubble function For every generic evolutionary scenario (a)-(c), there is no point on de Sitter-de Sitter orbits, for which H 2 = 0, with the exception of the initial point A for (a) and (b) scenarios, where the limit of H 2 can be equal to zero 4 . This means, that all scenarios represent the expanding universe (H > 0) at any finite time. Since H = 0 corresponds to the static universe, scenarios (a) for any n and (b) for −1 < n < 0 represent the initially static universe, which means that the scale factor a > 0 at t → −∞. In case of scenarios (b) for −2 < n ≤ −1 and (c) for any ξ it is a → 0 at t → −∞.  The evolution of physical quantities over the cosmological time t for the evolutionary scenario (a) has been presented in Figure 15. The left column corresponds to an orbit passing near the equilibrium point I, whereas the right column corresponds to an orbit passing near the equilibrium point C. In plots we see the occurrence of initial static state (w φ = −1, H = 0 and H = 0-from the side of t → −∞), and final de Sitter state (w φ = −1,Ḣ = 0 and H = const > 0-from the side of t → +∞), between which a non-exponential evolution takes place. For trajectories attracted by the point I it is possible for the universe to undergo a decelerated expansion (w φ > −1/3,Ḣ + H 2 =ä/a < 0) for a particular interval of time. The kinetic energy density Ω φ,kin is changing from +∞ to 0 during the evolution. Figure 16 shows the evolution of physical quantities over the cosmological time t for the evolutionary scenarios (b) (the left column) and (c) (the right column). Two de Sitter states emerge from both sides of t → ±∞ (w φ = −1,Ḣ = 0, H = const > 0). During evolutions (b) and (c), it appears the discontinuity of w φ , which takes place when the orbit intersects the equilibrium line F . As the solution is approaching this line, the equation of state parameter  w φ → −∞. After the solution exceeds line F , w φ is decreasing from the limit +∞; it appears the decelerated expansion (H 2 +Ḣ < 0) for certain time. The kinetic energy Ω φ,kin is changing from −∞ to 0 throughout the evolution.
In Figures 15-16, we can see the presence of the initial exponential expansion, stretching to t → −∞, which models the inflation in the early universe. The second exponential expansion, stretching to t → +∞, describes the late time accelerated expansion, on the condition we add the matter content of the current universe. The full model, taking explicitly the matter component into account, can be an interesting aim of further research.

Conclusions
In the paper, we have used bifurcation theory methods to study the sensitivity of the model under a variation of the parameters. To show the power of these methods we have considered cosmological models which provide the explanation of the acceleration conundrum and are considered in the context of the inflation.
From our investigation, the following conclusions can be derived: 1. The application of bifurcation theory allow us to distinguish some class of emergent cosmological models from Einstein static universe and de Sitter universe. The first of these cases takes place for the canonical scalar field with the value of the non-minimal coupling parameter ξ = 3 16 and the exponent appearing in the potential n > 0, while two other cases concern the phantom scalar field for ξ = 3 16 and −2 < n < 0 or for 3 16 < ξ < 1 4 and n = −2. Furthermore, we have indicated five cases of non-generic de Sitter-de Sitter evolutionary scenarios for either canonical or phantom scalar field.
2. In the phase space, there are present evolutional paths without the initial singularity, allowing to explain both the inflation and the accelerated phase of the expansion at the late times.
3. There are two types of initial states from which the Universe starts the evolution. In the first scenario, the Universe is emergent from the de Sitter state and in the second one, the Universe is emergent from the static universe.
4. From the bifurcation analysis, we have obtained a pair of the critical values of the parameters (ξ, n) which corresponds to bifurcation values.

5.
We have demonstrated that the methods of the bifurcation theory allows us to distinguish the regions of parameter values with the generic and non-generic evolution paths.