Cosmological dynamics and bifurcation analysis of the general non-minimally coupled scalar field models

Non-minimally coupled scalar field models are well-known for providing interesting cosmological features. These include a late time dark energy behavior, a phantom dark energy evolution without singularity, an early time inflationary universe, scaling solutions, convergence to the standard $\Lambda$CDM, etc. While the usual stability analysis helps us determine the evolution of a model geometrically, bifurcation theory allows us to precisely locate the parameters' values describing the global dynamics without a fine-tuning of initial conditions. Using the center manifold theory and bifurcation analysis, we show that the general model undergoes a transcritical bifurcation, which predicts us to tune our models to have certain desired dynamics. We obtained a class of models and a range of parameters capable of describing a cosmic evolution from an early radiation era towards a late time dark energy era over a wide range of initial conditions. There is also a possible scenario of crossing the phantom divide line. We also find a class of models where the late time attractor mechanism is indistinguishable from that of a structurally stable general relativity based model; thus, we can elude the big rip singularity generically. Therefore, bifurcation theory allows us to select models that are viable with cosmological observations.

Non-minimal coupled scalar field models are well-known for providing interesting cosmological features. These include a late-time dark energy behavior, a phantom dark energy evolution without singularity, an earlytime inflationary Universe, scaling solutions, convergence to the standard ΛCDM, etc. While the usual stability analysis helps us determine the evolution of a model geometrically, bifurcation theory allows us to precisely locate the parameters' values describing the global dynamics without a fine-tuning of initial conditions. Using the center manifold theory and bifurcation analysis, we show that the general model undergoes a transcritical bifurcation, predicting us to tune our models to have certain desired dynamics. We obtained a class of models and a range of parameters capable of describing a cosmic evolution from an early radiation era towards a late time dark energy era over a wide range of initial conditions. There is also a possible scenario of crossing the phantom divide line. We also find a class of models where the late time attractor mechanism is indistinguishable from a structurally stable general relativity-based model; thus, we can elude the big rip singularity generically. Therefore, bifurcation theory allows us to select models that are viable with cosmological observations.

I. INTRODUCTION
Non-minimal coupled scalar field models are often used to explain various cosmological observations. These models naturally arise from the quantum corrections to the scalar field theory and motivated by high energy physics such as superstrings and grand unified theories [1]. Further, these models provide a natural solution to the problem associated with the energy scale difference between inflation and the Universe's dark energy (DE) era [2].
Most of the cosmological model's governing equations are nonlinear and pose a severe impediment to extract exact analytical solutions. However, one can infer the global asymptotic behavior described by the cosmological equations using the advanced tools of dynamical systems. The main advantage is that we can represent the Universe's history geometrically. One can also predict the sensitivity of the solution to initial conditions. Dynamical system methods have been used extensively in cosmology; see [3][4][5][6][7][8][9][10][11][12][13][14][15] for relevant work and [16] for a comprehensive review.
The dynamical system of most cosmological models usually contains parameters. One can determine the system's global dynamics for fixed values of parameters using the formal stability analysis. On the other hand, to understand how the global dynamics changes with a change of parameters, the bifurcation theory plays a crucial role (see Refs. [17][18][19] for detailed information). The dynamical system's nonlinear nature usually leads to vital structures of the solutions, such as bifurcations and chaos. A more in-depth analysis of such forms is interesting from an observational perspective (e.g., see [20]).
One of the bifurcation theory's novelties is that we can use it to classify the Universe's evolution into two categories: generic and non-generic evolution [21]. While the former occurs for various solutions over a wide range of initial conditions, the latter corresponds to a particular solution for a given initial condition. The parametric relation associated with a * sjwomkhyllep@gmail.com † jibitesh@nehu.ac.in non-generic scenario forms a bifurcation boundary between regions of different generic cases in the parameter space. Nongeneric evolution is also exciting but requires fine-tuning of initial conditions. In some cases, generic evolution emerges from non-generic one in the form of bifurcation. Therefore, bifurcation theory can help extract a class of models describing the observed dynamical evolution irrespective of initial conditions for a wide range of parameters. To have a clear picture of how the bifurcation phenomenon depends on model parameters, one has to use bifurcation diagrams. These diagrams stratify the parameter space into different regions, each with distinct dynamical behavior. The steps involved in the bifurcation analysis are two-fold. Initially, we extract the range of parameters for both the generic and the non-generic evolution. Then, we analyze different qualitative behaviors that arise from each scenario. Another novelty of bifurcation theory is that it allows us to identify a structural stable model's emergence from a structurally unstable one. For instance, Szydlowski and Tambor showed that the notion of bifurcation and structural instability could be instrumental in detecting the emergence of the structurally stable ΛCDM model from the structurally unstable CDM model [22]. Kokarev further extended a similar analysis to various Friedmann-Robertson-Walker (FRW)-models [23]. Usually, structurally stable models are physically viable and hence fit with most observations. In a 2-dimensional system, Peixoto's theorem completely characterizes the structurally stable vector fields, which guarantee their generic behavior. Identifying structurally stable models is useful when the prediction of model parameters from the empirical analysis is unsettled. Thus, bifurcation theory helps in finding observationally viable models and further endows the usual stability analysis.
In most of the dynamical analysis for the non-minimal coupled scalar field models, the dynamical variables are constructed for a specific case of coupling or potential functions in flat or curved spacetime [12,[24][25][26][27][28][29][30][31]. However, the analysis for a general non-minimal coupled scalar field model will certainly help us to identify classes of viable models. The extension to a broad class of scalar field potentials and couplings might help us to connect the phenomenological models with some high-energy physical theories. Therefore, it will be sci-entific and economical to carry out the dynamical analysis for a broad class of coupling functions and potentials.
We find in the literature that bifurcation phenomena arise naturally in cosmological models. For instance, Ref. [32] shows that in FRW-models with perfect fluids and the cosmological constant, the expanding and contracting deSitter Universe arise as bifurcation. It is worth mentioning that interesting bifurcation scenarios were reported in the Randall-Sundrum braneworld model [33], interacting Veneziano ghost DE [34], Brans-Dicke model [35], non-minimal coupled scalar field model [31,36] etc. Recently, bifurcation scenarios and chaos were discussed in the context of Hořava-Lifshitz gravity [37], non-minimal coupled scalar field with Ratra-Peebles potential [21], interacting f (T ) gravity [38] and bulk viscous cosmology [39]. These recent work show that the study of bifurcation is important in cosmology, giving rise to interesting scenarios.
In the non-minimal coupled scalar field context, interesting bifurcation scenarios were reported for a specific coupling and potential function. For instance, Hrycyna et al. [31] obtained a particular bifurcation value of a coupling constant for the case of a constant potential. Then, Szydlowski et al. in [36] analyzed the phase space's structural stability of a specific coupling model for a broad class of potentials. They found that an exponential potential constitutes a structurally stable model. Using the bifurcation methods, Humieja et al. in [21] extract the conditions of model parameters under which a specific non-minimal coupling with Ratra-Peebles potential generically evolve from an early de Sitter to a late time de Sitter state. The analysis in [21] was performed in the absence of a matter component. We extend the analysis for general coupling and potential functions along with the matter component in the present work. To meet our objective, we consider a different choice of dynamical variables to encompass a broad class of models. By employing bifurcation methods, we obtain a class of models and pinpoint the range of parameters capable of describing a cosmic evolution over a wide range of initial conditions from an early radiation era towards a late time DE era. We also found that the system undergoes a transcritical type of bifurcation, which predicts how to tune our models to have certain desired dynamics.
The bifurcation theory's concrete tools have been used extensively in various fields. However, they have not been applied systematically in many cosmological systems, particularly for the non-minimal coupled scalar field. Thus, it is imperative to use bifurcation theory to identify a class of scalar field models describing some of the main generic cosmic evolution. Therefore, the present work serves as an introductory analysis for scalar field models required to test against interesting observational signatures.
The paper's order is as follows: In Sec. II, we briefly discuss the framework of a non-minimal coupled scalar field model. We follow this by a dynamical system analysis of a nonminimal coupled scalar field model for a broad class of coupling function and potential in Sec. III. In Sec. IV, we demonstrate the dynamics by an example using the quadratic coupling functions and the power-law form of potentials. Within this section, we perform the stability analysis of critical points in Sec. IV A and the discussion on bifurcation scenarios in Sec. IV B. Lastly, we summarized the work in Sec. V.

II. NON-MINIMAL COUPLED SCALAR FIELD MODEL
We consider a model of a non-minimal coupled scalar field and a barotropic fluid in the present work. Here the nonminimal coupled scalar field is playing the role of DE, while a barotropic fluid is the matter component of the Universe. The total action is given by [40][41][42] S = 1 where the integration is taken over a 4-dimensional Lorentzian curved spacetime manifold. In the above action, κ 2 is the gravitational constant, R is the Ricci scalar, g is the determinant of the spacetime metric g ab (a, b = 0, 1, 2, 3), F is the coupling function, V is the potential of a scalar field φ and L m is the matter Lagrangian. We have used the units where c = 1 and is the positive parameter having the dimension of length. For a fixed scalar field, the above action reduces to the case of general relativity (GR) with a potential playing the role of a cosmological constant. While the case F(φ ) = 1 of (1) corresponds to the minimal coupled scalar field, F(φ ) = φ 2 ω reduces to the Brans-Dicke gravity limit, with ω as the Brans-Dicke parameter [27,[43][44][45]. On varying the action (1) with respect to the metric g ab , one can obtain the modified Einstein's field equation as where G ab is the Einstein tensor, ≡ ∇ a ∇ a with ∇ a as the covariant derivative with respect to the metric and T ab is the matter energy-momentum tensor given by In the above equation, ρ m and p m are respectively the energy density and pressure of the barotropic fluid, and u a is a four-velocity vector of the fluid. One interesting feature of the action (1) is that the effective Newton's gravitational parameter depends on the coupling function F, i.e., on the scalar field as The negative values of G eff and hence of F indicates the ghost instability in the theory [46]. On varying the action with respect to the scalar field φ , we get where the notation (·) ,φ denotes a derivative with respect to φ . Note that the second term of (5) arises from the non-minimal coupling of scalar field to gravity. At a very large scale, consistent with the observed data, we assume the homogeneous and isotropic Universe whose evolution is determined by the scale factor a(t) associated with the FRW metric where t is the coordinate time and x, y, z are the Cartesian coordinates. Under this metric, the field equations (2) and the Klein Gordon equation (5) can be reduced to the ordinary differential equations 2FḢ In the above equations, w is an equation of state (EoS) defined by a relation p m = wρ m and the upper dot denotes derivative with respect to t. The scalar field describes DE and for simplicity, we shall consider for matter a single barotropic fluid with a constant w, constrained to be between 0 and 1. While a non-relativistic dust fluid corresponds to w = 0, relativistic radiation fluid corresponds to w = 1 3 . We note here that one needs to consider a two-fluid model containing radiation and dust fluids for a more phenomenologically interesting case. Further, on assuming the conservation of the matter energymomentum tensor i.e., ∇ a T ab = 0, under a metric (6), one can obtain the conservation equatioṅ To determine the energy density contribution of each component, we introduce the relative energy densities of the scalar field and that of the barotropic fluid, respectively as These energy densities are connected by the Friedmann constraint (7) as While we can identify the first two terms of (11) as the kinetic component of the relative energy density of the scalar field, the last term corresponds to the relative potential energy density component of the scalar field. From (13), we can define the matter domination as a scenario where Ω m ≈ 1 and Ω φ ≈ 0. Similarly, we can also define from (13), a kinetic dominated solution or potential dominated solution when the first two terms or the last term in (11) dominate over the others, respectively. As expected, the above quantities (11) and (12) reduce to the familiar relative energy densities of a scalar field and matter for the minimal coupled case (i.e., F(φ ) = 1) respectively. Further, the quantity (11) reduces to the relative energy density of the cosmological constant for a non-dynamical scalar field (i.e., the GR case). In the context of minimal coupling, the above relative energy densities are usually bounded within the interval [0, 1]. However, this is not necessarily true in non-minimal coupling due to coupling function F. Under the physical assumption ρ m ≥ 0 and an attractive gravitational force, F(φ ) > 0, the problem of negative Ω m does not arise. Nonetheless, due to the second term in the right-hand side of the equation (11) for Ω φ , there is also a possibility for Ω φ to be negative. Thus taking into account the above conditions, the relation (13) implies that in the matter domination, the dust fluid can be relatively overdense (i.e., Ω m > 1) compared to the corresponding ΛCDM case.
The cosmological equations (7) and (8) can be rewritten in the form where ρ eff is the effective energy density and p eff is the effective pressure of all the components which are respectively given by Using equation (7), the effective EoS of all the components w eff defined as p eff ρ eff is given by While for the accelerated behavior of the Universe, one requires the condition w eff < − 1 3 , super-accelerated Universe or phantom dominated Universe demands w eff < −1. It is worth noticing from (18) that in the GR limit, within the matter domination epoch, we have w eff = w and under the scalar field potential dominated epoch (i.e., cosmological constant epoch) w eff = −1.
The above equations (7)-(10) are complicated to solve analytically, yet, by recasting them into a dynamical system, one can still obtain important information on the characteristics of solutions. Therefore, in the next section, we shall analyze the dynamics of a general class of non-minimal coupling scalar fields using dynamical system techniques.

III. DYNAMICAL SYSTEM ANALYSIS
In order to qualitatively analyze the background cosmological dynamics of the present model, we shall convert the cosmological equations (7)-(10) into a dynamical system using the following set of normalized variables [16]: We note here that the chosen variables are well-defined for F > 0, i.e., attractive gravity, which is also free from any ghost instability, even though the case F < 0 may lead to physically interesting scenarios [47]. From the cosmological equations (7)-(10), we see that there are basically four variables H, φ ,φ , ρ m . As we have considered the usual Hnormalized variables, the variable H is being absorbed by other variables, so we are left with three variables [16]. Since the H-normalized variables are connected by the Friedmann constraint (7), the number of independent variables reduces to two. The extra variables λ F , λ V are introduced to monitor the overall effect of coupling function and potential on the dynamics. It is important to note here that the above choice of variables fails for static Universe H = 0. However, these variables are of physical interest as the energy density of each component can be easily tracked in terms of these variables. In this work, we shall focus on the case of an expanding Universe i.e., H > 0 as favored by various present observational data. Therefore, we can choose the above normalized variables without any extra concern. Employing the variables (19), the cosmological equations (7)- (10) can be re-written as the following dynamical system: where . The prime notation denotes the differentiation with respect to the number of e-folds N = ln a(t). We note that the above system (20)-(23) reduces to the minimal coupling case for λ F = 0 [48].
For the above dynamical system to represents an autonomous system of equations, we consider a class of coupling function F and potential V where Γ F , Γ V can be written as functions of λ F , λ V respectively [49]. If λ F = λ F (φ ) is invertible, then we can express φ as function of λ F . As Γ F is a function of φ , therefore, we can also express Γ F as a function of λ F . Similarly, one can express Γ V as a function of λ V . In general, the quantities Γ F , Γ V may not be a functions of variables λ F , λ V . In such a case, one has to consider the higher derivatives of the scalar field function [50] or consider new dynamical variables [51]. We note that the above system has an invariant submanifold y = 0, as y vanishes when y = 0. This submanifold corresponds to a scenario where the scalar field potential vanishes. Physically, it means that if there is no potential source, the scalar field will not evolve. Further, depending on F and V (hence in the form of Γ F and Γ V ), the system also contains λ F = 0, λ V = 0 as invariant submanifolds. Therefore, a global attractor (if exists) should lie at an intersection of all these invariant submanifolds [16]. On the other hand, the absence of a global attractor makes the application of bifurcation theory more appealing as the evolution depends on the values of parameters and initial conditions.
Using the dynamical variables (19), one can express various cosmological parameters viz., the relative energy density parameter of the scalar field (Ω φ ) and of matter (Ω m ), the EoS of the scalar field (w φ ) and the effective EoS (w eff ) respectively as Notably the coupling term Hφ F ,φ of (7) can steer the value of w φ to ±∞ during the matter domination epoch. However, this does not cause any physical singularity problem as the effective EoS w eff remains smooth and finite. The divergence behavior of w φ reduces as the value of λ F approaches zero i.e., as the model approaches the minimal coupling case. By imposing the physical constraint ρ m ≥ 0 on the relation (13), the dynamical variables (19) obey the constraint Hence, the phase space of the system is given by From the cosmological equations (7)-(10), one can solve the scale factor a(t) evaluated at the critical point by re-writing the equations in terms of the dynamical variables as where Integrating equation (30), we get where a i and t i are constants of integration. We recall that 0 < β < 1 corresponds to a decelerated expanding Universe, while β > 1 corresponds to an accelerated expanding Universe. In order to extract the dynamics of the above system, we will carry out the standard procedures of the dynamical system analysis [16]. The critical points (A 1 , A 2± , A 3 , A 4 ) of the system (20)-(23) for the general case of F and V are presented in Table I along with the corresponding values of the cosmological parameters. The corresponding eigenvalues of the perturbed matrix of each critical point are presented in Table II. The existence and stability of each critical point can be determined without specifying the potential and coupling function by treating λ F * and λ V * as parameters. Therefore, there are as many critical points of the system (20)-(23) as the number of parameters λ F * and λ V * . Note here that λ F * and λ V * TABLE I. Critical points of the system (20)- (23).
TABLE II. Eigenvalues of the critical points of the system (20)-(23) presented in Table I. Here: x k denotes the corresponding x-component of a critical point.
The quantities Γ F and Γ V denote the derivatives of Γ F and Γ V with respect to λ F and λ V respectively. In Tables I and II, we have The behavior of the system (20)-(23) might change dramatically due to a small change of the parameters emerges from potential and coupling. Consequently, it will lead to a change in the phase space's topological structure and gives rise to bifurcation. To have a general information on the effect of various parameters on the dynamics of a system (20)-(23), we present the bifurcation diagrams for each critical point in Fig.  1. These diagrams comprise of finite number of regions in the parameter space (λ F * , λ V * ). Each region in a bifurcation diagram corresponds to parameter values with distinct dynamical behavior [18]. Bifurcation curves separating different regions in a diagram corresponding to the specific relation between parameters in which the perturbed matrix evaluated at a critical point has at least one zero real part eigenvalue. When the perturbed matrix evaluated at a critical point has at least one zero real part eigenvalue, a critical point is said to be non-hyperbolic. Otherwise, it is hyperbolic. For a hyperbolic point, we can use linear stability analysis to determine the stability of a point. However, for non-hyperbolic point, one has to analyze beyond the linear stability analysis via sophisticated tools of dynamical systems such as the center manifold theory (see Ref. [19] for details). Thus along the bifurcation curves, one has to investigate the nature of points by the center manifold theory. However, we shall postpone such analysis to a concrete model in Sec. IV and refrain from the general case analysis as the equations involved are complicated and not very illuminating. In what follows, we summarize each critical point's nature and identify the possible bifurcation scenarios: • Point A 1 corresponds to a solution with a vanishing potential component. For this point, the exponent for the solution (31) is given by It can be checked that 0 < β A 1 < 1 for 0 ≤ w ≤ 1 and any choice of λ F * . Hence, this point corresponds to a decelerated expanding solution for any choice of model parameters even though the scalar field can possibly behave as the quintessence field (−1 < w φ < − 1 3 ). This point can be either stable or saddle depending on the values of λ F * , λ V * and w. For instance, as shown in Fig.  1(a) for w = 0 case, this point behaves as stable node in region I only if Γ F (λ F * ) > 0 and G(λ F * , λ V * ) > 0. Region III represents the region of stable node only if Γ F (λ F * ) < 0 and G(λ F * , λ V * ) < 0; otherwise all regions in the (λ F * , λ V * ) parameter space represent the regions of saddle for any choice of coupling function F and potential V . However, for the case of radiation (w = 1 3 ), it is always saddle. Further, for λ F * = 0, this point corresponds to a usual decelerated matter dominated solution of the minimal coupled case (Ω m = 1, w eff = w). • Points A 2± correspond to kinetic dominated solutions which exist for any choice of model parameters. In this case, we have which are both positive and less than one and hence, correspond to decelerated expansion. The bifurcation diagrams for these points are given in Figs. 1(b), 1(c). These points can be either unstable node or saddle depending on the form of potential and coupling function. For example, region I in Fig. 1(b) represents region of unstable node of point Fig. 1(c) represents the region of unstable node of Otherwise, they are saddle in nature for any choice of parameters. When (λ F * , λ V * ) = (0, √ 6) and (0, − √ 6), points A 2+ and A 2− correspond to stiff matter dominated solutions (Ω φ = 1, w eff = 1) of the minimal coupled scalar field.

• Point A 3 corresponds to a scaling solution and exists for
all values of the model parameters except when λ V * = 0. This point can either be a stable node or stable focus or behaving as a saddle. The bifurcation diagrams for this point for w = 0 case is given in Fig. 1(d). In this plot, when Γ F (λ F * ) > 0 and G(λ F * , λ V * ) > 0, region II represents region of stable focus, regions VI and X represent regions of stable node and the remaining regions represent the regions where this point behaves as saddle. However, when Γ F (λ F * ) < 0 and G(λ F * , λ V * ) < 0, region I represents region of stable focus, regions V and IX represent regions of stable node and the remaining regions represent the regions where this point behaves as saddle. It is important to note that as the parameters values change across the bifurcation curves C5 and C6 of Fig. 1(d), the property of the Universe changes from a decelerated scaling solution (w eff > − 1 3 , Ω m > 0) to a decelerated scalar field dominated solution (w eff > − 1 3 , Ω m = 0). For this point, the exponent β in (31) is given by Therefore, it represents a decelerated expansion when 0 < λ V * λ V * −λ F * < 3 2 (w + 1) and an accelerated expanding Universe when . When λ F * = λ V * , it corresponds to an effective cosmological constant behavior (w eff = −1), but it is unphysical as Ω m < 0. Further, from the bifurcation diagram of this point, we have checked that within the unstable (or saddle) accelerated regions of the parameter space, this point is unphysical.
• Point A 4 corresponds to a scalar field dominated solution (Ω φ = 1). This point disappears for values of It can either be a stable node or unstable node or saddle depending on the choice of model parameters (see Fig. 1(e)). For Γ F (λ F * ) > 0 and G(λ F * , λ V * ) > 0, while region I represents a region of unstable node, regions III and VIII represent regions of stable node. The remaining regions represent the regions where this point behaves as a saddle. For Γ F (λ F * ) < 0 and G(λ F * , λ V * ) < 0, region II represents a region of unstable node, regions IV and VII represent regions of stable node and the remaining regions represent the regions where this point behaves as a saddle. The corresponding exponent for the scale factor solution (31) is given by  Fig. 1(e)) to a stable node (some parts of regions VII and VIII of Fig. 1(e)). It is worth mentioning that this point corresponds to an effective cosmological constant behavior (w eff = −1) when λ F * = λ V * or λ F * = λ V * 2 . The solution corresponds to λ F * = λ V * 2 is of interest as it is identical to the de Sitter solution in GR (since φ is constant in time). From the bifurcation diagram (i.e., Fig. 1(e)), one could confirm that there is a possibility that this point is stable when λ F * = λ V * 2 . This point's stable nature explains the possible late time convergence behavior of the scalar-tensor theory towards GR. In other words, such a model converges towards a structurally stable GRbased model. For coupling function F(φ ) = 1 + ξ φ 2 and V = V 0 φ n , this condition is met for n = 4. Also when λ F * < λ V * < 2λ F * , this point behaves as a stable phantom-like attractor. A similar result for a particular potential case has been reported in [52]. Therefore, our present work provides a framework for choosing proper coupling and potential functions to get interesting dynamics for a broad class of non-minimal coupled scalar field models.
As there is no stable critical point lying on the intersection of invariant sub-manifolds y = 0, λ F = 0, λ V = 0, therefore, the above system does not contain any global attractor. Further, from the above analysis, one can see that the present model exhibits interesting solutions that can describe various cosmological eras of the Universe. For instance, the decelerated scalar field dominated solution can describe the early radiation era, the matter scaling solution that can describe the intermediate dark matter (DM), and the late time scalar field dominated solution describing the DE dominated era. Numerically, we observe that the local stability behavior of critical points changes for the coupling and potential parameters (Fig. 1), and hence the system as a whole undergoes bifurcation. Using the bifurcation diagrams of various critical points, we can classify the parameters for which the evolution corresponds to the generic evolution and non-generic evolution. Geometrically, when an orbit evolves from an unstable point and then settles in a stable point, it is called generic evolution. However, if the initial or final point is a saddle, it is called non-generic evolution. From the above analysis, we found that only class of models where the parameters λ V * , λ F * belong to the common region of region I of Fig. 1(b) and region VII of Fig. 1(e); and also in the common region of region II of Fig. 1(c) and region VIII of Fig. 1(e) can possibly lead to physically interesting generic evolutionary scenarios. These correspond to the evolution from a decelerated kinetic dominated unstable node point A 2+ /A 2− to an accelerated potential dominated stable point A 4 via a matter scaling solutions A 1 or A 3 . However, solutions evolving near an intermediate point A 3 correspond to unphysical solutions (Ω m < 0). Therefore, the sequence of generic cosmic viable evolution is given by A 2± → A 1 → A 4 . Thus, bifurcation diagrams allow us to classify the coupling and potential functions, describing interesting cosmological dynamics.
Further from the bifurcation diagrams of each point, we see that there is an occurrence of transcritical bifurcation between critical points. For example, critical points  1(e)). This type of bifurcation occurs when two critical points interchange their stability properties at the bifurcation curve [18]. We can summarize these bifurcation scenarios as follows:

Critical points A 3 and A 4 interchange a saddle and stable node behavior along the bifurcation curves
for fixed value of w. In particular, for w = 0, the above curves are represented geometrically by curves C5 and C6 of Figs. 1(d) and 1(e).
Note that one can prove the existence of transcritical bifurcation analytically by using the Sotomayor's theorem [17][18][19]. However, since the bifurcation parameters λ V * , λ F * do not appear explicitly on the system (20)-(23), therefore, we postpone the analytical proof to a concrete example of coupling function F and scalar field potential V . Even though we can determine the general properties without specifying the concrete model, to understand the cosmological applications of a general model and better investigate its dynamics, we must assume a concrete example. Therefore, in the next section, we consider a specific model and analyze its cosmological dynamics in detail.

A. Stability Analysis
Here, we shall examine the case where the non-minimal coupling is of the form F(φ ) = ξ φ 2 and the scalar field potential is of form V (φ ) = V 0 φ n (power-law). For this example, we have λ F = −2 ξ , λ V = −n ξ and hence they are constants. This particular example is inspired physically by the string-dilaton, Brans Dicke actions and several effective quantum field theories [53]. Mathematically, this form of coupling and potential agrees with the Noether symmetry approach of the Lagrangian given by (1) and hence, leads to physically interesting exact solutions [54,55]. This model is also compatible with the solar system constraint test [56]. Various cosmological data constrain the value of ξ to be 10 −2 at 95% confidence limit. However, a degeneracy between the value of ξ and the present Hubble constant H 0 allows a larger value of ξ [57,58].
Stability analysis for this particular example has been performed earlier in [27] where the main focus is on the hyperbolic points. However, a discussion on the condition for non-hyperbolicity of points has not been performed. The nonhyperbolic nature of the critical point is important to analyze the bifurcation scenarios. Since for this concrete example, λ F and λ V are fixed, the system (20)-(23) reduces to following two-dimensional system: y = y 2(6 ξ + 1) 6(8 ξ + w + 1) + 2 ξ x (−2(6 ξ − 3 w + 2) The physical phase space of the reduced system is The critical points for the above system are given in Table  III. We note that critical points B 1 , B 2± , B 3 Stable node/Unstable node/Saddle points A 1 , A 2± , A 3 , A 4 respectively for this concrete example of coupling function and potential. The existence and stability behavior of these critical points are similar to the general case (see Table III). The bifurcation diagrams exhibiting the stability regions of critical points are given in Fig. 2. It is worth noting that each critical point shows a non-hyperbolic behavior for different values of parameters (represented by black colored curves in Fig. 2). As analyzed in the general case (see Sec. III), critical points B 1 and B 3 coincide along the curve D1 in which both of them are non-hyperbolic. Therefore, in this case, it is sufficient to analyze the non-hyperbolic nature only for a point B 1 . The analysis of center manifold theory for this case is performed in the appendix A and it is found that point B 1 behaves as a saddle. For a detailed mathematical background on the center manifold theory, we refer to [19]. Further, points B 2+ and B 2− coincide with point B 4 and are non-hyperbolic along the curves D2 and D3 respectively. Also, critical points B 3 and B 4 are non-hyperbolic and coincide along the curves D4 and D5. Therefore, in each case, we shall analyze the non-hyperbolic property for a point B 4 only. The analysis performed in the appendix B reveals that point B 4 behaves as a saddle along each bifurcation curve. As the phase space (39) is in general not compact, for the sake of completeness, we analyze the nature of the system (37)- (38) at infinity by employing the Poincaré's projection method involving the following transformation [19]: The compactified phase space of the resulting system is therefore given by with R = 1 − x 2 r − y 2 r . The critical points near infinity i.e., x 2 + y 2 → ∞ correspond to points on the circle (x r , y r ) ∈ R 2 : x 2 r + y 2 r = 1 .
We should take due care in applying the Poincaré compactification method, resulting in a wrong phase space topology, and hence false properties of the solutions [59]. However, in this case, the above choice of the dynamical variables does not destroy the phase space's global structure since H > 0 [60]. The method's advantage is to capture the possible critical points hidden at infinity that cannot be detected by the finite analysis. On employing this method, we found that there are four critical points near infinity viz., (x r , y r ) = (±1, 0), (0, ±1) lying on the equator of the Poincaré's sphere. The analysis shows that points P ∞ 1 (1, 0), and P ∞ 3 (−1, 0) are saddle in nature. Points P ∞ 2 (0, 1), P ∞ 4 (0, −1) are non-hyperbolic and hence the stability needs to be checked with an extra effort by the center manifold theory. However, by using a quick numerical check, we found that point P ∞ 2 is saddle and P ∞ 4 is stable for various values of model parameters n, ξ and w = 0. We note that of all critical points at infinity, only a saddle point P ∞ 2 corresponds to an accelerated solution, but it is unphysical. Thus, critical points at infinity cannot describe either late time or early time behavior of the Universe and hence are not phenomenologically interesting. Therefore, here we do not present a detailed calculation of the analysis at infinity. In the next section, we discuss the possible bifurcation scenarios occurring at the nonhyperbolic condition of each critical point.

B. Bifurcation Scenarios
In this section, we shall discuss the occurrence of local bifurcation of the system (37)- (38) with respect to parameters ξ and n. Then we extract the condition on ξ and n under which the present model describes the generic evolution of the Universe. In the case of a two-dimensional system, we can completely characterize the structural stability by Peixoto's theorem. However, we cannot extend the theorem to a system of dimensions greater than two [19]. As a result of Peixoto's theorem, a structurally stable system guarantees an open dense subset of initial conditions leading to a generic evolution (see appendix C).
In general, the necessary condition for the occurrence of bifurcation of the system is the non-hyperbolicity of the critical point. However, in the two-dimensional system, according to Peixoto's theorem, the existence of non-hyperbolic critical points also implies the structural instability of the system [19]. As we have seen in the previous section, the system (37)-(38) contains non-hyperbolic points P ∞ 2 and P ∞ 4 on the Poincaré sphere, therefore, the vector field of the system is structurally unstable. Moreover, finite critical points can be non-hyperbolic for some values of ξ , n and w. In what follows, similar to the general case, we again prepare the bifurcation diagrams (Fig. 2) for each critical point in the (ξ , n) parameter space and then apply the Sotomayor's theorem (see the appendix C for the statement). The theorem's main aim is to analytically investigate and specify the types of bifurcation.
The local bifurcation diagram for point B 1 is given in Fig.  2(a). A topological change occurs as this point changes from stable node to saddle along the curve D1 via a non-hyperbolic saddle node.
Another bifurcation occurs for points B 2+ and B 2− along the curves D2 and D3 respectively, where both the points undergo an upheaval from unstable node to saddle via a non-hyperbolic saddle node (see Figs. 2(b), 2(c)).
Point B 3 changes its stability from a stable node to a saddle along the bifurcation curves D1, D4, D5 of Fig. 2(d). For a particular case, n = 2, this point exhibits an effective cosmological constant behavior and for n = 6, the point's dynamics change from an unaccelerated to an accelerated behavior.
Out of the above mentioned critical points, point B 4 is an interesting point which can explain the late time behavior of the Universe. This point undergoes a stability change from an unstable node to a saddle along the curves D2 and D3 of Fig.  2(e). Again, a change from a stable node to a saddle occurs when the point passes through the bifurcation curves D4 and D5 of Fig. 2(e). As this point can represent interesting late time Universe, using bifurcation diagram (Fig. 2(e)), we summarize the stability property of this point for different range of ξ , n, w as follows: Fig. 2(e)), or, Fig. 2(e)).
It is worth mentioning that independent of the value of ξ , this point describes a stable phantom attractor behavior for 2 < n < 4. Also, when n = 4, point B 4 belongs to a stable region of parameter space. Therefore, as discussed in Sec. III (i.e., model with the condition λ F * = λ V * 2 ) this point corresponds to a deSitter solution in GR for n = 4. Similar result also holds for model considered in [21] with biquadratic potential.
From the above discussion on the bifurcation diagrams, we found the possibility of transcritical bifurcation exhibited by the system (37)- (38). In what follows, we summarize the occurrence of bifurcation and mathematically analyze the transcritical bifurcation with ξ as the bifurcation parameter using Sotomayor's theorem. Proof: In this case, the bifurcation value is in which the points B 1 and B 3 coincide and the real part of the eigenvalue corresponding to common critical point vanishes.
The eigenvector corresponding to a simple eigenvalue λ = 0 (i.e., multiplicity is 1) of the Jacobian matrix Df(B 1 , ξ 0 ) evaluated at a point B 1 when ξ = ξ 0 is where T stands for transpose and f is the vector field of the system (37)- (38). Also, the eigenvector corresponding to a simple eigenvalue λ = 0 of the transpose of the Jacobian matrix is w = [0 1] T .
After, few simple algebraic calculations one could easily verify that when ξ = ξ 0 , we have and where vector f ξ denotes the partial derivative of f with respect to ξ . Note here that the right hand side of (44), (45) is non-zero for any admissible value of n and w, otherwise ξ is undefined or negative. Hence, by the Sotomayor's theorem, the system (37)-(38) undergoes a transcritical bifurcation when ξ = 3 2 1−w 2 n(3w−1)−6 (w+1) , which is represented by a curve D1 of Figs. 2(a), 2(d) for w = 0. In a similar manner, one could also verify the following statements.
for fixed w.
For w = 0, the above equation is represented by two branch curves D4 and D5 of Figs. 2(d), 2(e). The above bifurcation scenarios help us to properly separate the Universe's evolution into generic and non-generic evolution. In Table IV, we present the conditions satisfied by the model parameters to describe the Universe's possible generic cosmological evolution. This type of evolution is interesting as it can determine the initial phase and the final phase of the Universe's cosmological evolution for a wide range of initial conditions. Out of all the scenarios presented in Table IV, only scenarios II and IV are of cosmological interest as they could explain the late-time behavior of the Universe and fits with various CMB and BAO observational data [57]. In Fig.  TABLE IV. Conditions on parameters ξ , n, w for which dynamics of Universe governed by the system (37)-(38) undergoes a generic evolution.
Here, Q ± = 3 wξ +7 ξ ± √ 9 w 2 ξ 2 +66 wξ 2 +12 wξ +73 ξ 2 +12 ξ 2ξ . 3, we present the global phase space diagram for a particular choice of n and ξ which corresponds to a generic evolution of scenario II where the model exhibits a late-time acceleration (similar dynamics is exhibited by scenario IV). It is important to remark here that for the parameters within the range of generic scenario II, there is no bifurcation, so it is sufficient to present only one phase portrait for this scenario. Furthermore, in Sec. V, we will discuss the evolution of cosmological quantities described by this model for parameters corresponding to scenario II.

V. DISCUSSION AND CONCLUSION
In this work, we studied the global qualitative cosmological dynamics of a non-minimal coupled scalar field for a general class of coupling function and potential. We focused on the bifurcation analysis to investigate the effect of varying the model parameters on the global dynamics. The main objective of applying bifurcation theory is to determine the existence of generic evolutionary scenarios. It will help us identify the Universe's initial and final phases for a wide range of initial conditions. The bifurcation theory also allows us to understand how physics described by the phase space change with parameters.
The general model described by the system (20)-(23) contains different interesting cosmological solutions for different model parameters. For instance, the present model exhibits a scalar field dominated solution A 2± , resembling a stiff matter Universe or radiation Universe for some choice of coupling and potential parameters. The scalar field's sole contribution to an early radiation epoch is an interesting scenario missing in a minimal coupled canonical scalar field within the GR context. For some choice of coupling function and potential, the general model also shows the presence of DM-DE scaling solutions A 1 and A 3 describing an intermediate matter epoch. Lastly, the present model exhibits a late time acceleration of the Universe via a critical point A 4 . Thus, from the analysis presented in Sec. III, one can select a class of non-minimal coupled scalar field models describing a viable sequence of cosmic evolution. In particular, for models with λ F * = λ V * 2 , point A 4 corresponds to a late time de Sitter solution of the GR case. Hence, such a class of scalar-tensor theories converges towards a structurally stable GR-based model, i.e., the ΛCDM model. This is a common feature of the scalar-tensor gravity which has been verified numerically and analytically [61][62][63][64][65][66]. Therefore, the present analysis identifies a broad class of scalar-tensor models generically possessing this property. Further, it is possible that the point A 4 corresponds to a late time super-accelerated phase (w eff < −1) when λ F * < λ V * < 2λ F * . Hence, in the presence of a non-minimal coupling, the model can explain a super-accelerated phase (w eff < −1) without the need for the introduction of a phantom field. Depending on the potential and coupling function, point A 4 is a saddle in nature and hence can also describe the possible inflationary exit phenomenon. In Sec. IV, we illustrate the global dynamics in detail by considering a simple example corresponding to a quadratic coupling function and a power-law form potential. For this particular example, we find the range of parameters n and ξ describing a generic evolution scenario (see Table IV). In Fig.  4, we have plotted the evolution of cosmological parameters against the redshift z for a choice of parameters corresponding to one of the generic scenarios, i.e., scenario (II) which starts from an unstable decelerated solution (B 2± ) towards a deSitter like solution B 4 (a similar evolution is obtained for scenario IV). Recall that redshift z = −1 + 1 a where the present value of the scale factor taken to be unity, with z = 0 corresponds to the present Universe and z = −1 corresponds to its infinite future. It is worth mentioning here that only the parameters' values corresponding to the generic evolution of scenarios II or IV satisfy various observational constraints coming from Planck and BAO datasets [58]. Therefore, mathematically, stability analysis and bifurcation theory help us to locate various physically rich models. For instance, we can find the parameter's range which can generically describe the thermal history of the Universe starting from an early radiation domination (w eff 1 3 ) or stiff matter solution (w eff 1) represented by points B 2± towards a DE dominated solution B 4 (w eff −1) via a matter like solution B 1 (w eff 0) i.e., B 2± → B 1 → B 4 (see Figs. 3, 4).
The possibility that Ω m > 1 and w φ diverges at the onset of a matter-dominated era (as explained in paragraphs after Eqs. (13) and (27)) can also be confirmed from Fig. 4(a). The divergence of w φ , however, does not cause any issue to the evolution behavior of w eff , as this corresponds to the vanishing of scalar field energy density. The overdensity of matter component is not a surprise in cosmological models where interaction between different components occurs [67]. Such behavior is more visible by comparing a change in the behavior of Ω m and w φ to ξ (i.e., λ F ) from Figs. 4(a) and 4(c) (the range of divergence of w φ reduces and Ω m → 1 as ξ → 0). Interestingly, the behavior is consistent with the result of [56], that the background dynamics of the present model approach GR with the cosmological constant in the limit ξ → 0. For a generic scenario (II), there is also a possibility of crossing the phantom divide line and eventually the solution settles down towards a cosmological constant behavior (see Fig. 3). We can confirm it by taking a closer look at the late time evolution of w eff (see Fig. 4(b)). In Fig. 4, we choose initial conditions such that the Universe agrees with the current observational data i.e., Ω m ≈ 0.3, w eff ≈ −0.8 [68].
Our analysis shows that the non-minimal coupled scalar field model describes a rich cosmological history of the Universe. For instance, the model exhibits the transition: radiation → matter → DE and the possible crossing of phantom divide line but avoiding the big rip singularity. It is worth noting that only by using the stability analysis, one can determine such a transition. However, bifurcation tools help us locate the parameter values describing such dynamics without a fine-tuning of initial conditions. With the present choice of variables (19), this model can explain the graceful exit scenario and the late time DE era separately. Such a result is a common feature of many classes of scalar-tensor theories. Therefore, it is of inter-est to extend the analysis to the case where the scalar field is coupled non-minimal with gravity and matter, as discussed in [69]. Such discussion is beyond the scope of the present work.
, y → y and we apply to the dynamical system (37)-(38) (which we have not presented here). Then, the resulting system is then transformed into a standard form upon the introduction of new variables X,Y given by On employing these new variables, we can rewrite the corresponding dynamical system as where g 1 , f 1 have not been presented due to their length. Then by the center manifold theory, there exist a continuously differentiable function h : R → R defined by X = h(Y ) = a 2 Y 2 +a 3 Y 3 +O(4), where a 2 , a 3 ∈ R are determined from the quasi-linear equation where A = 0, B = 3 2 n(w−1)+2(w+1) n and D denotes the derivative with respect to Y . On substituting the expression of A, B, f 1 , g 1 and h in (A1), we obtain the values of a 2 and a 3 given by a 2 = 2 6(1−w 2 ) 3 nw−n−6 w−6 (w + 1) (n + 1) (3 nw − n − 6 w − 6) (1 − w) (3 w − 1) 2 (nw − n + 2 w + 2) , The following equation then gives the flow on the corresponding local center manifold i.e., This equation implies point B 1 is always saddle when ξ = 3(1−w 2 ) 2(3 nw−n−6 w−6) .
Appendix B: Center manifold theory analysis for critical point B 4 In this case, we analyze the stability behavior of the point B 4 when ξ = 6 n 2 −8n−20 (i.e., curves D2, D3 of Fig. 2) or ξ = 3(w+1) n 2 −3 nw−7 n−6 w−6 (i.e., curves D4, D5 of Fig. 2). Following the similar analysis as for the point B 1 , we obtained that the flow on the corresponding local center manifold in both the cases is given by

Appendix C: A brief introduction to bifurcation tools
Here, we present preliminaries and two landmark theorems required for understanding bifurcation analysis. For more details, reader can refer to Refs. [18,19,70].
Definition 1 Two dynamical systems are said to be locally topologically equivalent if there exists a homeomorphism (i.e., a continuous invertible function whose inverse is also continuous) mapping orbits of one system onto orbits of another system preserving the direction of time.
If the qualitative behavior remains topologically equivalent for all nearby vector fields, then the system or the vector field is said to be structurally stable. For a two dimensional system, Peixoto's theorem completely characterizes the structural stability of vector fields on a compact, two-dimensional manifold. To understand the landmark theorem, we present a definition of non-wandering points.
Definition 2 A point x on a manifold is a non-wandering point of the flow φ t defined by the vector field if for any neighborhood U of x and for any T > 0 there is t > T such that φ t (U) ∩U is a non-empty set.
Peixoto's Theorem Let f be a C 1 -vector field on a compact, two dimensional, differentiable manifold M. Then f is structurally stable on M if and only if (i) the number of critical points and cycles is finite and each is hyperbolic; (ii) there are no trajectories connecting saddle points; and (iii) the set of all non-wandering points consists of critical points and limit cycles only.
Recall that a limit cycle is an isolated closed path. By isolated, it means that neighboring trajectories are either spiral toward or away from a limit cycle. When the dynamical system depends on some parameters, the system's phase portrait also varies as parameters vary. Thus, either the phase portrait remains topologically equivalent, or its topology changes as parameter changes. The occurrence of topologically inequivalent phase portraits under a change of parameters is called a bifurcation. A parameters' value at which the topology changes is called a bifurcation value.
Different types of bifurcation can occur for a given dynamical system. One can classify different types of bifurcation using Sotomayor's theorem (see [18] for more details). Since, in our work, we obtained only transcritical bifurcation, in what follows, we state this theorem to determine the occurrence of transcritical bifurcation. Sotomayor's theorem for transcritical bifurcation Consider the systemẋ = f(x, µ) where the set of vector fields f equipped with the standard C 1 -norm 1 forms a Banach space such that f(x 0 , µ 0 ) = 0. Suppose the Jacobian matrix (A ≡ Df(x 0 , µ 0 )) has a simple eigenvalue λ = 0 with eigenvector v and w is an eigenvector of the transpose of the Jacobian matrix A T corresponds to the eigenvalue λ = 0. Then the above system experiences a transcritical bifurcation at the equilibrium point x 0 as the parameter µ varies through the bifurcation value µ = µ 0 , if the following three conditions hold: • w T f µ (x 0 , µ 0 ) = 0 and f µ denotes partial derivative of f with respect to µ.