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.


I. 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 parameter of the cosmic fluid. This behavior has been investigated in several cosmological theories, for example: hessence 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 Section II we make a general description of the model under study and we present the field equations. In Section II A 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 Section II B. In Section II B 1 we make the analysis of zero-curvature models with two scalar fields, one a linear function of the other. Furthermore, in section II C we present some improvements of the results from reference [31]. In particular, in Section II C 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 section II D we introduce an alternative of the variables used in [31] for flat FRW models that renders the phase space in compact form. In this new variables we find static universe configurations that cannot by discussed in the framework of [31]. This illustrates the applicability of this formulation in compact variables. Moreover, in Section III 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 Section IV. In particular we apply two different methods to determine the integrability, the symmetry analysis and the singularity analysis. Finally, in Section V we discuss our results and draw our conclusions.

II. 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 . (3) The scalar curvature of the 3-metric induced on the space-like surfaces orthogonal to e 0 = ∂ t is given by 3 R = 6k a 2 , k = 1, 0, −1, 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, the barotropic index of the background matter, γ, is assumed to take values on the interval 0 < γ < 2.
Introducing the following normalized variables the Friedmann equation becomes and the deceleration factor q ≡ −äa/ȧ 2 can be expressed in terms of the new variables as The new variables (5) can be used to rewrite the system (4) as a system of first-order ordinary differential equations (ODEs): where the primes denote derivative with respect to τ = log a 3 .

A. Invariant set and monotonic functions
By following the procedures in the Appendix A, we can find the invariant sets.
, it follows from Proposition 10 that the sign of 3 R and the sign of nx φ + mx ψ are invariant. From equations (8)(9)(10)(11) 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 nonvanishing 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 2dimensional set Γ = x ∈ R 3 : 0 < x 2 φ − x 2 ψ < 1, y = 0 is invariant under the flow. 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 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 Then, W 4 is positive and strictly monotone along orbits of the invariant set S 4 .
Remark 3. In [31] it was devised the monotonic func- variant 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. Therefore, some global results can be found from the linear (local) analysis of equilibrium points.
B. The equilibrium points The system (8-11) admits six types of equilibrium points (with both z and y non-negative) denoted by M, F, C ± , P, M S, and CS. In tables I and II 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 3dimensional 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 phase-space. 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.
The equilibrium point P denotes a zero-curvature power-law inflationary model if δ < 1/3, corresponding to the equilibrium point with the same label investigated in [31]. It is worth noticing that the results of [31] concerning to this equilibrium point applies here only in the case of zero-curvature models with a dust background. It is non hyperbolic if δ = 1 3 or δ = γ 2 . P is a local inflationary attractor if either 0 < γ < 2/3, δ < γ/2 or 2/3 < γ < 2, δ < 1/3. If 2/3 < γ < 2 and 1/3 < δ < γ/2 then it is no longer a local inflationary solution nor an attractor for non-zero curvature models. However in this case the eigenvalues are (+, −, −, −) and then, the equilibrium point has a 3-dimensional stable subspace. If 0 < γ < 2/3 and γ/2 < δ < 1/3 the eigenvalues of P are (−, −, −, +). Hence, although an inflationary solution, P is not longer a local attractor. In this case the stable subspace is 3-dimensional.
As stated in [31], point P is the local attractor for zero-curvature 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 zero-curvature models without matter.
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 III.
defined in the phase space given by In table IV it is presented some information of the equilibrium points O, C ± , T, P from [31]. We have the following cases [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 defined on the phase plane In the figure 2 are presented some orbits in the phase plane of the system (19)- (20) 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  (20). We use the notation δ = m 2 − n 2 . It is assumed 0 < γ < 2. in the case m = 0.5, n = 0.6 whence P is a phantomlike 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].
Proposition 5. Let be the vector field X given by (14)(15)(16) 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 (14-16) 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: By applying theorem 2, we have that there exists a coordinate transformation x → y, such that (14-16) has normal form (22)(23)(24) 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. Now we integrate the truncated system (22), (23) (24) up to order N with initial condition (y 1 (t 0 ), y 2 (t 0 ), y 3 (t 0 )) = (y 10 , y 20 , y 30 ).
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 (29). The other components are given by 2. Normal expansion to third order for C− In this section we show normal form expansions for the vector field (14-16) 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 (14)(15)(16) then, there exist a transformation to new coordinates x → z, such that (14)(15)(16), defined in a vicinity of x * , has normal form where Proof. First step: transforming to the Jordan form Using a linear coordinate transformation the system (14-16) is transformed to where x stands for the phase vector with with Second step: simplifying the quadratic part By the hypotheses the eigenvalues λ 1 = 0, its eigenvectors form a basis of R 3 . The linear operator The eigenvalues Λ m,i for the allowed m, i are To obtain the normal form of (34-36) 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 (41) is the coordinate transformation required in theorem 2. By applying this theorem we prove the existence of the required constant c 2 . Now, let us calculate the value of c 2 . By applying the transformation (42) the vector field (37) transforms to we have i.e., The vector fieldX 3 (y) introduced above has the coefficients:X Third step: simplifying the cubic part After the last two steps, the equation (37) 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 nonresonant terms of third order we will consider the coordinate transformation y → z given by where The transformation (48) 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, where δ is a real value small enough. Therefore, the dynamics of (34-36) restricted to the unstable manifold, is given, up to order O(|z| 4 ), by 0, 0) . Then, the origin is the past attractor for an open set of orbits of (34)(35)(36).
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 (51) we have the expression for z 1 : Substituting (53) in (50)(51)(52) we have that the center manifold of the origin is given by the graph:
Therefore, we can study the reduced system defined in the compact phase space where the prime means derivative and ε = ±1 appears due to (56) cannot be solved globally , and it corresponds to the sign of H, where ε = +1 corresponds to ever expanding universe and ε = −1 corresponds to ever contracting universes.
The point P referred in table IV is recovered by setting ε = +1.

T
The point T referred in table IV is recovered by setting ε = +1.

(a)
The static solutions S1 and S2 exists for n > 0, −n ≤ m ≤ n.
In figure 3 (a), (b), (c), and (d) and (e) are presented some orbits in the phase plane of the system (57), (58), and (59) for ε = +1, that resembles the results found in [31]. The static universe configurations cannot by found in the framework of [31]. This illustrates the applicability of this formulation in compact variables.

III. 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 an alternative set of dynamical variables, albeit non compact, 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.

A. Negative and zero curvature models
In this section we extent the analysis presented in [83], since now it is included the possibility of the crossing through H = 0.
Let us introduce the following set of normalized variables: (D, U φ , U ψ , W, Ω), given by 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 Ω = Ω m D 2 ≥ 0 and Ω de D 2 ≡ U 2 φ −U 2 ψ +W 2 ≥ 0. Ω m , Ω de , Ω k denotes the fractional energy densities of dark matter, dark energy and "curvature". Then, from (64) 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 (62) 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 (62) 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 (62) are where we have defined the new independent variable The ODEs (68) 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 (64) and from the requirement of non-negativeness of the energy densities of the sources, that Ω = 0 and U 2 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 the set Ω = 0, U 2 φ − U 2 ψ = 0.

The dynamical system (68) is invariant under the transformation of coordinates
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 (71).
Remark 8. The existence of the monotonic function (72) 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 V it is displayed the conditions for existence and the eigenvalues of the linearized system around each equilibrium point of the system (68). In table VI 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 ± SF and ± CS 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 ± SF coincide as δ → γ 2 + and ± CS and ± SF coincide as δ → 1 3 + . Additionally, + SF (resp., − SF ) 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 (71) it is sufficient to characterize the behavior of the equilibrium points + SF, + CS, and + M S. The dynamical character of their counterparts in the "negative" branch − SF, − CS, and − M S is determined by the symmetry (71).
The equilibrium point + SF 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, + SF represents a de Sitter solution, and it is the attractor in the invariant set D = 1. If δ < 0 and 2 3 < γ < 2 then, + SF is the future attractor (in the full phase space) for (ever) expanding models and it represents a phantom field. Using the symmetry (71) we find that, if δ < 0 and 2 3 < γ < 2 then, − SF is the past attractor for contracting models and it represents a phantom field. If 0 < δ < 1 3 and 2 3 < γ < 2, then + SF 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 − SF 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 + CS 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 − CS 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 + CS 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 + SF at the bifurcation value δ = 1 3 and that the stability is transferred to the equilibrium point + CS (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 [×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, ±SF and ±M S an to the equilibrium sets ±K±, are, in each case, the same as those displayed, but the first from the left. 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 + SF 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 + SF. These equilibrium points are such that n U φ + m U ψ = 0. Hence, by the structure of the function (72), for ever expanding models (i.e., for 0 < D < 1), the late time behavior is determined by the local behavior of either + SF or + CS or + M S. By making use of the symmetry (71) it is possible to determine the early time behavior of contracting models (i.e, for −1 < D < 0) from the local behavior of either − SF or − CS or − M S. Combining the above linear analysis, with the information we have using monotonic functions we have the following B. Positive curvature models (k = +1) In this section we reproduce the results found by one of us in [83]. Introducing the normalization factor we haveD (i.e., at a singularity). This means that it is not possible thatD vanishes at a finite time. Now, using the normalized variables From the Friedmann equation we find and by definition By the restrictions (75,76), 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 It is useful to express some cosmological parameters in terms of our state variables.
(Ω m , Ω de , Ω k ) = Ω , 1 −Ω, Observe that the system (81,82) is invariant under the transformation of coordinates (83) 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 (83).

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 Q 2 0 → 1 or nÛ φ +mŨ ψ → 0 orΩ → 0. In the table VII 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 ± CS 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 (83), 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 (83).
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: (b) The future attractors are:

IV. 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.

A. 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 (86) 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 (85) 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 , VII: equilibrium points in the quintom model with k = +1. We use the same notation as in table V. If γ = 2 3 , the set of equilibrium points, S, (corresponding to static solutions) reduces to two sets of equilibrium points analogous to V1,2 (see table V). When the flow is restricted to the invariant sets Q0 = ±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.
the Poisson bracket. Thus, we can refer that the gravitational field equations described by the Lagrangian (86) 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 : Because of the nonlinearity of the differential equations it is not possible to write the analytical solution in closed-form by using well-known functions. Hence, we continue our analysis with the singularity analysis where the solution can be written by using Laurent expansions.

B. 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 (14), (15) and (16) 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 (14), (15) and (16) 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 14), (15) and (16) 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 (14), (15) and (16) passes the singularity test and its algebraic solution is given by the Right Painlevé Series (102), (103), (103). From the latter we can infer that the leadingorder term (95), i.e. the dominant behavior close to the singularity is an unstable solution.
We perform the same analysis for the complete dynamical system (8), (9), (10) and (11), where we find similar results. In particular the leading-order behavior is with constraint equation and resonances s 1 = −1, s 2 = 0, s 3 = 0 and s 4 = 0. Hence, we have the following proposition.
Proposition 9. The dynamical system described the four first-order equations (8), (9), (10) and (11), passes the singularity test with leading-order terms given by (108). The algebraic solution is given by the Right Painlevé Series (111) where x = (x φ , x ψ , y, z) and constraint equation (110). From the form of the Laurent expansions we refer that the leading-order term describe a past attractor for the dynamical system.

V. 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 + SF , 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 + CS, + 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 + SF 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 − SF , 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 − CS, − 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 + SF 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 + SF ). 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 lower-dimensional 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 La-grangian admits three linear-independent 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.
Acknowledgments GL thanks to Department of Mathematics and to Vicerrectoría de Investigación y Desarrollo Tecnológico at Universidad Católica del Norte for financial support. AP acknowledges financial support of FONDECYT grant No. 3160121. JLM acknowledges to Consejo Nacional de Ciencia y Tecnología (CONA-CYT) of Mexico for financial support throughout the Phd program on Water Sciences and Technology, Grant No. 003497.

Appendix A: Invariant sets and monotonic functions
In the following we shall investigate some invariant sets that can be trivially identified and we construct properly monotonic functions defined on them.
Recall the definition and results: Definition 1 (Invariant set). Let S ⊂ R n be a set. S is called an invariant set under the vector field if y ∈ S implies x(τ, y) ∈ S (where x(0, y) = y) for all τ ∈ R. If we consider the property valid for τ ≥ 0 we say that S is positively invariant. On the other hand, if the property is valid for τ ≤ 0 we say that S is negatively invariant.
Proposition 10 (Proposition 4.1, [79] p. 92). Consider the autonomous vector field (A1) with flow g τ . Let be defined a C 1 function Z : R n → R which satisfies Z ′ ≡ ∇Z · X(x) = αZ where α : R n → R is a continuous function. Then, the subsets of R n defined by x ∈ R n |Z(x) 0 are invariant sets for the flow g τ .
Definition 2 (Definition 4.8 [79], p. 93). Let g τ (x) be a flow on R n , let S be an invariant set of g τ (x) and let Z : S → R be a continuous function. Z is monotonic decreasing (increasing) function for the flow g τ (x) means that for all x ∈ S, Z(g τ (x)) is a monotonic decreasing (increasing) function of τ.
The existence of a continuous monotone function on an invariant set S simplifies the orbits in S significantly. As states the following proposition.
Proposition 11 (Proposition 4.2 ( [79], pp 93)). The existence of a continuous monotone function on an invariant set S implies that S contains no equilibrium points, periodic orbits, recurrent orbits or homoclinic orbits.
Theorem 1 (Monotonicity Principle (Proposition A.I in [78] developed in collaboration with M Glaum, cited as Theorem 4.12 in [79])). Let g τ (x) be a flow on R n with S an invariant set. Let Z : S → R be a C 1 (R n ) function whose range is the interval (a, b) where a ∈ R ∪ {−∞}, b ∈ R ∪ {+∞}, and a < b. If Z is decreasing on orbits in S, then for all x ∈ S, ω(x) ⊂ {s ∈S\S|lim y→s Z(y) = b} and α(x) ⊂ {s ∈S \ S|lim y→s Z(y) = a}.
Appendix B: Normal Forms for vector fields Let X : R n → R n be a smooth vector field satisfying X(0) = 0. We can formally construct the Taylor expansion of x about 0, namely, X = X 1 + X 2 + . . . + X k + O(|x| k+1 ), where X r ∈ H r , the real vector space of vector fields whose components are homogeneous polynomials of degree r. For r = 1 to k we write Observe that X 1 = DX(0)x ≡ Ax, i.e., the matrix of derivatives.
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 + N r=1 w r (y) + O(|y| N +1 ), where J is the real Jordan form of A = DX(0) and w r ∈ G r , a complementary subspace of H r on B r = L A (H r ), where L A is the linear operator that assigns to h(y) ∈ H r the Lie bracket of the vector fields Ay and h(y):