Jacobi stability analysis of scalar field models with minimal coupling to gravity in a cosmological background

We perform the study of the stability of the cosmological scalar field models, by using the Jacobi stability analysis, or the Kosambi-Cartan-Chern (KCC) theory. In the KCC approach we describe the time evolution of the scalar field cosmologies in geometric terms, by performing a"second geometrization", by considering them as paths of a semispray. By introducing a non-linear connection and a Berwald type connection associated to the Friedmann and Klein-Gordon equations, five geometrical invariants can be constructed, with the second invariant giving the Jacobi stability of the cosmological model. We obtain all the relevant geometric quantities, and we formulate the condition of the Jacobi stability for scalar field cosmologies in the second order formalism. As an application of the developed methods we consider the Jacobi stability properties of the scalar fields with exponential and Higgs type potential. We find that the Universe dominated by a scalar field exponential potential is in Jacobi unstable state, while the cosmological evolution in the presence of Higgs fields has alternating stable and unstable phases. By using the standard first order formulation of the cosmological models as dynamical systems we have investigated the stability of the phantom quintessence and tachyonic scalar fields, by lifting the first order system to the tangent bundle. It turns out that in the presence of a power law potential both these models are Jacobi unstable during the entire cosmological evolution.

We perform the study of the stability of the cosmological scalar field models, by using the Jacobi stability analysis, or the Kosambi-Cartan-Chern (KCC) theory. In the KCC approach we describe the time evolution of the scalar field cosmologies in geometric terms, by performing a "second geometrization", by considering them as paths of a semispray. By introducing a non-linear connection and a Berwald type connection associated to the Friedmann and Klein-Gordon equations, five geometrical invariants can be constructed, with the second invariant giving the Jacobi stability of the cosmological model. We obtain all the relevant geometric quantities, and we formulate the condition of the Jacobi stability for scalar field cosmologies in the second order formalism. As an application of the developed methods we consider the Jacobi stability properties of the scalar fields with exponential and Higgs type potential. We find that the Universe dominated by a scalar field exponential potential is in Jacobi unstable state, while the cosmological evolution in the presence of Higgs fields has alternating stable and unstable phases. By using the standard first order formulation of the cosmological models as dynamical systems we have investigated the stability of the phantom quintessence and tachyonic scalar fields, by lifting the first order system to the tangent bundle. It turns out that in the presence of a power law potential both these models are Jacobi unstable during the entire cosmological evolution.

I. INTRODUCTION
A large number of cosmological observations, obtained initially from distant Type Ia Supernovae, have convincingly proven that the Universe has undergone a late time accelerated expansion [1][2][3][4]. In order to explain these observations a deep change in our paradigmatic understanding of the cosmological dynamics is necessary, and many ideas have been put forward to address them. The "standard" explanation of the late time acceleration is based on the assumption of the existence of a mysterious component, called dark energy, which is responsible for the observed characteristics of the late time evolution of the Universe.
On the other hand, the combination of the results of the observations of high redshift supernovae, and of the WMAP and of the recently released Planck data, indicate that the location of the first acoustic peak in the power spectrum of the Cosmic Microwave Background Radia- * Electronic address: bogdan.danila22@gmail.com † Electronic address: t.harko@ucl.ac.uk ‡ Electronic address: mankwongmak@gmail.com § Electronic address: kppraibo@kmitl.ac.th ¶ Electronic address: sorin@tokai.ac.jp tion is consistent with the prediction of the inflationary model for the density parameter Ω, according to which Ω = 1. The cosmological observations also provide strong evidence for the behavior of the parameter w = p/ρ of the equation of state of the cosmological fluid, where p is the pressure and ρ is the density, as lying in the range −1 ≤ w = p/ρ < −1/3 [5].
In order to explain the observed cosmological dynamics, it is assumed usually that the Universe is dominated by two main components: cold (pressureless) dark matter (CDM), and dark energy (DE) with negative pressure, respectively. CDM contributes Ω m ∼ 0.3 [6], and its introduction is mainly motivated by the necessity of theoretically explaining the galactic rotation curves and the large scale structure formation. On the other hand, DE represents the major component of the Universe, providing Ω DE ∼ 0.7. Dark energy is the major factor determining the recent acceleration of the Universe, as observed from the study of the distant type Ia supernovae [5]. Explaining the nature and properties of dark energy has become one of the most active fields of research in cosmology and theoretical physics, with a huge number of proposed DE models( for reviews see, for instance, [7][8][9][10][11]).
One interesting possibility for explaining DE are cosmological models containing a mixture of cold dark matter and quintessence, representing a slowly-varying, spatially inhomogeneous component [12]. From a theoret-ical as well as a particle physics point of view the idea of quintessence can be implemented by assuming that it is the energy associated with a scalar field Q, having a self-interaction potential V (Q). When the potential energy density V (Q) of the quintessence field is greater than the kinetic one, then it follows that the pressure p =Q 2 /2−V (Q) associated to the quintessence Q-field is negative. The properties of the quintessential cosmological models have been actively considered in the physical literature (for a recent review see [13]). As opposed to the cosmological constant of standard general relativity, the equation of state of the quintessence field changes dynamically with time [14]. Alternative models, in which the late-time acceleration can be driven by the kinetic energy of the scalar field, called k−essence models, have also been proposed [15].
Scalar fields φ that are minimally coupled to gravity via a negative kinetic energy, can also explain the recent acceleration of the Universe. Interestingly enough, they allow values of the parameter of the equation of state with w < −1. These types of scalar fields, known as phantom fields, have been proposed in [16]. For phantom scalar fields the energy density and pressure are given by ρ φ = −φ 2 /2 + V (φ) and p φ = −φ 2 /2 − V (φ), respectively. The interesting properties of phantom cosmological models for dark energy have been investigated in detail in [17][18][19]. Recent cosmological observations show that at some moment during the cosmological evolution the value of the parameter w may have crossed the standard value w = −1, corresponding to the general relativistic cosmological constant Λ. This cosmological situation is called the phantom divide line crossing [19]. In the case of scalar field models with cusped potentials, the crossing of the phantom divide line was investigated in [18]. Another alternative way of explaining the phantom divide line crossing is to model dark energy by a scalar field, which is non-minimally coupled to gravity [18]. Scalar fields are also assumed to play a fundamental role in the evolution of the very early Universe, playing a major role in the inflationary scenario [20,21]. Originally, the idea of inflation was proposed to provide solutions to the singularity, flat space, horizon, and homogeneity problems, to the absence of magnetic monopoles, as well as to the problem of large numbers of particles [22,23]. However, presently it is believed that the most important feature of inflation is the generation of both initial density perturbations and of the background cosmological gravitational waves. These important cosmological parameters can be determined in many different ways, like, for example, through the study of the anisotropies of the microwave background radiation, the analysis of the local (peculiar) velocity galactic flows, of the clustering of galaxies, and the determination of the abundance of gravitationally bound structures of different types, respectively [23].
In many inflationary models the dynamical evolution of the early Universe is driven by a single scalar field, called the inflaton, with the inflaton rolling in some underlying self-interaction potential [20][21][22][23]. One common approximation in the study of the inflationary evolution is the slow-roll approximation, which can be successfully used in two separate contexts. The first situation is in the study of the classical inflationary dynamics of expansion in the lowest order approximation. Hence this implies that the contribution of the kinetic energy of the inflaton field to the expansion rate is ignored. The second situation is represented by the calculation of the perturbation spectra. The standard expressions deduced for these spectra are valid in the lowest order in the slow roll approximation [24]. Finding exact inflationary solutions of the gravitational field equations for different types of scalar field potentials is also of great importance for the understanding of the dynamics of the early Universe. Such exact solutions have been found for a large number of inflationary potentials. Moreover, the potentials allowing a graceful exit from inflation have been classified [25].
Hence, the theoretical investigation of the scalar field models is an essential task in cosmology. Among the various methods used to study the properties of scalar fields the methods based on the applications of the mathematical formalism of the qualitative study of dynamical systems is of considerable importance.
The usefulness of dynamical systems formulation of physical models is mainly determined by their powerful predictive power. This predictive power is essentially determined by the stability of their solutions. In a realistic physical system, due to the limited precision of the measurements, some uncertainties in the initial conditions always exist. Therefore a physically meaningful mathematical model must also offer detailed and useful information on the evolution of the deviations of the possible trajectories of the dynamical system from a given reference trajectory. Hence an important requirement in mathematical modelling is the understanding of the local stability of the physical and cosmological processes. This information on the system behavior is as important as the understanding of the late-time deviations. The global stability of the solutions of the dynamical systems described by systems of non-linear ordinary differential equations is analyzed in the framework of the well-known mathematical theory of Lyapounov stability. In this mathematical approach the fundamental quantities that measure exponential deviations from a given trajectory are the so-called Lyapunov exponents [26,27]. It is usually very difficult to determine the Lyapounov exponents analytically for a given dynamical system, and thus one must resort to numerical methods. On the other hand, the important problem of the local stability of the solutions of dynamical systems, described by ordinary differential equations, is less understood.
Cosmological models have been intensively investigated by using methods from dynamical systems and Lyapounov stability theory [28][29][30][31][32]. In particular, phase space analysis proved to be a very useful method for the understanding of the cosmological evolution. When studying the evolution of cosmological models, the dynamical equations can be represented by an autonomous dynamical system, described by a set of coupled -usually strongly non-linear -differential equations for the physical parameters. This representation allows the study of the Lyapunov stability of the model, without explicitly solving the field equations for the basic variables. Furthermore, the importance of the Lyapunov analysis is related to the fact that stationary points of the dynamical system correspond to exact or approximate analytic solutions of the field equations. Thus the dynamical system formulation provide a useful tool for obtaining exact or approximate solutions of the field equations in cosmologically interesting situations.
Even that the mathematical methods of the Lyapounov stability analysis are well established, the study of the stability of the dynamical systems from different points of view is extremely important. The comparison of the results of the alternative approach with the corresponding Lyapunov exponents analysis can provide a deeper understanding of the stability properties of the system. An alternative, and very powerful method for the study of the systems of the ordinary differential equations is represented by the so-called Kosambi-Cartan-Chern (KCC) theory, which was initiated in the pioneering works of Kosambi [33], Cartan [34] and Chern [35], respectively. The KCC theory was inspired and influenced by the geometry of the Finsler spaces (for a recent review of the KCC theory see [36]). From a mathematical point of view the KCC theory is a differential geometric theory of the variational equations for the deviations of the whole trajectory with respect to the nearby ones [37]. In the KCC geometrical description of the systems of ordinary differential equations one associates a non-linear connection, and a Berwald type connection to the system of equations. With the use of these geometric quantities five geometrical invariants are obtained. The most important invariant is the second invariant, also called the curvature deviation tensor, which gives the Jacobi stability of the system [36][37][38][39]. The KCC theory has been applied for the study of different physical, biochemical or technical systems (see [38][39][40][41][42].
An alternative geometrization method for dynamical systems, with applications in classical mechanics and general relativity, was proposed in [43] and [44], and further investigated in [45,46]. The Henon-Heiles system and Bianchi type IX cosmological models were also investigated within this framework. In particular, in [45] a theoretical approach based on the geometrical description of dynamical systems and of their chaotic properties was developed. For the base manifold a Finsler space was introduced, whose properties allow the description of a wide class of physical systems, including those with potentials depending on time and velocities, for which the Riemannian approach is unsuitable.
It is the purpose of the present paper to consider a systematic investigation of the Jacobi stability properties of the flat homogeneous and isotropic general relativistic cosmological models. By starting from the standard Friedmann equations we perform, as a first step in our analysis, a "second geometrization" of these equations, by associating to them a non-linear connection, and a Berwald connection, respectively. This procedure allows to obtain the so-called KCC invariants of the Friedmann equations. The second invariant, called the curvature deviation tensor, gives the Jacobi stability properties of the cosmological model. The KCC theory can be naturally applied to systems of second order ordinary differential equations. The Friedmann equations can be formulated as second order differential equations, similarly to the Klein-Gordon equation describing the scalar field. Therefore the KCC theory can be applied to matter and scalar field dominated cosmological models. We obtain the general condition for the Jacobi stability of scalar fields, which is described by two inequalities involving the second and the first derivative of the scalar field potential, the energy density of the field, as well as the time derivative of the field itself. The geodesic deviation equations describing the time variation of the deviation vector are also obtained. As an application of the developed formalism we investigate the stability properties of the scalar field cosmological models with exponential and Higgs type potentials, respectively. It turns out that the exponential potential scalar field is Jacobi unstable during its entire evolution, while the time evolution of the scalar field cosmological models with Higgs potential show a complicated dynamics with alternating stable and unstable Jacobi phases. The Jacobi stability properties of the Higgs type models are determined by the numerical value of the ratio of the self-coupling constant and the square of the mass of the Higgs particle.
The Lyapunov stability properties of the scalar field cosmological models are usually investigated by reformulating the evolution equation as a set of three first order ordinary differential equations. In order to apply the KCC theory to such systems they must be lifted to the tangent bundle. From mathematical point of view this requires to take the time derivative of the first order equations, so that their "second geometrization" can be easily performed. We consider in detail the Jacobi stability properties of the phantom quintessence and tachyon scalar field cosmological models. We study in detail the Jacobi stability condition of these models, and we find that they are Jacobi unstable during the entire expansionary cosmological evolution.
The present paper is organized as follows. We review the basic ideas and the mathematical formalism of the KCC theory in Section II. The Jacobi stability analysis of the homogenous isotropic flat cosmological models by using the second order formulation of the dynamics is performed in Section III. We consider both the cases of the matter dominated and scalar field dominated cosmological models. As an application of the developed formalism we investigate in detail the Jacobi stability of the scalar fields with exponential potential, and Higgs poten-tial, respectively. The Jacobi stability of the first order dynamical system formulation of scalar field cosmological models is considered in Section refsect4, in which the KCC geometrization of the phantom quintessence and tachyonic scalar field models is analyzed in detail. We discuss and conclude our results in Section V. The KCC geometric quantities giving the geometric description of the phantom quintessence and tachyon scalar field cosmologies are presented in Appendix A and Appendix B, respectively.

II. KOSAMBI-CARTAN-CHERN (KCC) THEORY AND JACOBI STABILITY
In the present Section we briefly summarize the basic concepts and results of the KCC theory (for a detailed presentation see [36] and [37]).

A. Dynamical systems as paths of a semispray
In the following we denote by M a real, smooth ndimensional manifold, and by T M its tangent bundle. On an open connected subset Ω of the Euclidian (2n + 1) dimensional space R n × R n × R 1 we introduce a 2n + 1 dimensional coordinates system x i , y i , t , i = 1, 2, ..., n, where x i = x 1 , x 2 , ..., x n , y i = y 1 , y 2 , ..., y n and t is the time t. The coordinates y i are defined as We assume that the time t is an absolute invariant, and therefore the only admissible change of coordinates will bẽ Definition [47]. A deterministic dynamical systems is a formal set of rules that describe the evolution of points in some set S with respect to an external, discrete, or continuous time parameter t running in another set T .
In a rigorous mathematical formulation, a dynamical system is a map [47] φ : which satisfies the condition φ(t, ·)•φ(s, ·) = φ(t+s, ·) for all times t, s ∈ T . In order to model realistic dynamical systems or physical processes additional structures must be added to the above definition. In many cases the equations of motion of a dynamical system can be derived from a Lagrangian L via the Euler-Lagrange equations, where F i , i = 1, 2, ..., n, is the external force. For a regular Lagrangian L, the Euler-Lagrange equations given by Eq. (4) are equivalent to a system of second-order differential equations [48] where each function is called a semispray [49][50][51]. The functions G i (x i , y i ) are the local coefficients of the semispray, and they are defined on domains of local charts. In the particular case in which the coefficients Conversely, for any system of ordinary differential equations of the form (7), which is globally defined, the functions G i define a semispray on T M [50,51].
More generally, one can start from an arbitrary system of second-order differential equations of the form (5), where no a priori given Lagrangian function is assumed, and study the behavior of its trajectories by analogy with the trajectories of the Euler-Lagrange equations.

B. The KCC geometrization of dynamical systems
To associate a geometrical structure to the dynamical system defined by Eqs. (5), we introduce first a nonlinear connection N on M , with coefficients N i j , defined as [48] The nonlinear connection can be understood in terms of a dynamical covariant derivative ∇ N [47]: for two vector fields v, w defined over M, we introduce the covariant derivative ∇ N as For N i j (x, y) = Γ i jl (x)y l , Eq. (9) reduces to the definition of the covariant derivative for the special case of a linear connection.
For a non-singular coordinate transformations given by Eq. (2), the KCC-covariant differential of an arbitrary vector field θ i (x) on the open subset Ω ⊆ R n × R n × R 1 is defined as [36,38] For θ i = y i we obtain where the contravariant vector field ǫ i on Ω is called the first KCC invariant.

The curvature deviation tensor
Let us now vary the trajectories x i (t) of the system (5) into nearby ones according tõ where |η| is a small parameter and ξ i (t) are the components of some contravariant vector field defined along the path x i (t). Substituting Eqs. (12) into Eqs. (5) and taking the limit η → 0 we obtain the variational equations [37][38][39][40] By using the KCC-covariant differential we can write Eq. (13) in the covariant form where we have denoted and we have introduced the Berwald connection G i jl , defined as [36][37][38][39][40]48] The tensor P i j is called the second KCC-invariant, or the deviation curvature tensor, and Eq. (14) is called the Jacobi equation, respectively. In both Riemann or Finsler geometry, when the system (5) describes the geodesic equations in the given geometry, Eq. (14) represents the Jacobi field equation.
The trace P of the curvature deviation tensor is obtained as The third, fourth and fifth invariants of the system (5) are defined as [37] The third invariant P i jk can be interpreted as a torsion tensor. The fourth and fifth invariants P i jkl and D i jkl are called the Riemann-Christoffel curvature tensor, and the Douglas tensor, respectively [36,37]. In Berwald spaces these tensors always exist. Generally, they can be used to describe the geometrical properties of systems of secondorder differential equations.

C. The definition of Jacobi stability
In many physical applications the behavior of the trajectories of the dynamical system (5) in a vicinity of a point x i (t 0 ) is of major interest. In the following for simplicity we take t 0 = 0. The trajectories x i = x i (t) can be considered as curves in the Euclidean space (R n , ., . ), where ., . defines the canonical inner product of R n . As for the deviation vector ξ we assume that it satisfies the natural initial conditions We describe the focusing tendency of the trajectories around t 0 = 0 as follows: if ||ξ (t)|| < t 2 , t ≈ 0 + , the trajectories are bunching together. On the other hand, if ||ξ (t)|| > t 2 , t ≈ 0 + , the trajectories are dispersing [36][37][38][39]. Alternatively, the focusing tendency of the trajectories can be described in terms of the deviation curvature tensor: the trajectories of the system of second order differential equations (5) are bunching together for t ≈ 0 + if and only if the real parts of the eigenvalues of the deviation curvature tensor P i j (0) are strictly negative. On the other hand the trajectories are dispersing if and only if the real parts of the eigenvalues of the tensor P i j (0) are strictly positive [36][37][38][39].
Based on the above intuitive considerations we can define the Jacobi stability for a second order system of differential equations as follows [36][37][38][39]: Definition: If the system of second order differential equations (5) satisfies the initial conditions with respect to the norm ||.|| induced by a positive definite inner product, then we call the trajectories of (5) as Jacobi stable if and only if the real parts of the eigenvalues of the deviation tensor P i j are strictly negative everywhere. Otherwise, the trajectories are called Jacobi unstable.
The focussing/dispersing behavior of the trajectories of a system of second order ordinary differential equations near the origin is represented in Fig. 1. In the important two dimensional case the curvature deviation tensor can be written in a matrix form as with the eigenvalues given by The eigenvalues of the curvature deviation tensor can be obtained as solutions of the quadratic equation A very powerful algebraic method to obtain the signs of the eigenvalues of the curvature deviation tensor is represented by the Routh-Hurwitz criteria [52]. According to these criteria all of the roots of the polynomial P (λ) are negatives or have negative real parts if the determinants of all Hurwitz matrices det H j , j = 1, 2, .., n, are strictly positive. For n = 2, corresponding to the case of Eq. (21), the Routh-Hurwitz criteria takes the simple form The curvature properties along a given geodesic are described by the eigenvalues of the deviation curvature tensor λ ± , which are invariant functions on the tangent space. Moreover, once they are known, we can introduce two quantities that can characterize the way the geodesic explores the base manifold. They are defined as the (half) of the Ricci curvature scalar along the flow, κ, and the anisotropy θ, given by and respectively.
In the case of a three-dimensional dynamical system with n = 3, the characteristic equation of the matrix of the curvature deviation tensor becomes Therefore the conditions of the Jacobi stability for a three-dimensional system of second order ordinary differential equations can be formulated as follows:

III. JACOBI STABILITY ANALYSIS OF ISOTROPIC MATTER DOMINATED AND SCALAR FIELD COSMOLOGIES
In the present Section we use the KCC approach for the study of the dynamical properties of the matter dominated and scalar field cosmologies. We explicitly obtain the non-linear and Berwald connections, and the deviation curvature tensors. The eigenvalues of the deviation curvature tensor are also obtained, and we study their properties in the equilibrium points of the matter dominated model. The study of the sign of the eigenvalues allows us to obtain the Jacobi stability properties of the fixed points of the matter field dominated cosmological models. Next, we proceed to a detailed analysis of the scalar field cosmologies in the framework of the KCC theory. A full "second geometric" description is introduced, and the time evolution of the relevant physical and geometrical parameters is obtained. In particular, the nature of the Jacobi stability is analyzed in detail. In the present study we restrict our analysis to the case of the flat Friedmann-Robertson-Walker metric, given by where a is the scale factor.
A. Jacobi stability of matter dominated cosmological models For a Universe filled with pressureless dust and radiation only, the cosmological expansion is described mathematically by the Friedmann equations, which take the well-known form where a dot denotes the derivative with respect to the time t. In Eqs. (31) and (32) H =ȧ/a denotes the Hubble function, ρ m is the baryonic matter energy density, ρ r is the energy density of the radiation, while Λ is the cosmological constant.

Friedmann equations as an autonomous dynamical system
In order to reformulate the cosmological evolutions equations as a dynamical system, we need to introduce first the density parameters (Ω m , Ω r , Ω Λ ) of the matter, radiation and cosmological constant, defined as The density parameters satisfy the normalization relation As the basic variables (x, y) in the phase space we adopt the quantities x ≡ Ω r , and y ≡ Ω Λ , respectively [31]. Then the density parameter of the mater is given by Ω m = 1 − x − y, and the range of the variables (x, y, Ω m ) is 0 ≤ x ≤ 1, 0 ≤ y ≤ 1, and 0 ≤ Ω m ≤ 1. To describe the cosmological dynamics we define the physically significant phase space as Φ = {(x, y) : x + y ≤ 1, 0 ≤ x ≤ 1, 0 ≤ y ≤ 1}. Next we take the time derivatives of x and y with respect to the time t, and, after introducing the new time variable τ = ln a(t), we can formulate the Friedmann equations as an autonomous dynamical system given by [31] dx dτ The critical points of the system (35) and (36) in the phase space region defined by Ψ are obtained by solving the algebraic equations x(1 − x + 3y) = 0 and (3 + x − 3y)y = 0, respectively, and are given by The Lyapunov stability properties of the system follows from the study of the Jacobian matrix The critical point of the system (35) and (36) have a clear cosmological interpretation. Thus, the critical point P dS = (0, 1), with Ω Λ = 1, has a(t) ∝ exp Λ/3t , and is associated to an accelerated de Sitter type expansion, being a future attractor [32]. The critical point P r = (1, 0), with Ω r = 1 and a(t) ∝ √ t, corresponds to the radiation-dominated era in the cosmological evolution of the Universe, and is a source point, or a past attractor. Finally, the critical point P m (0, 0), with Ω m = 1 and a(t) ∝ t 2/3 , corresponds to the decelerating matterdominated phase of the cosmological expansion. It turns out that P m is a saddle critical point [32].

The KCC geometrization and the Jacobi stability of the Friedmann equations
In order to apply the KCC theory to the cosmological dynamical system given by equations (35) and (36), we first relabel the variables as x ≡ x 1 and y ≡ x 2 . We also denote y 1 = dx 1 /dτ and y 2 = dx 2 /dτ , respectively. Hence we obtain Next, we take the derivative with respect to τ of Eqs. (39) and (40), respectively, thus obtaining the following lift on the tangent bundle of the cosmological dynamical system, By comparison with Eqs. (5) we obtain immediately Hence the components of the non-linear connection (N ) = N j i associated to the matter dominated cosmological dynamical system in the presence of a cosmological constant are obtained as For the components of the deviation tensor P = P j i we obtain where is the Hessian of f , and H g is the Hessian of g. Explicitly, the curvature deviation tensor for the Friedmann cosmological dynamical system can be obtained as By evaluating P at the critical points gives where

Jacobi stability of the critical points of the matter dominated cosmological models
In order to obtain the Jacobi stability of the critical points of the cosmological model described by Eqs. (35) and (36), we need to compute the numerical values of the curvature deviation tensor at the critical points. At the critical point P dS = x 1 = 0, x 2 = 1 the curvature deviation tensor takes the form and has the eigenvalues λ 1 = 4 and λ 2 = 9. Hence it follows that the critical point P dS of the standard ΛCDM cosmological model is Jacobi unstable. For the critical point P r = x 1 = 1, x 2 = 0 , we find with the corresponding eigenvalues of the curvature deviation tensor obtained as λ 1 = 16 and λ 2 = 1/4. Hence, the critical point P r is also Jacobi unstable. Finally, for the last critical point P m = (0, 0) we obtain with the eigenvalues λ 1 = 9, λ 2 = 1/4. Hence, we obtain the result that all the critical points of the dynamical system equivalent to the cosmological Friedmann equations are Jacobi unstable.

B. The cosmological evolution equations in the presence of minimally coupled scalar fields
Let us consider a rather general class of scalar field models, minimally coupled to the gravitational field, for which the Lagrangian density in the Einstein frame reads where R is the curvature scalar, φ is the scalar field, V (φ) is the self-interaction potential and κ = 8πG/c 4 is the gravitational coupling constant, respectively. In the following, we use natural units with c = 8πG = h = 1, and we adopt as our signature for the metric (+1, −1, −1, −1), as is common in particle physics. For a flat FRW scalar field dominated Universe the evolution of a cosmological model is determined by the system of the field equations and the evolution equation for the scalar field where the overdot denotes the derivative with respect to the time-coordinate t, and the prime denotes the derivative with respect to the scalar field φ, respectively. By substitutingȧ/a from Eq. (53) into Eqs. (54) and (55) we can reformulate the dynamics of the scalar field cosmological models in terms of two second order non-linear ordinary differential equations, given bÿ andφ respectively.

The non-linear and Berwald connections, and the KCC invariants of the scalar field cosmological models
In the following we introduce a new notation for the dependent variables a and φ, and for their time derivatives, respectively, as The cosmological dynamics of scalar field dominated Universes can be formulated as a second order differential system, given by two second order differential equations of the form From Eqs. (56) and (57) it follows immediately that and respectively. Therefore we first obtain the components of the non-linear connection as For the non-zero components of the Berwald connection, defined as we obtain The components of the first KCC invariant of the minimally coupled scalar field model are obtained as and respectively. The components of the curvature deviation tensor for minimally coupled scalar field cosmological models are given by For the trace of the curvature deviation tensor we obtain while for χ = P 1 1 P 2 2 − P 1 2 P 2 1 = P 1 1 P 2 2 we have Therefore if the conditions P 1 1 + P 2 2 < 0 and P 1 1 P 2 2 − P 1 2 P 2 1 > 0 are simultaneously satisfied, the scalar field cosmological model is Jacobi stable. These conditions allows us to formulate the following Jacobi stability condition of isotropic and homogeneous scalar field cosmological models. If the parameters V (φ), ρ φ , p φ ,φ of a homogeneous scalar field in an isotropic flat FRW geometry simultaneously satisfy the conditions the corresponding cosmological model is Jacobi stable, and Jacobi unstable otherwise. For the variational differential equations determining the deviation vector ξ, we obtain and respectively. In the case of the isotropic cosmological scalar field models the third, fourth and fifth KCC invariants, as defined by Eqs. (18), are identically equal to zero.

Applications: the case of the scalar field with exponential self-interaction potential
As an application of the KCC geometrization of the scalar field cosmological models we will consider the case of a scalar field with an exponential self-interaction potential of the form where V 0 and λ are constants. The cosmological equations describing the time evolution of this scalar field model areä respectively. By introducing a new time variable τ = √ V 0 t, it turns out that the system of Eqs. (81) and (82) can be written as The system of Eqs. (83)-(84) must be integrated with the initial conditions a(0) = a 0 ,ȧ(0) = a 0 H 0 , φ (0) = φ 0 , andφ (0) =φ 0 , respectively.
The variations of the scale factor and of the scalar field of the exponential potential scalar field filled Universe are presented, for different values of the parameter λ, in Fig. 2. In order to numerically solve the field equations we have used the initial conditions a(0) = 0, φ(0) = 0.001,ȧ(0) = 1, andφ(0) = −0.001, respectively.
As one can see from the figures, the scale factor is a monotonically increasing function of time, while the scalar field also increases during the cosmological evolution. The time variation of the scalar field potential is presented in Fig. 3.
The scalar field potential is a monotonically decreasing function of time, which tends, in the large time limit, to zero. The time variations of the KCC invariant 2κ = P 1 1 +P 2 2 and χ = P 1 1 P 2 2 −P 1 2 P 1 2 are represented in Figs. 4. In order for the considered model of the Universe be stable, the invariants must simultaneously satisfy the conditions 2κ < 0 and χ > 0, respectively. As one can see from the Figures, from the chosen set of parameters, these conditions are not satisfied for any interval of time during the cosmological evolution. Therefore it follows that during its entire evolution an exponential potential scalar field Universe is in a Jacobi unstable state. Finally, in Figs. 5 we present the time variation of the components of the deviation vector ξ i , i = 1, 2.
The component ξ 1 of the deviation vector increases exponentially in time, indicating that the trajectories do diverge exponentially near the origin. Thus, this result also confirms the presence of the Jacobi instability for the exponential potential scalar field cosmology.

Scalar fields with Higgs potential
As a next case of the investigation of the Jacobi stability of a cosmological scalar field model we consider that the Universe is filled with a Higgs-like field, with self-interaction potential given by  where V 0 is a constant, and M 2 < 0 is related to the mass of the Higgs boson by the relation m H = V ′′2 , where v 2 = −M 2 /λ gives the minimum of the potential. The Higgs self-coupling constant λ ≈ 1/8 [53], a value inferred based on the determination of m H from accelerator experiments. By introducing a new dimensionless time variable τ = M t, the basic equations determining the cosmological evolution are given by where v 0 = V 0 /M 2 and η = λ/4M 2 > 0, respectively. The time variations of the scale factor of the Higgs field filled Universe and of the scalar field are represented in Fig. 6. To numerically integrate the gravitational field equations, we have fixed the value of the constant v 0 to v 0 = 0.005, and we have adopted as initial conditions the values the initial conditions a(0) = a 0 ,ȧ(0) = a 0 H 0 , φ (0) = φ 0 , andφ (0) =φ 0 , respectively.
The Universe filled with a Higgs type scalar field is expanding, with the scale factor monotonically increasing in time. In the initial phases of the expansion the scale factor can be approximated by a linearly increasing function of time. The Higgs scalar field φ keeps a constant value in the initial stages of the evolution, followed by a rapid increase associated with an oscillatory behavior with a decreasing amplitude, associated with the energy dissipation. The variation of the Higgs potential is represented in Fig. 7.
After an initial phase in which the potential is constant, rapid oscillations with decreasing amplitude do follow, and the potential decays in time. The variation of the KCC invariants 2κ and χ is represented in Fig. 8.
Jacobi stability of the cosmological model with Higgs type scalar field requires the conditions 2κ < 0 and χ > 0 to be simultaneously satisfied. The Jacobi stable/unstable regions formed during the cosmological expansion of scalar field cosmologies with Higgs selfinteracting potential are presented in Fig. 9.
In Figs. 10 we present the time variation of the components of the deviation vector ξ i , i = 1, 2.

IV. JACOBI STABILITY ANALYSIS OF THE PHANTOM QUINTESSENCE AND TACHYONIC COSMOLOGICAL MODELS
In the present Section we will consider the Jacobi stability analysis of the standard dynamical system formulation of scalar field cosmological models. In this formulation the Friedmann equations are rewritten in the equivalent form of a three-dimensional first order dynam-  ical system. Geometrically, we can describe the solutions of a first order dynamical system as a flow ϕ t : D ⊂ R n → R n , or, more generally, ϕ t : D ⊂ M → M, where M is a smooth n-dimensional manifold. The canonical lift of ϕ t to the tangent space T M is given byφ t : ∂t . Therefore, in order to apply the KCC theory to first order dynamical systems we must first lift the equations to the tangent bundle. Mathematically, this is equivalent to simply take the time derivative of the dynamical system.
A. First order dynamical system formulation of quintessence and phantom quintessence scalar field cosmological models In the following we assume that the energy density and pressure of the scalar field can be generally represented as where ζ = ±1. The case ζ = +1 corresponds to the quintessence fields, while ζ = −1 describes the phantom scalar fields. We assume that the Universe is filled by ordinary matter with energy density ρ m and pressure p m , and scalar fields. Each component of the model is characterized by their dimensionless density parameters Ω m and Ω φ , defined as and satisfying the constraint Ω m + Ω φ = 1. The cosmological equations describing this Universe model take the form respectively. In the following we assume that the matter obeys a linear barotropic equation of state of the form p m = (γ m − 1) ρ m , where γ m is constant, and 1 ≤ γ m ≤

From the second Friedmann equation we immediately obtain the relation
As basic variables in the first order dynamical description of cosmological dynamics we introduce the quantities x and y, defined as In these variables the density parameter of the scalar field is given by Ω φ = x 2 +y 2 , while the density parameter of the matter can be written as Ω m = 1 − x 2 + y 2 . Moreover, the energy density and pressure of the scalar field are p φ = 3H 2 ζx 2 + y 2 and p φ = 3H 2 ζx 2 − y 2 , respectively. Instead of the ordinary time variable t we introduce the new independent variable τ , defined as τ = ln a(t), giving d/dt = H (d/dτ ). Then, by taking the time derivative of x and y we obtain after some simple transformations [32] dx dτ The above system is not a closed system since one more equation involving V ′ /V is still lacking. Therefore we introduce a third variable z, defined as z = V ′ /V . By taking its derivative with respect to the time τ we obtain where Γ = V V ′′ /V ′2 . In order to lift this dynamical system to the tangent bundle we relabel the coordinates as (97) Then, by taking the derivatives with respect to τ of the first order cosmological dynamical system we obtain the equivalent second order system The KCC geometric quantities (non-linear connection, deviation curvature tensor, and geodesic deviation equations) describing the Jacobi stability properties of the phantom quintessence cosmological model are presented in Appendix A.

B. The phantom scalar field with power law potential
We assume that the potential of the phantom quintessence scalar field is given by a power law function, so that V (φ) = V 0 φ α , where α is a constant. Then it follows immediately that z = V ′ /V = α/φ, Γ = (α − 1)/α = constant, and Γ ′ (z) = 0. An important cosmological indicator is represented by the parameter of the total equation of state of the matter w, which for dust is defined as The variations of the density parameter of the phantom quintessence field Ω φ and of the parameter of the total equation of state are represented, as a function of the dust (γ m = 1) cosmological matter density parameter in Fig. 11. The initial conditions used to integrate the cosmological dynamical system are x(0) = 0.28, y(0) = 0.47, and z(0) = 0.1, respectively.
In this model the Universe starts its evolution in the large time limit from a matter dominated phase, with Ω m = 1. During the cosmological expansion the role of the scalar field becomes dominant, and the Universe reaches the present day in a state of accelerated expansion, with Ω φ = 0.75, and Ω m = 0.25, respectively. On the other hand, the de Sitter phase with w = −1 is reached only for vanishing ordinary matter density, when the Universe is fully dominated by the phantom quintessence field. It is also important to note that the cosmological evolution is basically independent on the numerical values of the exponent α in the scalar field potential.
The conditions of the Jacobi stability of the phantom quintessence cosmological model in its three-dimensional dynamic system representation are given by the four conditions that must be satisfied by the quantities (Σ, Φ, Ψ, Ω), which are functions of the components of the deviation curvature tensor, and which are presented in Eqs. (26)- (29), for different values of α. The time variations of (Σ, Φ, Ψ, Ω) are represented in Fig. 12.
As one can see from Fig. (12), the Jacobi stability condition Ω > 0 is not satisfied during the entire time evolution of this cosmological model. Therefore phantom quintessence scalar field cosmological models with power law potential are Jacobi unstable for all time intervals.
C. Jacobi stability of tachyon field cosmological models Tachyon scalar fields have been proposed as possible candidates to explain both inflation and the late acceleration of the Universe [54,55]. It is also possible for the tachyonic field to trigger an the inflationary expansion, and at a later time generate a non-relativistic fluid, which could account for the existence of dark matter. A tachyonic scalar field is described by the Lagrangian L = −V (φ) 1 −φ 2 , and it can be described in terms of an effective energy density and pressure, given by The Friedmann equations are given by  The parameter of the equation of state of the tachyon field is given by By adopting again for the matter the linear barotropic equation of state p m = (γ m − 1) ρ m , 1 ≤ γ m ≤ 2, from the Friedmann equation we obtain first To formulate the cosmological evolution equation as a dynamical system we introduce a set of variables (τ, x, y, z) defined as In these variables the density parameters of the scalar field and of the matter are given by Then, by taking the derivative of x and y with respect to τ we obtain where Γ(z) = (2/3)V V ′′ /V ′2 . By lifting the cosmological field equations to the tangent bundle we obtain The KCC geometric quantities (non-linear connection, deviation curvature tensor, and geodesic deviation equations) describing the Jacobi stability properties of the tachyon scalar field cosmological model are presented in Appendix B.

Power law potential tachyonic scalar field cosmological models
In the following we assume again that the potential of the tachyonic scalar field is power law type, with V (φ) = V 0 φ α , where V 0 and α are constants. This choice immediately gives Γ = (2/3) (α − 1) /α. The parametric dependence of the density parameter Ω φ of the tachyonic scalar field with power law potential on the matter energy density Ω m is represented in the left panel of Fig. 13. The parametric variation of the matter density parameter Ω m as a function of the parameter of the total equation of state of matter w is shown in the right panel of Fig. 13.
In order to numerically integrate the dynamical system corresponding to the tachyon scalar field model we have used the initial conditions x(0) = 0.28, y(0) = 0.45, and z(0) = 0.1, respectively, and we have assumed that the ordinary matter in the Universe is in the form of zero pressure dust, with γ m = 1. There is a linear relation between Ω φ and Ω m . The Universe starts its cosmological evolution in a matter dominated phase, with Ω m = 1, Ω φ = 0, with decelerating expansion. The energy density of the tachyon field increases in time, and the Universe enters in a de Sitter phase, with the matter density becoming negligibly small. The parameter of the total equation of state satisfies the condition w < 0 during the entire cosmological evolution. It is also interesting to note that the time evolution of this model is basically independent on the numerical values of the exponent α in the power law potential. The time variations of (Σ, Φ, Ψ, Ω) describing the Jacobi stability of the tachyon scalar field cosmological model are represented in Fig. 14.
As one can see from Fig. (14), the Jacobi stability condition Ω > 0 is not satisfied during the entire cosmological evolution period of this scalar field model. Therefore it follows that tachyon scalar field cosmological models with power law potential are Jacobi unstable for all time intervals. This result is independent on the numerical values of the exponent α in the tachyon scalar field potential.

V. DISCUSSIONS AND FINAL REMARKS
In the present paper we have investigated the Jacobi stability properties of the scalar field cosmological models by using the KCC theory, which represents a powerful mathematical method for the analysis of dynamical systems. Scalar field cosmological models represent a nontrivial testing object for studying non-linear effects in the framework of general relativity. From a mathematical point of view the Jacobi (in)stability represents a natural generalization of the (in)stability of the geodesic flow on a differentiable manifold, endowed with a Riemannian or Finslerian type metric to a non-metric setting. The KCC theory can be applied to scalar field cosmological models that can be formulated mathematically as sets of second order ordinary non-linear differential equations. Then the geometric invariants associated to this system (nonlinear and Berwald connections), and the deviation curvature tensor, as well as its eigenvalues, can be explicitly obtained. The time evolution of the components of the deviation vector can also be obtained by explicitly solving the geodesic deviation equations.
The Jacobi stability, and its theoretical foundation, the KCC theory, offers an alternative approach to the "classical" Lyapunov approach, by investigating the deviations of the entire trajectory of the cosmological evolution equations with respect to the nearby ones under the effects of a small perturbation. In the framework of general relativity we may call the applications of the KCC theory to the study of the gravitational fields as a "second geometrization", in which already geometric quantities are supplemented by additional geometric structures. Hence general relativistic cosmological models can be described in geometric terms originating from their dynamical system structure, with these new geometric structure fully determined by the underlying Riemannian geometry, and the physical properties of the scalar fields (their self-interaction potential). The stability properties of the perturbations of a given trajectory describing the cosmological evolution are determined by the properties of the curvature deviation tensor, a geometric quantity constructed from the connections (non-linear and Berwald) associated to the dynamical system describing the cosmological evolution. It is important to note that the KCC theory can be directly applied to systems of second order differential equations, which can be interpreted geometrically as the paths (or geodesics) associated to a semispray. In investigating the Jacobi stability of cosmological models we have followed two approaches. Since the cosmological evolution equations (the Friedmann equations) are second order differential equations, the KCC theory can be naturally and directly applied to study the stability of the cosmic evolution. As a first step one obtains the two non-zero components of the non-linear connection, with the N 1 2 component depending on the product of the scale factor and of the time derivative of the field, while the N 2 2 component depends on the energy density of the scalar field, as well as of the scalar field potential. After obtaining the components of the deviation curvature tensor we have formulated the general condition of the stability of the scalar field cosmological models, which is determined by two inequalities involving the second and the first derivative of the scalar field potential, the energy density of the field, and the time variation of the scalar field itself.
As an applications of the developed formalism we have considered two scalar field models, both being relevant for the study of both the early and the late stages of the cosmological evolution. The first case we did consider is the scalar field with exponential potential. We have studied in detail the KCC geometric properties of this model. It turns out that the Jacobi stability condition, which can be expressed in terms of the components of the deviation curvature tensor is not satisfied during the cosmological evolution, and that the Universe described by the exponential potential scalar field is in a Jacobi unstable state. This result is independent on the numerical values of the parameter λ, describing the properties of the potential, and can also be inferred from the behavior of the components of the deviation vector ξ i , with ξ 1 diverging exponentially in time. As a second case we have considered the case of the Higgs type potential. For this potential the KCC geometric quantities show a complex behavior. After a period in which the scalar field and the potential are almost constant, the field starts to oscillate, with the amplitude of the oscillations decreasing in time. This behavior of the Higgs field is also reflected in the behavior of the components of the deviation curvature tensor, which are also some oscillating functions. The Jacobi stability of this cosmological model strongly depends on the numerical value of the parameter η = λ/4M 2 . For small values of η, the Universe evolves between successive Jacobi stable and unstable states. With the increase of the numerical value of η the time intervals in which the Universe is Jacobi stable decrease quickly, and for large values of η the Universe is in a Jacobi unstable state dur-ing its entire cosmological evolution.
As a second approach for the study of the Jacobi stability of scalar field cosmologies we have considered the first order dynamical system formulation of scalar field evolution equations. In this approach by introducing a new set of variables, expressed in terms of the square root of the potential, the time derivative of the scalar field, and the Hubble function, respectively, the Friedmann equations in the presence of scalar fields can be reformulated as a first order dynamical system, consisting of three highly non-linear ordinary differential equations. In order to apply the KCC theory this dynamical system must be lifted to the tangent bundle, and formulated as a second order differential system. We have analyzed, by using this approach, two specific scalar field models, the phantom quintessence and the tachyon scalar field with power law potentials. It turns out that for this choice of the potential both scalar field models are Jacobi unstable. The power law potential gives a very simple form for the function Γ(z), which takes constant values during the cosmological evolution. This situation is similar to the case of the exponential potential, and leads to a significant simplification of the mathematical formalism.
We have started our study of the applications of the KCC theory to cosmological problems with the investigation of the standard mater dominated cosmological models in their dynamical system formulation. We have studied the Jacobi stability of the critical points of different models, and we have shown that they are Jacobi unstable. A full comparison between the Jacobi and Lyapunov properties of the critical points for second order systems was given in [38] and [39], and hence we will not discuss in detail this relation. However, this study of the critical points of matter dominated cosmological models also shows the fundamental differences between the Lyapunov stability and KCC theories: while Lypunov stability is mostly restricted to the study of critical points, the KCC theory has the potential of investigating the deviation of the full trajectory during the entire period of the cosmological evolution. Therefore we can consider a Lyapunov stability analysis of steady states (called the linear analysis), and a "Lyapunov type" stability analysis of the whole trajectory (the KCC or Jacobi stability analysis), and these two methods are complementary but distinct to each other.
The KCC theory also introduces the first set of KCC invariants ǫ i , i = 1, ..., n, giving the contravariant KCC derivative of the vector field y i . The first KCC invariant can be interpreted as an external force. We did not study in detail the time evolution of the first KCC invariant, since its properties are not directly related to the stability issues that were our main points of interest.
In the present paper we have performed a stability analysis of the scalar field cosmological models, in which we have considered a description of the deviations of the whole trajectories of the differential system describing the cosmological dynamics, and we have provided some basic theoretical and computational tools for this study.
Further investigations of the Jacobi stability properties of cosmological models may provide some methods for discriminating between different evolutionary scenarios, as well as for the better understanding of some other fundamental processes, like, for example, structure for-mation, that played an essential role in the evolution of our Universe.
Conflict of interest The authors declare that there is no conflict of interest regarding the publication of this paper.
The geodesic deviation equations are obtained in the form