Revisiting dynamics of interacting quintessence

We apply the tools of the dynamical system theory in order to revisit and uncover the structure of a nongravitational interaction between pressureless dark matter and dark energy described by a scalar field $\phi$. For a coupling function $Q = -(\alpha d\rho_m/dt + \beta d\rho_\phi/dt )$, where t is the cosmic time, we have found that it can be rewritten in the form $Q = 3H (\alpha \rho_m + \beta (d\phi/dt)^2 )/(1-\alpha +\beta)$, so that its dependence on the dark matter density and on the kinetic term of the scalar field is linear and proportional to the Hubble parameter. We analyze the scenarios $\alpha=0$, $\alpha = \beta$ and $\alpha = -\beta$, separately and in order to describe the cosmological evolution we have calculated various observables. A notable result of this work is that, unlike for the noninteracting scalar field with exponential potential where five critical points appear, in the case studied here, with the exception of the matter dominated solution, the remaining singular points are transformed into scaling solutions enriching the phase space. It is shown that for $\alpha \neq 0$, a separatrix arises modifying prominently the structure of the phase space. This represents a novel feature no mentioned before in the literature.


Introduction
Recent cosmological observations indicate that our universe is currently undergoing an accelerated expansion phase. This has been confirmed by a wide variety of astronomical and cosmological data which includes measurements of high red-shift supernovae Ia (SNIa) luminosity, temperature anisotropies of Cosmic Microwave Background (CMB), Baryon Acoustic Oscillations (BAO) and Large Scale Structure (LSS) among others [1,2,3,4,5,6,7,8,9,10]. To explain such a late time acceleration in the context of general relativity it is necessary to assume the existence of a mysterious component with negative pressure broadly known as dark energy (DE). In the Lambda Cold Dark Matter (ΛCDM) model, the dark energy is described by a cosmological constant with equation of state (EoS) parameter ω = −1, which accounts approximately for the 70% of the total energy content of the universe [11,12,13]. It is still necessary to introduce an additional component dubbed Cold Dark Matter (CDM), which is postulated in order to increase the amount of structure formation needed to be in agreement with cosmological observations and represents around 25% of the cosmic inventory . scalar-tensor theories of gravity represent the prototypical way in which deviations from GR are modeled (see Refs. [22,23,24] for reviews). As an example, in Brans-Dicke gravity one introduces an additional scalar mode besides the metric tensor replacing the gravitational coupling G N by a point-dependent scalar field [25]. Alternatively, in the so-called f (R, G) theories the Lagrangian is a general function of the Ricci scalar R or the Gauss-Bonnet term G in the Jordan frame [26]. This gives rise to field equations with fourth-order derivatives and GR is recovered after the simplest choice of the function f (R, G) ∝ R. As a consequence of introducing an arbitrary function there is a lot of freedom to explain the observed data. Additionally, it is known that actions involving a finite number of power laws of curvature corrections and their derivatives can be considered as low-energy approximations to strings or supergravity theories giving rise to the so-called extended theories of gravity [27]. Finally, in the Dvali-Gabadadze-Porrati (DGP) braneworld model one assumes the existence of a five-dimensional (5D) Minkowski spacetime of infinite volume within which ordinary four-dimensional (4D) Minkowski spacetime is embedded. It is precisely the presence of additional dimensions that realizes cosmic acceleration through the leakage of gravity into the extraspace at cosmological scales. This latter model, however, is plagued by ghost instabilities that cast doubts upon its validity [28,29].
Alternatively, the lack of knowledge on the nature of the dark sector has motivated several approaches to unveil their physical properties. One of the simplest scenarios is assuming the existence of a minimally coupled scalar field φ with a self-interacting potential V (φ). This model arises from theories of gravity such as scalar-tensor theories and in the low-energy limit of string theories and has been the subject of interest due its ability to explaining various stages of the universe evolution [30,31,32,33]. The canonical scalar field dubbed quintessence resembles to the inflaton scalar field which was first proposed to explain the inflationary scenario which provides solutions to some issues of the big bang cosmology such as the initial singularity, flatness, horizon, homogeneity problems and the absence of magnetic monopoles [34]. Compared to other scalar fields such as k-essence, phantom and quintom, quintessence represents the simplest scenario without having theoretical problems such as the appearance of propagating ghost modes and Laplacian instabilities [35]. Its dynamical behavior is characterized by the equation of state parameter ω φ = P φ /ρ φ , where P φ and ρ φ denote its pressure and energy density respectively. For physically relevant cosmological scenarios the parameter is located into the interval −1 ≤ ω φ ≤ −1/3, where ω φ = −1 corresponds to the cosmological constant model. Quintessence models can be classified in two classes, freezing and thawing, depending on whether the equation of state decreases towards −1 or departs from it [36].
There is also the possibility that dark energy might interact with dark matter through a nongravitational cou-pling Q which is usually introduced at the level of the cosmological field equations. This represents an energy flow between the dark components and the sign of Q determines the direction of the energy transfer: for Q > 0 the matter fluid is giving energy to the dark energy fluid and vice versa for Q < 0. Notice that because of our current lack of knowledge about the nature of these two components, it would be imprudent to discard a nongravitational interaction between them. Although this kind of models was first proposed in order to alleviate the cosmic coincidence problem, it was found that they also improve predictions on LSS, BAO, CMB anisotropies, galaxy clusters and H(z) data among other cosmological and astrophysical experiments [37,38,39,40,41,42]. A wide variety of theoretical and phenomenological interacting scenarios have been proposed and investigated in the literature (see Ref. [43,44] for reviews and references therein). To name a few, theoretical aspects such as the possibility to construct an interacting Lagrangian from which the interaction term can be derived is analyzed in [45]. In [46] the authors study specific models of this class where they showed that cosmic chronometers and Type Ia supernovae data have a preference for interacting Quintessence models that lower H 0 relative to ΛCDM. In [47] physical limits on the equation of state parameter of the DE component non-minimally coupled with DM are examined in light of the second law of thermodynamics and the positiveness of entropy. The study of the growth of cold dark matter density perturbations in the nonlinear regime is performed in [48], and in [49] is shown that if the interaction between a quintessence field and cold dark matter is purely by momentum exchange, this generally leads to a dark energy sound speed that deviates from unity. Recently, assuming the dark energy component as a quintessence scalar field with Lagrangian function modified by the quadratic generalized uncertainty principle, in [50] the authors investigate the behaviour of solutions of the field equations for some interacting models of special interests in the literature. Even though current cosmological data are compatible with such energy transfer models, the evidence so far is not completely conclusive [51,52].
Since both the quintessence scalar field cosmology and the interacting dark energy models exhibit interesting phenomenological features, in the present work we perform a phase-space and stability analysis of the interacting scenario with exponential scalar potential and pressureless dark matter. Additionally, we compute some cosmological relevant quantities such as the dark energy density parameter, dark matter density parameter and the deceleration parameter. For the interacting term corresponding to a linear combination of the time derivatives of dark matter and scalar dark energy densities, Q = −(αρ m + βρ φ ), we analyzed the special cases depending on the value of the coupling parameters α = 0, α = β and α = −β separately. We will adopt dynamical system techniques which allow us to compute the equilibrium points and we focus on the attractor solutions that can give rise to late time acceleration. If the attractor solution exists, the evolution of several models corresponding to a wide range of initial conditions converges towards an unique asymptotic behavior.
The structure of this work is as follows. In section 2 we review the basic equations governing the cosmology of the interacting dark sector scenario. In section 3 we introduce the master equations for the dynamical analysis and compute helpful cosmological parameters. We establish the interacting quintessence model and we discuss the adequate choice of the variables of the phase space. In section 4 the dynamical system is solved, the fixed points are determined and their respective stability is analyzed. We will draw conclusions and discuss future perspectives in section 5. In appendix A we report the regions of existence of the critical points for the different cases studied in this work and in appendix B we will explicitly show the calculations around the conservation equations for the case of a nongravitational interacting scenario.

Cosmological equations
In this section, we briefly introduce the dynamics of the cosmic components for a non gravitational interacting model. We consider only two components in the cosmic inventory: the quintessence scalar field representing the dark energy and the cold dark matter described by a pressureless barotropic perfect fluid. Let us assume a flat Friedmann-Lemaître-Robertson-Walker(FLRW) metric with line element ds 2 = −dt 2 + a(t)[dr 2 + r 2 (dθ 2 + sin 2 θ dϕ 2 )], (1) where a(t) is the scale factor and t is the cosmic time. As a consequence of the interaction between the dark sector constituents the gravitational field equations become Here the dot represents the derivative with respect to the cosmic time, κ 2 ≡ 8πG, with G the gravitational coupling constant, φ is the scalar field, ρ m denotes the dark matter energy density and ρ φ and P φ represent the energy density and the pressure of the scalar field respectively The function V (φ) is a self-interaction scalar potential and in this work we assume an exponential potential of the form V (φ) = V 0 exp[−κλφ], where V 0 > 0 is a constant with dimensions of mass and λ is a dimensionless constant.
This corresponds to the simplest example of quintessence scalar field and can be easily justified from high-energy phenomenology [53]. The cosmological dynamics of the exponential potential is captivating because of the appearance of accelerated solutions which can be employed to explain both the inflationary stage and the late time dynamics [54,55]. Besides, the exponential potential has the interesting property of generating tracking solutions, i.e., for an appropriate choice of the parameter λ, the quintessence field evolves like radiation during the radiation-dominated era, and like matter during the matter-dominated era. This family of cosmological models have been extensively discussed in relation with early time inflation, high energy physics and late time accelerated scenarios [56,57,58,59,60,61,62,63,64].
Assuming the existence of an additional non gravitational interaction Q which is introduced at the level of the cosmological field equationṡ where the sign of Q determines the direction of the energy transfer: for Q > 0 the matter fluid is giving energy to the dark energy fluid and vice versa for Q < 0. It is important to mention that, in order to satisfy the local energy conservation, equations (6) and (7) are not independent due to the Bianchi identities as we show in appendix B. In this work we will restrict the discussion to quintessence models where dark energy is assumed to be a scalar field with self-interacting exponential potential, however, we can assume the existence of more complicated potential terms [65,66,67,68,69,70]. For dynamical systems studies of interacting dark energy as a perfect fluid we can find in the literature [71,72,73,74,75]. Finally, the evolution for the scalar field is given bÿ

Phase space variables
The application of the dynamical system theory is specially useful when one deals with scalar-field cosmological models [76,77,78,79]. It should be mentioned that, as far as we know, reference [80] is a pioneering study on the application of dynamics systems to cosmology. From the dynamical systems tools one may obtain very useful information on the asymptotic dynamics of the system which is characterized by: i) source critical points which may be pictured as past attractors, ii) saddle equilibrium configurations that attract the phase space orbits in one direction but repel them in another direction, iii) attractor solutions to which the system evolves for a wide range of initial conditions, or iv) limit cycles, among others.
In order to trade the system of second order equations (7) and (8) by a system of autonomous ordinary differential equations one has to choose a suitable set of variables.
In general, there are many possible ways to achieve this task. The most common one is to consider the normalized variables introduced in [81]: Here we are assuming that only expanding cosmologies arise: H ≥ 0 (with Y ≥ 0). The constraint (2) written in terms of the set of normalized variables takes the form where Ω m and Ω φ are the density parameters usually defined for the dark matter and scalar field respectively. Hence, the physically meaningful phase space corresponds to the region Using this set of variables (9) and the Friedmann constraint (10) the system of equations which governs the dynamics reduces to the following set of autonomous equations where the prime denotes the derivative with respect to the logarithm of the scale factor. It is important to mention that, in general, the system of equations (12) -(13) is not closed unless the coupling function Q can be expressed in terms of the variables (9). In this work, we consider the coupling function which was incorrectly studied in [82] and [83] as we show at the end of this section. Replacing (14) in (6) -(7) and after some algebra, the system of equations takes the form We can algebraically solve forρ φ andρ ṁ and we can write the interaction term (14) into the form Then we have shown that the interaction function (14) is equivalent to an interaction term Q lineally proportional to the Hubble parameter and a lineal combination of the dark matter density and the kinetic term of the scalar field.
In order to proceed to the phase-space analysis it is necessary to compute the functionsḢ/H 2 and Q in terms of variables X and Y . From equation (3) in addition with (4) and (5), we obtaiṅ (20) Finally, with the help of (9), (19), (20), we derive the first order dynamical system from equations (12) and (13) as The choice α = β = 0 represents the non-interaction scenario studied in [81]. For the case β = 0 it can be shown that making the transformation 3α 1−α → α we can reproduce the model II analyzed in [84].
The deceleration parameter which is defined as one of the geometrical parameters through which the dynamics of the universe can be quantified is depicted by while the effective equation of state parameter ω ef f ≡ Ptot ρtot can be written as Finally, we notice that, in the case α = 0 and the limit X → 0, Y → Y 0 with Y 0 = ±1, the last term of the right hand side of (21) diverges and therefore the system of equations (21)- (22) does not satisfy the fundamental existence and uniqueness theorem for nonlinear Ordinary Differential Equations Systems because it is not continuously differentiable (See page 74 of the reference [85]).
this shows that the vertical line X = 0 is a separatrix in the compact phase space (11), namely, the dynamics of the region X > 0 is completely disconnected causally from the region X < 0. Actually, in a neighbourhood of X = 0 there exist two different trajectories with the same initial or end condition: in the case α/ is the initial condition for two different trajectories which depart from it (for example, see bottom panels of Figs. 6 and 9; by the contrary, in the case α/(1 − α + β) < 0, the point is the end point for two different trajectories (for example, see top panels of Figure  6 and Figure 9. As it was previously mentioned, the interaction coupling (14) was studied in references [82], [83] where a corresponding mistaken dynamical system was analyzed providing wrong results in the critical points found and the subsequent stability analysis; specifically, the first of equations (9) of reference [82] and equation (18) of reference [83] are wrong because, in both, the interaction term is missing (which provides precisely the term that is proportional to the inverse of the variable X in the right hand side of equation (21) in this work).

Dynamical analysis, critical points, and stability
This section is devoted to analyze the cosmological dynamics of the system of cosmological equations (7) and (8) by the system of autonomous ordinary differential equations (ODE-s) (21) and (22) in the formẋ = f (x). Here x is called a point in the phase space and f corresponds to the column vector of the autonomous equations. A critical (or equilibrium) point x c , is a point in the phase space that satisfies the condition f (x c ) = 0. In order to determine the stability properties of the system we expand around x c as x = x c +u, with u the column vector of the perturbations. Therefore, for each critical point we expand the perturbation equations up to first order asu = IM u, where the matrix IM contains the coefficients of the perturbation equations. Finally, the eigenvalues of IM are evaluated for all critical points in order to determine its type and stability.
For the dynamical system defined by the equations (21) and (22), there are six critical points, these are reported in Table 1. We have defined the following equations: To determine the existence of critical points A + , A − , B + and B − we use the constraint given by equation (11). Note that the X-component for the critical points C + and C − is antisymmetric under λ → −λ, this is X(α, β, −λ) = −X(α, β, λ) while the second restriction 0 ≤ X 2 + Y 2 ≤ 1 holds. For the critical points A + , A − , B + and B − it is easy to find their stability since there is no dependence on the parameter λ. On the other hand, for the critical points C + and C − , this represent a more complicated task and for this reason we analyze only some special cases in order to simplify the analysis.

Scenario α = 0
For a vanishing coupling constant α, the interacting kernel reduces to while the functions ∆ and Γ take the form In this scenario we have found five critical points (the critical points A + , A − , B + and B − reported in Table 1, reduce to O, D + and D − ) reported in Table 2. The existence conditions, stability, acceleration, and ω ef f are reported in Table 3. The region of existence of the critical point E + is the region 1 reported in appendix A and we show this in Figure 1. The region of existence of the critical point E − is region 2 reported in appendix A and this region is not bounded in the parameter λ, for this reason, we only report part of this and we show in Figure 2 the region considered.
Here, q < 0 means that there is acceleration and we can see in Table 3 that only for the critical point E − we have acceleration and this is shown in Figure 2, where below the black dotted curve we have acceleration and in the other case we have deceleration.
The critical points of the dynamical system for the choice α = 0, as well as their stability properties, are Table 1. Critical points of the dynamical system. Table 2. Critical points of the dynamical system for α = 0. Table 3. The physically meaningful critical points of the autonomous system for the case α = 0.

Point Existence
Stability Saddle in all region i in Figure 1 No  listed and briefly discussed below. For a couple of illustrative scenarios see Figure 3.
(i) Point O: The matter dominated solution exists for a coupling constant β = −1 and it is independent of the specific form of the self-interacting potential. Here the effective equation of state parameter vanished (ω ef f = 0) and therefore there is no acceleration (q = 1/2). For −1 < β < 1 this point behaves as saddle, otherwise it is unstable. In the non-interacting scenario this solution behaves always as saddle, therefore the chance that this point can be related to an origin of some trajectories in the phase space is due to the presence of a nongravitational interaction.
(ii) Point D + : The dark energy scaling solution exists for all values of the parameter λ and for 0 ≤ β < 1. The special case β = 0 denotes an universe dominated by the scalar field kinetic energy (Ω m = 0, X = 1 and Y = 0) and the limit β → 1 corresponds to a matter dominated universe (Ω m → 1, X → 0 and Y = 0). The effective EoS parameter is depicted by ω ef f = 1−β 1+β ∈ (0, 1], and then the decelerated parameter is non-negative corresponding to a decelerated solution. For λ < 6 1−β 2 the solution is a past attractor and for λ > 6 1−β 2 it behaves as saddle and therefore it cannot be a late-time state of the universe. (iii) Point D − : The scaling solution does not depend on the parameter λ, but it is still required that 0 ≤ β < 1. The limit β = 0 corresponds to a stiff matter universe (Ω m = 0, X = −1 and Y = 0), and β → 1 denotes a matter dominated universe (Ω m → 1, X → 0 and Y = 0). The equation of state parameter is always positive ω ef f = 1−β 1+β ∈ (0, 1], therefore the solution is decelerated. This point behaves as saddle for λ < − (iv) Point E + : Exists for the region 1 reported in the appendix. It behaves always as saddle and thus it cannot attract the universe at late times. For β = 0 the matter energy density vanished (Ω m = 0) while X = λ √ 6 and Y = 1 − λ 2 6 , which corresponds to the scalar field dominated universe. For β = 1 3 and λ = 27 4 the universe has the components Ω m = 1 2 , X = 1 2 and Y = 0. In general, Ω m , X, and Y never vanish simultaneously, despite this, the scalar field kinetic energy never dominates.
(v) Point E − : The scaling solution exists for the region 2 reported in the appendix. This point can be unstable (node and spiral), a centre or either stable (node and spiral). This solution is accelerated if the parameters lie in the region λ < 2−4β (1+β) and −1 < β ≤ 0 and then it can be the late-time state of the universe. As an example, for β = − 1 2 and λ = 1 we have the following quantities X ≈ 0.187, Y ≈ 0.939, Ω m ≈ 0.082 and q = −0.270. In the non-interacting scenario this point is either a stable node or a stable spiral.
In this case, the interaction is where α = β ≡ ζ. The critical points are reported in Table  4, and the ∆ and Γ functions are: The existence, stability, acceleration (q < 0) and ω ef f are reported in Table 5. The region of existence of the critical points H + and H − are called region 3 and region 4 respectively, this is reported in the appendix A. The region of existence of the critical point H + is shown in Figure 4. For the critical point H − , the region of existence is not bounded in λ, for this reason, the stability of the region is shown in Figure 5. For Table 7 we define w 1 (ζ) and w 2 (ζ): The critical points of the dynamical system for the choice α = β, as well as their stability properties, are listed and briefly discussed below. For a couple of illustrative scenarios see Figure 6.
(iv) Point G − : This solution exists for all values of the parameter λ, and 0 ≤ ζ < 3 − 2 √ 2 is required. The case ζ = 0 corresponds to a scaling solution, while ζ = 0 is related to a matter dominated solution Ω m = 1, where the dynamical variables X and Y identically vanished. The cosmological parameters read Ω m → 0.58, X → −0.64 and Y = 0 in the limit ζ → 3 − 2 √ 2, therefore the scalar field never dominates. The equation of state parameter is located into the interval ω ef f ∈ [0, 0.41), and hence the solution is always decelerated.   6. Phase portrait of the autonomous system of ODE (21) and (22) for the case α = β and the specific choice (β = −0.01, λ = 2) -top panels-and (β = 0. − 5, λ = 2.5) -bottom panels-. Here the gray line denotes the separatrix X = 0 which exists whenever Y = 1 and divides the regions where the kinetic term of the scalar field is positive or negative. It can be seen that the trajectories end at it for α/(1 − α + β) < 0 or depart from it if α/(1 − α + β) > 0. Due to an optical illusion one could assume that some trajectories can cross the separatrix, and then the behavior of the trajectories around the separatrix is shown in the right-hand side. In the first scenario, Point F+ and Point F− denote unstable nodes, Point G+ and Point G− correspond to saddle nodes, Point H+ describes an accelerated stable solution, which can be of cosmological interest. For the latter scenario, Point F+, Point F−, Point G+ and Point G− vanish, leaving Point H+ describing a saddle node and Point H− which is associated to a stable spiral. This scenario is decelerated.

Scenario α = −β
For this case the interaction is where α = −β ≡ η, same interaction was studied in [83] but again with wrong equations. The critical points are reported in Table 6. The ∆ and Γ functions are: The existence, stability, acceleration (q < 0) and ω ef f are reported in Table 7. The region of existence of the critical points K + and K − are reported in the appendix A label by region 5 and region 6 respectively. For Table 7 we define the following functions where δ = 1 − 8η + 30η 2 − 72η 3 + 81η 4 , furthermore (η − 2η 2 ) = 0.
The critical points of the dynamical system for the choice α = −β, as well as their stability properties, are Table 6. Critical points of the dynamical system for α = −β ≡ η. Table 7. The physically meaningful critical points of the autonomous system for the case α = −β ≡ η.
(i) Point I + : Exists for η ≤ 0, depending on the λ value there is a saddle point or unstable point, so that their phenomenological properties remain the same independently of the potential. For the special case η = 0 the universe is dominated by the scalar field kinetic energy (X = 1, Ω m = 0 and Y = 0). When η is very large negative is also dominated by scalar field kinetic energy (X → 0.70, Y = 0 and Ω m = 0.5 . For this critical point there is no acceleration.
(ii) Point I − : Exists for non-positive η, depending on the λ value there is a saddle point or unstable point, so that their phenomenological properties remain the same independently of the potential. The scenario η = 0 the universe is dominated by the scalar field kinetic energy (X = −1, Y = 0 and Ω m = 0). When η is very large negative is also dominated by scalar field kinetic energy (X → −0.70, Y = 0 and Ω m = 0.5 . For this critical point there is no acceleration. (iii) Point J + : Exists for η ≥ 0 (η = 1 2 ), there is a saddle point, stable point or unstable point, depending on the λ value. For the special case when η = 0 the universe is matter dominated (X = 0, Ω m = 1 and Y = 0). For the case when η tends to positive infinity we have X → 0.70, Ω m → 1 2 and Y = 0. There is no acceleration for this crit-ical point J + .
(v) Point K + : Exists for the region 5 reported in the appendix A. One part of this region is shown in Figure 7. For f 3 (η) < λ < f 2 (η) and η ≤ 0 this is a saddle point.
For λ > f (η) and η > 1 this point is unstable and is very interesting to note that there is a region where is accelerated as can be seen in Figure 7.

Conclusion
In this work we have performed a dynamical analysis in a spatially flat, homogeneous and isotropic spacetime for a nongravitational interaction scenario between pressureless dark matter and a quintessence scalar field with selfinteracting exponential potential. Here we have considered a coupling function which is of cosmological interest, namely, Q = −(αρ m + βρ φ ), and we have studied the following scenarios α = 0, α = β and α = −β, separately. We have found that this form of the interacting term can be rewritten as in equation (19), so that its dependence on the matter density and on the kinetic term of the scalar field is lineal and directly proportional to the Hubble parameter Q = 3H αρm+βφ 2 1+β−α . In order to describe the cosmological evolution for each solution we have calculated various observables such as the effective equation of state parameter, the DE and DM dimensionless density parameters, the effective EoS parameter and the deceleration parameter. For every case, we have found the existence of singular points which can be related to relevant epochs in the history of the universe such as the matter dominated solution, the stiff matter universe, the scalar field dominated solution and the scaling scenario. Besides, it is necessary to mention a relevant aspect that has been overlooked in previous works [82,83,84]: the existence of a separatrix whenever the coupling α is a non-vanishing constant. This completely modifies the structure of the phase portrait since divides it into two regions which are causally disconnected according to the sign of the kinetic term of the scalar field. In fact, at point X = 0 whenever Y 0 = ±1, the dynamical system fails to be continuously differentiable and therefore the system of cosmolog- Fig. 9. Phase portrait of the autonomous system of ordinary differential equations (21) and (22) for the case α = −β and the specific choice (β = 2, λ = 2.5) -top panels-and (β = 0.25, λ = 1.75) -bottom panels-. The gray line corresponds to the separatrix X = 0 which exists for Y = 1. It can be seen that the trajectories end at it for α/(1 − α + β) < 0 or depart from it if α/(1 − α + β) > 0. The behavior of the trajectories around the separatrix is shown in the right-hand side. In the upper scenario, Point I+ and Point I− correspond to unstable nodes and Point K− denotes a saddle point. It is not relevant from cosmological considerations since there is no stable solutions. In the lower case, Point J+ and Point J− behave as saddle and Point K− corresponds to an accelerated stable solution which can be the late time state of the universe.
ical equations does not satisfy the fundamental existence and uniqueness theorem for nonlinear ordinary differential equations. Therefore there can be two different trajectories with the same initial or final condition. This represents one of the main results of our analysis.
In the case where the nongravitational interaction is absent α = β = 0, we recovered the results of standard quintessence. Moreover, when the energy transfer between dark components is present we found that some critical points may survive or they may completely disappear depending on the model parameters. We have found that except for the matter dominated solution, the remaining critical points correspond to scaling solutions where neither dark energy nor dark matter dominates. This represents a modification with respect to the non interacting scenario analyzed in [81] where the authors described the phase space of the universe with four critical points: (i) the ordinary dark matter dominated solution (always saddle), (ii) the kinetic energy dominated scenario described by a stiff fluid EoS which can be saddle or unstable, (iii) the scalar field dominated solution representing a saddle or stable node and (iv) the scaling solution which is always stable (node or spiral).
For the case α = 0 and β = 0, there are five critical points denoted by Point O, Point D + , Point D − , Point E + and Point E − . Critical point O exists for β = −1 and behaves as unstable node or saddle. Point D + and Point D − have slightly more complicated stability conditions and denote as scaling solution whenever 0 < β < 1. The limit β → 1 corresponds to the matter dominated solution and the stiff matter solution is recovered for β = 0. This represents a modification to the mistaken results found in [82,83], where the authors determine that these solutions always correspond to saddle points. Finally, although Point E + is always saddle, the true richness of the system is found at point E − , since according to the values of the parameters β and λ the solution can be a stable node or a stable spiral, as well as an unstable node or unstable spiral o even a centre point for the case of one o more vanishing eigenvalues. It is important to mention that all solutions, except for Point E − , are decelerating solutions.
The second scenario α = β, possesses six critical points, namely Point F + , Point F − , Point G + , Point G − , Point H + and Point H − . Here, Point F + and Point F − can be saddle or unstable and thus they cannot be the latetime state of the universe. These represent decelerating solutions where the quintessence tracks the dark matter behavior. The limit Q → 0 corresponds to the ordinary stiff matter scenario. Point G + and Point G − represent non-accelerating solutions which exist for non-negative η. For η = 0, the scaling solution behaves as stable node and then can be relevant at late-times. The matter dominated solution is retrieved in the limit Q → 0 and it is always saddle. Point H + is always saddle and denotes only a transient epoch in the cosmological history. Finally, Point H − is the only accelerated solution and the stability properties indicate that this can be a stable node or a stable spiral, as well as an unstable node or unstable spiral o even a centre point.
Finally, for α = −β, there are six critical points, namely Point I + , Point I − , Point J + , Point J − , Point K + and Point K − . Here, Point I + and Point I − can be saddle or unstable and thus they cannot be the late-time state of the universe. These represent decelerating solutions where the quintessence tracks the dark matter behavior. The limit Q → 0 corresponds to the ordinary stiff matter scenario. Point J + and Point J − represent non-accelerating solutions which exist for non-negative η. For η = 0, the scaling solution behaves as stable node and then can be relevant at late-times. The matter dominated solution is retrieved in the limit Q → 0 and it is always saddle. Point K + is always saddle and denotes only a transient epoch in the cosmological history. Finally, Point K − is the only accelerated solution and the stability properties indicate that this can be a stable node or a stable spiral, as well as an unstable node or unstable spiral o even a centre point.
We close this work by mentioning that, even though the existence of a nongravitational interaction between the dark energy and dark matter does not generate the appearance of new critical points, it does greatly modify the stability of the solutions at background level. There could still be the case that the interaction could leave their signatures on observables related to cosmological perturbations such as the density fluctuations and the power spectrum. Although such an investigation lies beyond the scope of the present paper, it could be interesting to investigate the relevance of this interacting scenario by confronting with cosmological observations. Region 6:

B Energy-momentum conservation
In this Appendix we will show the calculations around equations (6) and (7). In the standard cosmological model dark matter and dark energy are considered to be uncoupled with separately conserved energy-momentum tensors. General covariance requires the conservation of their sum, so that ∇ ν G µν = ∇ ν (T m µν + T φ µν ) = 0.
The standard way of coupling two interacting matter components consists of adding a nonvanishing current F µ = F µ (ρ m , ρ φ , u α , ∇ α u α , ∇ α ρ m ) to the right-hand side of the conservation equations, such that which guarantees the overall energy-momentum conservation. The projection onto the orthogonal 4-velocity u µ defines the interacting term Q where the four-velocity u µ satisfies the condition u µ u µ = −1, and therefore u ν ∇ µ u ν = 0. For a flat Friedmann-Lemaître-Robertson-Walker (FLRW) background spacetime metric and a pressureless matter component, equations (44) take the forṁ