The past and future dynamics of quintom dark energy models

We study the phase space of the quintom cosmologies for a class of exponential potentials. We combine normal forms expansions and the center manifold theory in order to describe the dynamics near equilibrium sets. Furthermore, we construct the unstable and center manifold of the massless scalar field cosmology motivated by the numerical results given in Lazkoz and Leon (Phys Lett B 638:303. arXiv:astro-ph/0602590, 2006). We study the role of the curvature on the dynamics. Several monotonic functions are defined on relevant invariant sets for the quintom cosmology. Finally, conservation laws of the cosmological field equations and algebraic solutions are determined by using the symmetry analysis and the singularity analysis.


Introduction
Recent observations coming from type I-a Supernovae (SNIa), large scale structure (LSS) formation, and the cosmic microwave background (CMBR) Radiation [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16], indicates the universe is expanding with an accelerated rate without been settled until now a definitive model to account for this phenomena. One of the current alternatives models for such acceleration of the expansion is the yet unknown form of matter named dark energy. Several models of dark energy use a single scalar field, e.g., the so called phantom field and the quintessence field [17,18]. Furthermore, one can consider models with two fields, in particular combining quintessence and phantom fields one can built up the so-called quintom models . These models admits the crossing of the boundary w = −1 with w being the equation of state parama e-mail: genly.leon@ucn.cl b e-mail: anpaliat@phys.uoa.gr c e-mail: jorge.morales@cimat.mx eter of the cosmic fluid. This behavior has been investigated in several cosmological theories, for example: h-essence cosmologies [28,29]; holographic dark energy [46][47][48][49]; inspired by string theory [33]; by considering spinor matter [37,38]; for arbitrary potentials [40][41][42][43][44]; using isomorphic models consisting of coupled oscillators [50], Pais-Uhlenbeck oscillator [51]; inspired in scalar tensor theories [52][53][54][55][56] as well as in modified theories of gravity [57,58]. Furthermore, interacting quintom in nonminimal coupling has been studied in [59], and quintom with nonminimal derivative coupling was studied in [60]. An extension of the quintom scenario, in which the scalar fields are coupled through a kinetic interaction, have been studied in [61,62] while an interacting quintom model was recently constructed by the application of the generalized uncertainty principle (GUP) in the scalar field Lagrangian [63]. In particular, the second scalar field of the quintom model in [63] corresponds to higher-order derivatives given by the GUP. Furthermore, one of the particularly exciting solutions of quintom type models is that they can provide a classically stable bouncing solution. This was first pointed out in [61] and two important extensions were made in [62] and [63], which discussed the Lee-Wick and Horndeski realizations, respectively. The cosmology of the quintom model with exponential potentials was investigated in [30] and [31,32]. More specifically, [31,32] the potential accounts for the interaction between the conventional scalar field and the phantom field. In contrary with [30] where no such interaction was assumed. In [31] was found that the interaction term dominates over the mixed terms of the potential. This simplification has important consequences such that the existence of scaling attractors excludes the existence of phantom attractors, as a difference with the results in [32]. Some of these results were extended in [40], for arbitrary potentials. In the recent years still there an special interest in this kind of models with the phantom crossing behavior [67][68][69][70][71][72][73]. In line with this interest, we stress in this paper a deep insight about the model proposed in [31].
The plan of the paper it follows. In Sect. 2 we make a general description of the model under study and we present the field equations. In Sect. 2.1 are presented several monotonic functions defined on invariant sets of the phase space using the results summarized in the Appendix A. The equilibrium points of the model for any curvature choice using H -normalization are discussed in Sect. 2.2. In Sect. 2.2.1 we make the analysis of zero-curvature models with two scalar fields, one a linear function of the other. Furthermore, in Sect. 2.3 we present some improvements of the results from reference [31]. In particular, in Sect. 2.3.1 we discuss the normal expansion with undetermined coefficients of arbitrary order that is used to construct the unstable and center manifold up to fourth order in [74] for flat FRW metric. We explicitly state in the Appendix B the main techniques to use, following the approach of [75][76][77]. In Sect. 2.4 we introduce an alternative of the variables used in [31] for flat FRW models than renders the phase space in compact form. Moreover, in Sect. 3 we propose several normalization procedures for investigating models with flat and non flat spatial curvature, that allows to circumvent some technical problems that arises, for example when the crossing through H = 0, where H is the Hubble factor, happens. The integrability of the quintom cosmological model is studied in Sect. 4. In particular we apply two different methods to determine the integrability, the symmetry analysis and the singularity analysis. Finally, in Sect. 5 we discuss our results and draw our conclusions.

The quintom model
In this work we consider the so-called two-field quintom with Lagrangian [31]: where the scalar field φ represents quintessence and ψ represents a phantom field, and we include a co-moving perfect fluid in the gravitational action. We use the Friedmann-Robertson-Walker (FRW) line element given in spherical coordinates by: ds 2 = −dt 2 + a(t) 2 dr 2 1 − kr 2 + r 2 dθ 2 + r 2 sin 2 θ dϑ 2 . ( The curvature of spatial slices is given by: where k identifies the three types of FRW universes: closed, flat, and open, respectively. The field equations derived for the line element (3) in the context of General Relativity, are given by where H =ȧ (t) a(t) denotes the Hubble function, while the dot denotes differentiation with respect to variable t. Moreover, for the barotropic index of the background matter, γ, we assume that it takes values only in the interval 0 < γ < 2.
We introduce the following normalized variables the Friedmann equation becomes For convenience we shall write q ≡ −äa/ȧ 2 (the deceleration factor) in terms of the new variables as The new variables (63) can be used to rewrite the system (5) as a system of first-order ordinary differential equations (ODEs): The primes denote derivative with respect to τ = log a 3 .

Invariant set and monotonic functions
By following the procedures in the Appendix A, we can find the invariant sets.

By using the equations
, it follows from Proposition 10 that the sign of 3 R and the sign of nx φ + mx ψ are invariant. From equations (9)- (12) it is obvious that any combination of signs for y and z is an invariant set for the flow.
Defining on this sets the function Z = x 2 φ − x 2 ψ + y 2 − 1 that has directional derivative along the flow given by Z = α Z where α is a continuous real-valued non-vanishing function on R 3 given by α = x 2 φ − x 2 ψ − y 2 . Hence, by applying Proposition 10, the desired result follows. In the same way it is proved that y = 0 is an invariant set of the flow by defining α = 1 + q − 3 mx φ + nx ψ . Combining the former results we can prove that the 2-dimensional set Notice that the equilibrium points O and C ± studied in [31] are in the boundary of .
Thus, W 1 is positive and strictly monotone along orbits of the invariant set S 1 . The boundary of S 1 is the invariant set Thus, W 2 is positive and strictly monotone along orbits of the invariant set S 2 . Hence, the boundary of S 2 is the invariant setS 2 \S 2 = {x : nx φ + mx ψ = 0} ∪ {x : z = 0}. Let be any initial state W 3 is positive monotonic decreasing (increasing) in the invariant set S 3 for non-zero curvature models provided γ > Then, W 4 is positive and strictly monotone along orbits of the invariant set S 4 .S 4 \S 4 = {x : Let be any initial state x 0 ∈ S 4 such that O(x 0 ) remains bounded. By Theorem 1 we find α(x 0 ) ⊂ C := {x :

Remark 3
In [31] it was devised the monotonic function , defined in the invariant set of zero-curvature models with dust matter. There it was argued that the orbits initially near the hyperbolas (representing massless scalar field cosmologies) tends asymptotically to the flat FRW universe. On the other hand, from the monotonic function W 4 we can conclude that (for zero-curvature models) the orbits which initially near the hyperbolas tends asymptotically to either the matter scaling (MS) solutions (equilibrium points like T ) or (in the absence of scaling solutions) to powerlaw solutions (equilibrium points like P).

Remark 4
The previous gallery of monotonic functions can be used to ruled out, according to Proposition 11, the existence of homoclinic orbits, recurrent orbits and more complex behavior in several invariant sets. The dynamics is dominated by equilibrium points and heteroclinic orbits joining it. Table 1 Location, existence and some properties of the equilibrium points for z ≥ 0, and y ≥ 0. We use the notation δ = m 2 − n 2

The equilibrium points
The system (9)-(12) admits six types of equilibrium points (with both z and y non-negative) denoted by M, F, C ± , P, M S, and C S. In Tables 1 and 2 we offer some partial information of these equilibrium points like the location, conditions for existence, dynamical character and some additional information. For zero-curvature models (with dust background, i.e., γ = 1), the equilibrium point M S and P reduces respectively to the points T and P, investigated in [31]. The equilibrium point M represents the Milne universe, it has negative curvature. It is no hyperbolic if γ = 2 3 . F denotes the flat FRW universe and it is no hyperbolic if γ = 2 3 or γ = 2. For γ < 2/3 (resp. γ > 2/3) the eigenvalues of M (resp. F) are (−, −, +, +) and the eigenvalues of F (resp. M) are (−, −, +, −). In this case both equilibrium points act locally as saddles. The 3-dimensional stable subspace of F (resp. M) intersects the 2-dimensional unstable subspace of M (resp. F) at the z axis. This means that the orbits initially at this axis are asymptotic to the past to the Milne (resp. flat FRW) universe and to the flat FRW (resp. Milne) universe toward the future.
The equilibrium sets given by the hyperbolas C ± were investigated in [31], the same analysis done in that reference is valid in our more general context. In the former section we conjecture that they are local sources for the quintom phasespace. In fact, if δ > 1 then part of the hyperbola is a local attractor to the past, and part of it is a saddle. If δ < 1 then the entire hyperbola is a local attractor to the past.
As stated in [31], point P is the local attractor for zerocurvature models provided there are not scaling regimes. As we have seen, for some region in the parameters space the equilibrium points P and M S exist in the same phase portrait with P having eigenvalues (−, −, −, +) (i.e, a saddle) and M S being a local attractor (i.e., all its eigenvalues having positive real part). If 0 < γ ≤ 2/3, 1/3 < δ < 1 or 2/3 < γ < 2, γ/2 < δ < 1 the eigenvalues of P are (+, −, −, +). In the latter case, the stable subspace of P is 2 dimensional and then it is the local attractor for zerocurvature models without matter.

A flat FRW model with one scalar field a linear function of the other
The dynamics is governed by the equations: defined on the phase plane: The conditions for existence and the dynamical behavior of the equilibrium points are presented in Table 3.
In Fig. 1 it is shown some orbits of the system (13) . Fig. 1 Some orbits in the phase plane of the system (13)- (14) for five choices of the parameters represents the physical portion of the phase plane which is either half ellipse or half hyperbola.

Flat FRW geometry and dust matter
For the flat FRW geometry and dust matter dynamics is governed by the vector field [31]: defined in the phase space given by Here the prime denotes differentiation with respect to a new time variable τ = log a 3 , where a is the scale factor of the space-time. The deceleration factor q ≡ −äa/ȧ 2 can be written as In Table 4 is is presented some information concerning the equilibrium points O, C ± , T, P reported in [31]. The stability of the system is as follows [31]: Case (i) If m < n 2 + 1/2, the point P is an stable node and T does not exists. Case (ii) For n 2 + 1/2 < m ≤ n 2 + 4/7, the point T is a stable node and P is a saddle point. Case (iii) For n 2 + 4/7 < m < √ 1 + n 2 , point T is a spiral node and the point P is a saddle. Case (iv) For m > √ 1 + n 2 the point T is a spiral node whereas the point P does not exist.
For 3 R = 0, z = 0, that is, the quintom cosmological model (1) with potential (2) in the spatially flat geometry and without any other matter source, the dynamical system is reduced to Table 4 Location, existence and deceleration factor of the equilibrium points for m > 0, n > 0 and y > 0. We use the notation δ = m 2 − n 2 (from Ref. [31]) defined on the phase plane In the Fig. 2 are presented some orbits in the phase plane of the system (20)- (21) for three choices of the parameters. This figures depict the typical behavior of the quintom model with exponential potential and no other matter component (vacuum quintom), where the P is always an stable node. The shadowed region corresponds to the physical region. The region w < −1 is reached provided m 2 − n 2 < 0, so the crossing does indeed occur in the case m = 0.5, n = 0.6 whence P is a phantom-like solution. For m = 0.6, n = 0.5 the equation of state remains above w = −1 and P is a quintessence-like solution, and for m = n = 0.2, w = −1, such that P represents a de Sitter solution.

Normal expansion up to arbitrary order
In this section we obtain the normal form expansion for an arbitrary point lying on the curve C − [74]. (15)

Proposition 5 Let be the vector field X given by
∈ Z, then, there exist constants a r , b r , c r , r ≥ 2 (non necessarily different from zero) and a transformation of coordinates x → y, such that (15)-(17) has normal form defined in neighborhood of y = (0, 0, 0).
Proof Applying a linear coordinate transformation the linear part of X is reduced to Let L J be the linear operator that assigns to h(y) ∈ H r the Lie bracket of the vector fields Jy and h(y): Applying this operator to monomials x m e i , where m is a multiindex of order r and e i basis vector of R 3 , we find The eigenvectors in B for which m,i ≡ (m · λ) − λ i = 0 form a basis of B r = L J (H r ) and those such that m,i = 0, associated to the resonant eigenvalues, form a basis for the complementary subspace, G r , of B r in H r .
Since λ 1 = 0, and λ 3 / ∈ Z by hypothesis, the resonant equations of order r has a unique solution: Then, By applying Theorem 2, we have that there exists a coordinate transformation x → y, such that (15)-(17) has normal form (23)- (25) where a r , b r and c r are some real constants. The values of of all these constants can be uniquely determined, up to the desired order, by inductive construction in r.

If
Then, the orbit approach the origin as τ → −∞ provided λ 3 > 0. 2. If y 10 = 0, we obtain the quadratures The y 1 -component of the orbit passing through (y 10 , y 20 , y 30 ) at τ = τ 0 with y 10 = 0 is obtained by inverting the quadrature (30). The other components are given by

Normal expansion to third order for C −
In this section we show normal form expansions for the vector field (15)-(17) defined in a vicinity of C − expressed in the form of Proposition 6 (see [74]) Proposition 6 Let be the vector field X given by (15) ∈ Z, then, there exist a transformation to new coordinates x → z, such that (15)-(17), defined in a vicinity of x * , has normal form where

Fist step: transforming to the Jordan form
Using a linear coordinate transformation the system (15)-(17) is transformed to where x stands for the phase vector with 2 . Second step: simplifying the quadratic part By the hypotheses the eigenvalues To obtain the normal form of (35)-(37) wee must look for resonant terms, i.e., those terms of the form x m e i with m and i such that m,i = 0 for the available m, i. Only one term of second order is resonant : The required function to eliminate the non-resonant quadratic terms is given by The quadratic transformation with h 2 defined as in (42) Since − L (2) J h 2 (y) + X 2 (y) = X (1,0,1),3 y 1 y 3 e 3 , we have y = Jy + X (1,0,1),3 y 1 y 3 e 3 +X 3 (y) + O(|y| 4 i.e., c 2 = X (1,0,1 The vector fieldX 3 (y) introduced above has the coefficients: Third step: simplifying the cubic part After the last two steps, the Eq. (38) is transformed to As we will see later on, there is only one term of order three which is resonant: (2,0,1),3 = 0 → c 3 z 2 1 z 3 e 3 . Therefore, as in the last step, in order to eliminate non-resonant terms of third order we will consider the coordinate transformation y → z given by where where (1,2,0 (49) is the required by Theorem (2). By using this theorem we prove the existence of the required constant c 3 . To find its we must to calculate where e i , i = 1, 2, 3, is the canonical basis in R n . Then, Observe that the transformation h 3 does not affect the value of the coefficient of the resonant term of order r = 2. Then, the result of the proposition follows.
If λ − 3 > 0, the unstable manifold of the origin is, up where δ is a real value small enough. Therefore, the dynamics of (35)-(37) restricted to the unstable manifold, is given, up to order O(|z| 4 0, 0, 0) . Then, the origin is the past attractor for an open set of orbits of (35)- (37).
For λ − 3 > 0, the center manifold of the origin is, up to order O(|z| 4 ), W c loc (0) = {(z 1 , z 2 , z 3 ) ∈ R 3 : z 2 = z 3 = 0, |z 1 | < δ} where δ is a real value small enough. Or as a graph in the original variables defined by By taking the inverse, up to fourth order, of (52) we have the expression for z 1 : Substituting (54) in (51)-(53) we have that the center manifold of the origin is given by the graph:

Compact variables formulation
The Friedman equation for k = 0 can be written as It is clear that when H 2 + 1 6ψ 2 + 2 3 V (φ, ψ) = 0 we obtain for nonnegative potential the trivial solution. Since this solution is not relevant in our set up we consider H 2 + 1 6ψ whereD They satisfy That is, we have five variables and two constraints. Therefore, we can study the reduced system defined in the compact phase space where the prime means derivative and ε = ±1 appears due to (57) cannot be solved globally for , and it corresponds to the sign of H , where = +1 corresponds to ever expanding universe and = −1 corresponds to ever contracting universes. The equilibrium points/lines have coordinates (X φ , X ψ , Y ): 1. O ± := (0, 0, 0). For = +1 we obtain a representation of the point O referred in Table 4. 2. C ± := (±1, X ψ , 0) referred in Table 4. Table 4 is recovered by setting = +1.  Table 4 is recovered by setting = +1.
In Fig. 3a-d are presented some orbits in the phase plane of the system (58), (59), (60) for the four choices of the parameters chosen in [31] and ε = +1.
The static universe configurations cannot by discussed in the framework of [31]. This illustrates the applicability of this formulation in compact variables. This case is shown in Fig 3e.

Alternative normalization
Due the negative sign of the kinetic energy of the phantom field, it is possible that H can be zero in a finite time. If so, the variables in [31] are not well defined as H → 0. They can diverge in a finite parameter time, whence the resulting system is not a dynamical system. To avoid this difficulty we will use a set of dynamical variables similar to those introduced in [82] section III.E, instead of the ones used in [31], for investigating negative (and zero curvature models). For investigating positive curvature models we shall make use of similar variables to those defined in [82] section VI.A.

Negative and zero curvature models
In this section we extent the analysis presented by one of us in [83], since now it is included the possibility of the crossing through H = 0.
Let us introduce the following set of normalized variables: This choice of phase-space variables renders the Friedmann equation as By definition −1 ≤ D ≤ 1. Then, From the requirement of the non-negativeness of the energy densities of quintom and matter, follows de , k denotes the fractional energy densities of dark matter, dark energy and "curvature". Then, from (65) follows that is bounded as de is. However, we cannot guaranteed that the variables U φ , U ψ , W are bounded. In this sense, the dynamical system constructed from the variables (63) is unbounded (then, non-compact). Hence, we cannot guaranteed that the system admits both future and past attractors. Numerical experience (see for instance [31]) supports the conjecture that past and future attractors can exists depict unboundedness.
For D = 0 (i.e., H = 0) the variables (63) are related with those in [31] via The equation of state (EoS) parameter, w, and the deceleration parameter, q, can be stated in terms of our variables as where The evolution equations for (63) are where we have defined the new independent variable The ODEs (69) define a flow on the state space Observe that the variable D does not decouples from the other variables. Hence, the system is of higher dimension (5-dimensional) that the system we can obtain by using the variables defined in [31] (4-dimensional in the case of negative curvature).

Remark 6
In the case D = 0 we find from (65) and from the requirement of non-negativeness of the energy densities of the sources, that = 0 and U 2 φ − U 2 ψ + W 2 = 0. Then, Equilibrium points of quintom model with k = 0, −1. We use the notations δ = m 2 − n 2 and λ ± = [×s] means multiplicity s. The subscripts on the label have the following significance. The left subscript gives the sign of D (denoted by = ±1) and indicates whether the corresponding model is expanding (+) or contracting (−). The right subscript gives the sign of U φ (i.e., the sign ofφ) and is denoted by the sign ±. When the flow is restricted to the invariant sets D = ±1, the eigenvalues associated to the equilibrium points ± M, ± F, ± S F and ± M S an to the equilibrium sets ±K± , are, in each case, the same as those displayed, but the first from the left Label Coordinates Existence Eigenvalues From the last inequality we can conclude that the submanifold D = 0 is not invariant for models with W > 0, but acts as a membrane . For massless scalar field (MSF) cosmologies, the submanifold D = 0 is invariant and is equivalent to (69) is invariant under the transformation of coordinates Thus, it is sufficient to discuss the behavior in one part of the phase space. The dynamics in the other part being obtained via the transformation (72).

Remark 7
Let be defined in the interior of the phase space the function This function is a monotonic function in the regions 0 < D < 1 and −1 < D < 0 for > 0 and n U φ + m U ψ = 0. From the expression of M we can see immediately that either D 2 → 1 or → 0, or n U φ + m U ψ → 0 or n U φ + m U ψ → +∞ (which means that U φ or U ψ or both diverge) or k → 0 asymptotically. Particularly, M is a monotonic decreasing function for models with 0 < D < 1 (i.e., for ever expanding models) and takes values in the interval (0, +∞). Since M → 0 as → 0 or as n U φ + m U ψ → 0 (since its denominator is always bounded), then, the future asymptotic state of the system for ever expanding models corresponds to models with = 0 or to models where n U φ + m U ψ = 0. Since M → +∞ as n U φ + m U ψ → +∞ or as D → 1 or as k → 0, then the past asymptotic state of the system for ever expanding models corresponds to models where at least one of the variables U φ or U ψ diverges, or to models with H large, or to flat models.
The asymptotic dynamics for contracting models can determined from the asymptotic dynamics of ever expanding models by means of the coordinate transformation (72).

Remark 8
The existence of the monotonic function (73) in our state space rules out any periodic, recurrent, or homoclinic orbit in the interior of the phase space. Therefore, some global results can be found from the linear (local) analysis of equilibrium points.
In the Table 5 it is displayed the conditions for existence and the eigenvalues of the linearized system around each equilibrium point of the system (69). In Table 6 we give the cosmological parameters associated to the equilibrium points. In the following we will characterize them.
The sets of equilibrium points V 1,2 and ± K ± and the isolated equilibrium points ± M are located in the invariant set of MSF cosmologies without matter. The isolated equilibrium point ± F is located in the invariant set of MSF cosmologies with matter.
The arcs ± K ± denotes cosmological models dominated by dark energy ( de → 1) particularly by its kinetic energy. The DE mimics a stiff fluid solution. The stable (resp., unstable) manifold of − K ± (resp., + K ± ) is at least 3-dimensional in the full phase space.
The equilibrium points ± M denotes the Milne's universe. They are non hyperbolic for γ = 2 3 . The equilibrium points ± F represent the flat FRW universe. They are non-hyperbolic if γ = 2 3 or γ = 2. For γ > 2 3 , both ± M and ± P acts locally as saddles. The unstable (resp., stable) manifold of − M (resp. + M) is 3-dimensional in the full phase space. The equilibrium point + M (resp. − M) is unstable (resp. stable) to perturbations in D and in W. The stable (resp., stable) manifold of − F (resp. + F) is 3-dimensional. Whereas, the unstable (resp. stable) manifold − F (resp. + F) is 2-dimensional in the full phase space. If we restrict the dynamics to the invariant sets with D 2 = 1 the dimensionality of all the above invariant sets reduces in one unity.
The isolated equilibrium points ± S F and ± C S denotes scalar field dominated solutions and curvature scaling solutions, respectively. They are located in the invariant set W > 0, = 0. The equilibrium points ± M S (belonging to the invariant set W > 0, > 0) represents flat matter scaling solutions. Observe that the equilibrium points ± M S and ± S F coincide as δ → γ 2 + and ± C S and ± S F coincide as δ → 1 3 + . Additionally, + S F (resp., − S F) coincides with a point in the arc + K + (resp., − K − ) as δ → 1 − . These values of δ where equilibrium points coincide corresponds to bifurcations.
Due the symmetry (72) it is sufficient to characterize the behavior of the equilibrium points + S F, + C S, and + M S. The dynamical character of their counterparts in the "negative" branch − S F, − C S, and − M S is determined by the symmetry (72).
The equilibrium point + S F is non hyperbolic in the full phase space if either δ = 1 3 or δ = γ 2 or δ = 0. However, if we restrict our attention to the flow in the invariant set D = 1, the eigenvalue 2δ plays no role in the dynamics (it appears due the "extra" dimension, D). In fact, if δ = 0, the eigenvalues of the system restricted to D = 1 are (−2/3, −1[×2], −γ ). Then, + S F represents a de Sitter solution, and it is the tractor in the invariant set D = 1. If δ < 0 and 2 3 < γ < 2 then, + S F is the future attractor (in the full phase space) for (ever) expanding models and it represents a phantom field. Using the symmetry (72) we find that, if δ < 0 and 2 3 < γ < 2 then, − S F is the past attractor for contracting models and it represents a phantom field. If 0 < δ < 1 3 and 2 3 < γ < 2, then + S F represents a solution dominated by a quintessence field, and it is the future attractor for the flow restricted to the invariant set D = 1. If 0 < δ < 1 3 and 2 3 < γ < 2, then − S F represents a solution dominated by a quintessence field, and it is the past attractor for the flow restricted to the invariant set D = −1.
The equilibrium point + C S represents a solution with negative curvature. It is non hyperbolic if δ = 1 3 . It is a saddle point for the full phase space, since it have always at least a positive eigenvalue. It is never a source (all the eigenvalues are not simultaneously with positive real parts) for the class of ever expanding models. By using the symmetry argument we find that − C S is never a global attractor (or sink) for contracting models. But, if we restrict our attention to the flow in the invariant set D = 1 we find that in this case + C S is the future attractor provided δ > 1 3 and 2 3 < γ < 2. In other words, one can say that if 2 3 < γ < 2 the spatial curvature destabilizes + S F at the bifurcation value δ = 1 3 and that the stability is transferred to the equilibrium point + C S (of course, this argument is valid only in the invariant set D = 1).
The equilibrium point + M S represents a cosmological solution where DE tracks DM (it means that the DE EoS parameter scales as the DM EoS parameter). It is non hyperbolic for γ = 2/3 or γ = 2. It is in general a saddle point. However, when the dynamics is restricted to the invariant set D = 1 it can be the future attractor (for ever expanding models) provided 0 < γ < 2 3 and δ > γ 2 . From our hypothesis, this is not a region of physical interest. However, if we do not restrict the value of γ to be > 2 3 the above fact can be interpreted as follows. If 0 < γ < 2 3 the matter destabilizes + S F at the bifurcation value δ = γ 2 and that the stability is transferred to the equilibrium point + M S. Otherwise, whenever 0 < γ < 2 3 and δ < γ 2 the late time attractor in the invariant set D = 1 is the equilibrium point + S F. These equilibrium points are such that n U φ + m U ψ = 0. Hence, by the structure of the function (73), for ever expanding models (i.e., for 0 < D < 1), the late time behavior is determined by the local behavior of either + S F or + C S or + M S. By making use of the symmetry (72) it is possible to determine the early time behavior of contracting models (i.e, for −1 < D < 0) from the local behavior of either − S F or − C S or − M S. Combining the above linear analysis, with the information we have using monotonic functions we have the following (a) The past attractors are: 3. For 0 < D < 1 (a) The past attractor is + K ± if nU ψ ±m 1 + U ψ 2 < 1.

Positive curvature models (k = +1)
For investigating positive curvature models we shall make use of variables similar to those defined in [79] section VI.A. We reproduce the results found by one of us in [80]. Let us introduce the normalization factor Observe that (i.e., at a singularity). This means that it is not possible that D vanishes at a finite time. Let us introduce the following normalized variables (Q 0 ,Û φ ,Û ψ ,Ŵ ,ˆ ), given by From the Friedmann equation we find and by definition By the restrictions (77), (78), the state variables are in the state spacê As before, this state space is not compact. Let us introduce the time coordinate D has the evolution equation This equation decouples from the other evolution equations. Thus, a reduced set of evolution equations is obtained.
There is also an auxiliary evolution equation Table 7 equilibrium points in the quintom model with k = +1. We use the same notation as in Table 5. If γ = 2 3 , the set of equilibrium points, S, (corresponding to static solutions) reduces to two sets of equilibrium points analogous to V 1,2 (see Table 5). When the flow is restricted to the invariant sets Q 0 = ±1, the eigenvalues associated to the equilibrium points ±F , ±Ŝ F and ±M S an to the equilibrium sets ±K± , are, in each case, the same as those displayed, but the first from the left

Label
Coordinates: Existence Eigenvalues It is useful to express some cosmological parameters in terms of our state variables.
Observe that the system (80), (81) is invariant under the transformation of coordinates Thus, it is sufficient to discuss the behavior in one part of the phase space. The dynamics in the other part being obtained via the transformation (82).

Remark 9 The function
is monotonic in the regions Q 0 < 0 and Q 0 > 0 for Q 2 0 = 1, nÛ φ +nÛ ψ = 0,ˆ > 0. Hence, there can be no periodic orbits or recurrent orbits in the interior of the phase space. Furthermore, it is possible to obtain global results. From the expression N we can immediately see that asymptotically In the Table 7 it is summarize the location, existence conditions and the eigenvalues of the linearized system around each equilibrium point. In the following we will characterize the dynamical behavior of the cosmological solutions associated with them.
The equilibrium points ±K , ±F , ±Ŝ F and ±M S represents flat FRW solutions. They corresponds to the same cosmological solutions as the not hatted ones. We submit the reader to the early sections for the physical interpretation of them. Although ±Ĉ S represent a curvature scaling solution, it is different from ± C S in the sense that it represents a positive curvature model.
By the above discussion, we will focus here on the dynamical character of the equilibrium points, not in its physical interpretation. As before, due the symmetry (82), we will characterize the equilibrium points corresponding to the "positive" branch. The dynamical behavior of the equilibrium points in the "negative" branch is determined by the transformation (82).
Additionally, +Ŝ F (resp., −Ŝ F) coincides with a point in the arc +K+ (resp., −K− ) as δ → 1 − . These values of δ where equilibrium points coincide corresponds to bifurcations. Combining the above linear analysis, with the information we have using monotonic functions we have the following The past attractor are: The future attractors are:

Integrability and algebraic solution
In this section we examine the integrability for the gravitational field equations of the quintom model by using two methods for the study of the integrability of dynamical systems. In particular we apply the symmetry analysis, we determine point transformations which leaves invariants the differential equations which provide also conservation laws. The second method that we assume is the singularity analysis where we are able to write the algebraic solution for the cosmological quintom model of our consideration.

Symmetry analysis
The symmetry analysis has an important role for the study of the integrability and the determination of analytical solutions in gravitational theories [84][85][86]. There is a plethora of applications in the literature of the symmetry analysis in cosmological studies, in the dark energy models and in the modified theories of gravity, for instance see [87][88][89][90][91] and references therein. Consider now the point-like Lagrangian for the quintom cosmological model (1) with potential (2) in the spatially flat geometry and without any other matter source The field equations are the Euler-Lagrange equations of the later point-like Lagrangian plus the Hamiltonian constraint which is the first Friedmann equation. By following the Noether symmetry analysis, for a recent review on the approach and more details for the method we refer the reader to [90], we find that the Lagrangian (84) admits the following two Noether symmetries, except the trivial autonomous symmetry, where the corresponding Noetherian conservation laws The three conservation laws, H, 1 and 2 are linear independent, and in-involution, that is { 1 , 2 } = 0, { 1 , H} = 0 and { 2 , H} = 0, where {, } denotes the Poisson bracket. Thus, we can refer that the gravitational field equations described by the Lagrangian (85) form a Liouville integrable dynamical system [92].
Notice that the above conservation laws can be written in terms of the phase plane variables as such that at the invariant set {x ∈ R 2 : nx φ + mx ψ = 0} we have the trivial charge 1 = 0, whereas in the invariant set {x ∈ R 2 : x φ −x ψ m+n = 1} we have the trivial charge 2 = 0. Because of the nonlinearity of the differential equations it is not possible to write the analytical solution in closedform by using well-known functions. Hence, we continue our analysis with the singularity analysis where the solution can be written by using Laurent expansions.

Singularity analysis
Except from the symmetry analysis, singularity analysis is another approach to study the integrability of a given dynamical system. In singularity analysis integrability is not concerned with the display of explicit functions but with the demonstration of a specific property, such that: the existence of a Laurent series for each of the dependent variables. The series may not be summable to an explicit form, but does represent an analytic function apart from any singular points.
Moreover, the essential feature of this Laurent series is that it is an expansion about a particular type of movable critical point, a pole.
While the symmetry analysis is related with the existence of closed-form solutions or solutions expressed in terms of well-known first-order differential equations. There is a choice of the type of definitions for the symmetries which can be applied in the symmetry analysis, where different conclusions can be occurs. On the other hand, singularity analysis is straightforward in its application as it does not offer so many choices. Both methods are complementary [93].
The main steps for the application of the singularity analysis for a given differential equation are described by the ARS algorithm (from the initial of Ablowitz, Ramani and Segur) [94][95][96]. The latter algorithm has been used in various cosmological models to determine the integrability and write algebraic solutions, for instance see [97][98][99][100][101][102] and references therein. For a demonstration on the application of the ARS algorithm we refer the reader in the review of Ramani et al. [103]. Consider the dynamical system (15), (16) and (17 ) which describes the evolution of the universe in a spatially flat FRW universe. We search for singular behaviour of the form where t 0 is an arbitrary constant and denotes the position of the singularity. Thus we find the unique leading-order behavior to be with constraint equation From the latter it follows that two of the coefficients of the leading-order terms are arbitrary constants. In particular are the two constants of integrations, recall that t 0 is the third integration constant. Hence, we can conclude that equations (15), (16) and (17 ) form an integrable dynamical system. However, in order to verify it we proceed with the next steps of the ARS algorithm. In order to determine the resonances, we substitute in the dynamical system 15), (16) and (17) the following expressions in which ε 1 , ε 2 , ε 3 are infinitesimal parameters. We follow the steps described in [103] for the determination of the resonances in systems of differential equations we find the three resonances to be s 1 = −1 , s 2 = 0 and s 3 = 0.
Resonance s 1 = −1 is important to exist in order the singularity analysis to succeed. In particular it is related with the movable singularity. On the other hand, the values zero of the other two resonances tell us that two of the coefficient terms for the leading-order behavior are integration constants of the problem. Moreover, because s 2 , s 3 = 0, the solution of the dynamical system is expressed by Right Painlevé series, that is, where the first coefficient constants x φ 1 , x ψ 2 and y 1 are calculated to be Thus we can conclude that the dynamical system (15), (16 ) and (17) passes the singularity test and its algebraic solution is given by the Right Painlevé Series (101), (102), (102). From the latter we can infer that the leading-order term (94), i.e. the dominant behavior close to the singularity is an unstable solution.
Proposition 9 The dynamical system described the four first-order equations (9), (10), (11) and (12), passes the singularity test with leading-order terms given by (107). The algebraic solution is given by the Right Painlevé Series where x = x φ , x ψ , y, z and constraint equation (109).
From the form of the Laurent expansions we refer that the leading-order term describe a past attractor for the dynamical system.

Concluding remarks
We discussed some results concerning the asymptotic dynamics of quintom cosmologies with exponential potential. Using in a combined way several tools of the dynamical systems theory such as the linear stability analysis, Normal Forms calculations, Center manifold and Invariant Manifold Theorems, by developing Monotonic functions and implementing numerical solutions, one can find information bout the early time and late time behavior and at the intermediate stages of the evolution. We have divided the analysis of the negative curvature and zero curvature models from the analysis of positive curvature models.
For negative or zero curvature models the physical behavior of a typical quintom cosmology as as follow. For ever expanding cosmologies, near the big-bang a typical model behaves like a flat FL model with stiff fluid (i.e. the dark energy mimics a stiff fluid) represented by the equilibrium set + K + or by + K − (parameterized by the real value U ψ ) depending if nU ψ + m 1 + U ψ 2 or nU ψ − m 1 + U ψ 2 is less than 1. If 2 3 < γ < 2 and δ < 0 the late time dynamics is determined by + S F, i.e. the model is accelerating, close to flatness ( k → 0) and dominated by phantom dark energy ( de → 1). This means, that typically, ever expanding quintom model crosses the phantom divide (it EoS parameter takes values less than −1). The intermediate dynamics is will be governed to a large extent by the fixed points + C S, + M S, and + M, which have a lower dimensional stable manifold. For H large and positive (i.e., in the invariant set D = 1), the late time dynamics is determined by the equilibrium point + S F provided δ < 1 3 which represents quintessence (−1 < q < 0, i.e. −1 < w < − 1 3 ) or a phantom field (q < −1, i.e. w < −1) if δ > 0 or δ < 0 respectively. If δ = 0 it represents a de Sitter cosmological model. Curvature dominates the late time evolution in the invariant set D = 1 whenever 2 3 < γ < 2 and δ > 1 3 . On the other hand, for contracting models, the typical behavior is, in one sense, the reverse of the above, that is, if 2 3 < γ < 2 and δ < 0 the early time dynamics is determined by − S F, i.e. the model is accelerating, close to flatness ( k → 0) and dominated by phantom dark energy ( de → 1). The intermediate dynamics is will be governed to a large extent by the fixed points − C S, − M S, and − F, which have a lower dimensional stable manifold. For H large and negative (i.e., in the invariant set D = −1), the early time dynamics is determined by the equilibrium point + S F provided δ < 1 3 . Curvature dominates the early time evolution in the invariant set D = 1 whenever 2 3 < γ < 2 and δ > 1 3 . A typical model behaves at late times like a flat FL model with stiff fluid (i.e. the dark energy mimics a stiff fluid) represented by the equilibrium set − K + or by − K − (parameterized by the real value U ψ ) depending if nU ψ + m 1 + U ψ 2 or nU ψ − m 1 + U ψ 2 is greater than −1.
The physical interpretation of the equilibrium points corresponding to positive curvature models was given in [83]. The physical behavior of a typical closed quintom cosmology is as follows [83]. For ever expanding cosmologies, near the big-bang a typical model behaves like a flat FRW model with stiff fluid represented by the equilibrium set +K+ or by +K− , depending on the choice of the values of the free parameters m, n and x ψ . If δ < 1/3 and 0 < Q 0 < 1 the late time dynamics is determined by +Ŝ F, (with the same physical properties as + S F). The intermediate dynamics will be governed to a large extent by the fixed points +Ĉ S, +M S and +F which have the highest lower-dimensional stable manifold. For flat models (i.e., in the invariant set Q 0 = 1), the late time dynamics is determined by the equilibrium point +Ŝ F provided δ < 1/2 or +M S provided δ > 1/2 .
For contracting models, the typical behavior is, in one sense, the reverse of the above. If δ < 1/3 and −1 < Q 0 < 0 the early time dynamics is determined by −Ŝ F. The intermediate dynamics will be governed to a large extent by the fixed points −Ĉ S, −M S, −F which have the highest lowerdimensional stable manifold. For flat models (i.e., in the invariant set Q 0 = −1), the early time dynamics is determined by the equilibrium point −Ŝ F (or −M F ) provided δ < 1/2 (δ > 1/2). A typical model behaves at late times like a flat FRW model with stiff fluid (i.e. the dark energy mimics a stiff fluid) represented by the equilibrium set −K+ or by −K− depending on the choice of the values of the free parameters m, n and x ψ . Furthermore, we applied the symmetry analysis and the singularity analysis to determine the integrability of the quintom model. By using the Noether point symmetries in the case of the flat FRW geometry without any matter source we found that the point-like Lagrangian admits three linearindependent conservations laws which are in involution. Therefore, the field equations form a Liouville integrable dynamical system. However, because of the nonlinearity of the field equations that analytical solution can not be written in close-form expression.
On the other hand, by using the singularity analysis we were able to prove the integrability either when matter source exists and also the spatial curvature is not zero. We wrote the analytic solution of the field equations by using Painlevé Series by starting from a unstable singular solution.
The two methods, symmetries and singularities, are complementary, while in a future analysis will be of special interest to determine new analytical solutions for other kind of potentials.
We have the following result Theorem 2 (theorem 2.3.1 in [77]) Given a smooth vector field X(x) on R n with X(0) = 0, there is a polynomial transformation to new coordinates, y, such that the differential equation x = X(x) takes the form y = Jy + (B2)