Phase space analysis and singularity classification for linearly interacting dark energy models

In this paper, applying the Hartman–Grobman theorem we carry out a qualitative late-time analysis of some unified dark energy-matter Friedmann cosmological models, where the two interact through linear energy exchanges, and the dark energy fluid obeys to the dynamical equation of state of Redlich–Kwong, Modified Berthelot, and Dieterici respectively. The identification of appropriate late-time attractors allows to restrict the range of validity of the free parameters of the models under investigation. In particular, we prove that the late-time attractors which support a negative deceleration parameter correspond to a de Sitter universe. We show that the strength of deviation from an ideal fluid for the dark energy does not influence the stability of the late-time attractors, as well as the values of all the cosmological parameters at equilibrium, but for the Hubble function (which represents the age of the universe). Our analysis also shows that a singularity in the effective equation of state parameter for the dark energy fluid is not possible within this class of models.


Introduction
Late-time interactions between dark energy and dark matter do not violate current observational constraints [1,2]. In particular, it has been proposed that energy flows between the two dark components of the Universe can alleviate the "coincidence problem": why do we live in a special epoch a e-mail: Mohsen13.8b@gmail.com b e-mail: danielegregoris@libero.it c e-mails: khurshudyan@yandex.com; khurshudyan@ustc.edu.cn of the evolution of our Universe in which the amounts of dark energy and dark matter are of the same order of magnitude? [3][4][5]. Furthermore, they may mitigate as well the Hubble tension [6], and they have been investigated in light of the 21-cm line excess at cosmic dawn by one of us [7]. On the other hand, the picture of the dark energy fluid as a cosmological constant term entering the Einstein equations is problematic due to many reasons: it violates the causality principle, its adiabatic speed of sound is ill-defined, there is no theoretical framework in terms of elementary particles theories which can account for its physical properties with in particular its cosmological value being different from the one estimated from vacuum energy by 120 orders of magnitude [8]. These considerations have been making impossible to develop appropriate technologies for a direct detection of dark energy, making the "dark energy problem" to be considered the most urgent open question to answer in a recent survey conducted among the fifty most prominent cosmologists [9]. In fact, the idea itself that dark energy should be a physical component of our Universe has been constantly challenged suggesting that it may indeed constitute just an interpretative aspect of the current standard model of cosmology which, for example, neglects the gravitational role that astrophysical objects like galaxies, clusters, voids, and filaments have on the large scale evolution of the Universe [10][11][12][13].
For addressing the limits of the description of the dark energy as a cosmological constant, a number of fluid models based on non-ideal equations of state have been proposed and tested against observational datasets with the twofold purpose of accounting for its thermodynamical properties and modifying the value of its adiabatic speed of sound in such a way that it will not violate the causality principle [14][15][16][17][18][19]. Moreover, a time evolution of the equation of state for dark energy does not violate current observations [20][21][22], and various possibilities for its parametric reconstruction have been investigated [23][24][25]. Exact solutions to the field equations become scarce when the picture of dark energy with a cosmological constant is abandoned, and thus dynamical system techniques constitute a valuable tool for predicting the late-time state of the evolution of the models. For example, they have been employed in Galileon gravity [26], Horava-Lifshitz gravity [27], modified teleparallel models [28], Einstein-Aether model [29], inflationary models [30], for the study of the evolution of instabilities in a Friedmann universe [31], and by one of us to various modified Chaplygin gas scenarios [32] in a Friedmann cosmology including interacting cases [33], ghost energies regimes [34], and frameworks involving evolving both G and [35].
As a subsequent step, dynamical system techniques are usually adopted as well for studying the singularity properties of the cosmological model in hand [36][37][38][39]. In fact, the results known under the name of singularity theorems prove that under certain conditions on the matter content of the Universe, a singular state characterized by the divergence of certain physical quantities can be reached even in a finite amount of time during the evolution [40][41][42][43]. However, the mathematical demonstration of such results are usually based on the modeling of dark energy with a cosmological constant, and the roles of nonideal fluids have not been fully accounted for.
In this paper we investigate the asymptotic late-time evolution of a Friedmann universe whose matter content is a mixture of dark energy and pressure-less dark matter adopting three descriptions proposed in [44], but allowing also for interactions between these two constituents modeled as a linear flow of dark energy, or dark matter, or of their sum, respectively. We will adopt dynamical system techniques which allow us to compute the attractor equilibrium points in terms of elementary functions showing explicitly that they can account for an accelerated expanding Universe, and how this physical requirement permits to restrict the range of validity of the free parameters entering our models. The singularity type of these physical attractors will be presented. Then, the nine cases (three possible modelings for dark energy times three possible modelings for dark interactions) will be compared and contrasted with each other in light of their late-time properties.
The plan of the work is as follows: In Sect. 2, we introduce three different cosmological interacting models between dark energy and dark matter with three different nonideal equations of state for dark energy and provide a formalism to transform their Friedmann equations into an autonomous dynamical system. In Sect. 3 we perform the phase-space analysis for each interacting model and provide a summary of the obtained results commenting on the patterns which can be recognized. Section 4 classifies the possible singularities arising in these models by discussing the corresponding cosmological implications. Finally, in Sect. 5 we summarize our results and make our conclusion. In Appendix A a brief review will be provided about the approach of dynamical system theory, and the dynamical autonomous equations for each model that we consider will be derived. Moreover, Appendix B is devoted for the study of transition between decelerating and accelerating epochs for each of our model.

Basic equations of our cosmological model
In this section we will introduce the notations that we will follow throughout the paper and review the basic properties of the Friedmann cosmology under investigation. We adopt units such that 8π G = 1 = c. ρ de denotes the dark energy density, ρ dm the dark matter density, p de the dark energy pressure, while dark matter is assumed to be pressure-less. e is the number of Nepero. Let H =ȧ a be the Hubble function, H 0 its present value, a the scale factor of the universe and a 0 its present value. The dimensionless energy-matter parameters on which we will base our dynamical system analysis are: An over dot stands for differentiation with respect to the coordinate time t, and Q will be the interaction term between dark matter and dark energy. The parameter for the effective dark energy equation of state is constrained by current astrophysical observations as [20][21][22]: We assume that the geometrical properties of the Universe can be accounted for by the flat Friedmann metric [45] The matter content is described by the stress-energy tensor T μν = (ρ + p)u μ u ν + pg μν , where u μ = δ μ t is the fourvelocity of the reference observer. Introducing the Einstein tensor G μν , the equations governing the dynamics of the model are the Einstein field equations G μν = T μν and the Bianchi identities T μν ;ν = 0, in which a semicolon denotes a covariant differentiation. Thus, the basic equations can be reduced to the Friedmann equation, the acceleration equation and the equations for the conservation of dark energy and dark matter: In terms of the dimensionless matter parameters (1), the Friedmann equation can be cast into the constraint and we can observe that its form does not depend on the modeling of the interaction term coupling dark energy and dark matter that we can choose. Furthermore, physical requirements impose the energy densities to be non-negative, and thus we have the bounded variables 0 ≤ x ≤ 1 and 0 ≤ z ≤ 1 for which dynamical system techniques can be used [36], while H is the appropriate monotonic function for the Friedmann model. The particular cases x = 1 and z = 1 correspond to a dark energy and a dark matter dominated universe respectively. The deceleration parameter reads as [46,47]: which must be negative for cosmological meaningful models [20][21][22]. Now we need to derive also the remaining evolution equations in terms of appropriate variables for exploiting the techniques provided by the theory of dynamical systems. For this purpose, let a prime denote differentiation with respect to the number of e-folds of the Universe N = ln a [48]. We compute the equations of the model using the chain rule The evolution eq. of the Hubble function does not depend on the modeling of the interaction term and it reads as: The evolution equations for the matter parameters become from which we can see that the compatibility condition x = −z which follows from the Friedmann Eq. (8) automatically applies implying that we can develop the dynamical system analysis in terms of x or z without any difference [36]. In our analysis we will focus our attention on the following three phenomenological types of interactions between dark matter and dark energy [49][50][51][52][53][54][55]: with the physical constraints on the free parameters −1 ≤ b ≤ 1 and −1 ≤ γ ≤ 1. By considering both positive and negative signs in the model parameters, we allow both the dark energy to flow into dark matter and viceversa; moreover the strength of such flows are proportional to the amounts of dark energy, or dark matter, or a combination of the two, respectively, for the three cases. We choose to consider these types for the interaction terms in the dark sector because we assume that the interaction represents a small contribution in the whole evolution of the energy budget of the Universe, and in fact these terms can be interpreted as a first order Taylor expansion. Moreover, they are based on the assumption that the propagator of the dark particles is energy dependent. These choices have been shown to be a valid tool for alleviating the coincidence problem [56]. Letting y = y(H, x) be the equation of state for the dark energy fluid we can derive the second differential equation entering the dynamical system: The two equations which will constitute our dynamical systems are thus provided by (x , y ) and the equilibrium points for which (x = 0, y = 0) will be (x eq , y eq ). We stress the fact that the Hubble function H is not a good variable for applying the algorithm of the theory of dynamical systems since it is not dimensionless and it is not necessarily bounded: it must be eliminated from the dynamical equations by inverting the equation of state for the dark energy as H = H (x, y). Then, the dynamical system analysis will let us estimate this cosmological parameter as H eq = H (x eq , y eq ). In the current paper we will consider different equations of state for dark energy which are known under the name of Redlich-Kwong [57], Modified Berthelot [58], and Dieterici [59] which read respectively as: We are interested in such equations of state for the modeling of the dark energy fluid because of their applicability in cosmology [44]. We refer the reader to the appendix of such paper for a review of their microscopic foundation and of their main thermodynamical properties. Here we just want to mention that the positive parameter α quantifies the deviations from a perfect fluid behavior expressed in units of the present day critical density 3H 2 0 because in all the above cases we reduce to a linear equation of state p ∼ ρ in the limit α → 0. On the other hand, the parameter β, which must be negative for accounting for a dark-energy-like fluid (as easily recognized from the limit at small α), can be removed in favor of cosmologically meaningful quantities as respectively, where q 0 denotes the present-day value of the deceleration parameter q = −ä ȧ a 2 and x 0 = αx(z = 0), and z 0 = z(z = 0),z denoting the redshift. These equations of state have been proposed as improved and more realistic approaches to fluidodynamics than the Van der Waals equation of state still accounting for the internal forces acting between the molecules constituting a gas. They have been used for modeling liquid-vapor phase transitions with the most important difference from the Van der Waals theory being that the attractive force term has been allowed to be temperature-dependent avoiding the oscillations in the isotherm curves below the critical temperature, and they constitute an analytical alternative to the use of virial expansions [60,61]. For applying dynamical system techniques in cosmology, these equations of state must be recast as: From the equations of state (24)-(26) we can compute respectively: in which the double sign corresponds to the cases of expanding or contracting universes respectively, and where W (X) denotes the Lambert function which is the inverse of Xe X . Once evaluated at the equilibrium point, the requirement of having a positive argument for the square root will impose a further restrictions on the numerical values of the free parameters of the models we are considering. We note that in the evolution equation for y below, H will appear always squared, so that its sign does not affect our analysis.

Phase space analysis
In this section we perform phase space analysis, then reporting the values of the equilibrium points, exhibiting the restrictions they impose on the parameters of the models α, β, b, and γ , and discussing their stability and whether they are cosmologically acceptable for the nine combinations equation of state for the dark energy -interaction term. Together with a basic review on dynamical system theory, the autonomous dynamical systems equation for each of these combinations have been derived in the Appendix A.
3.1 Redlich-Kwong: interaction term Q 1 The dynamical system to consider is formed by the two Eqs.
(A7) and (A9). In this case we get three mathematical equilibria (x eq , y eq ) which read as: However only the first point is relevant for cosmology because the second and third imply a divergent and a vanishing H eq respectively. For such a physical attractor we get: (H eq , z eq , q eq , w eq ) From the values of such physical quantities we can restrict − 1 5 ≤ b ≤ 0, where the case b = 0 corresponds to a single fluid pure dark energy pictured as a cosmological constant universe. Moreover the relation (b − β − 1)(b + β − 1) < 0 must hold. Since the second factor is always negative for our choice of the model parameters, we can further restrict The Jacobian matrix associated to the dynamical system when specified to such equilibrium point reads as: implying Thus, taking into account the range for the free parameters of the model, we get DetJ > 0, which means that our equilibrium point is a a node (spirals are not possible because (TrJ) 2 4 − DetJ > 0), and we can understand that it is sta- An example of the phase portrait for such a system is displayed in Fig. 1, panel a, for the numerical choice of the free parameters as b = −0.1, γ = 0.5, β = −1.2, and α = 1.5; the equilibrium point is denoted with a circle.

Modified Berthelot: interaction term Q 1
The dynamical system to consider is formed by the two Eqs. (A7) and (A10). In this case of the two mathematical attractors (x eq , y eq ): we must exclude the second because it delivers a zero Hubble function at equilibrium. For the physically relevant equilibrium we get: (H eq , z eq , q eq , w eq ) From the values of such physical quantities we must restrict − 1 5 ≤ b ≤ 0, where the case b = 0 corresponds to a single fluid pure dark energy pictured as a cosmological constant universe. Moreover, the relation The Jacobian matrix associated to the dynamical system when specified to such equilibrium point reads as: implying Thus, taking into account the range for the free parameters of the model, we get DetJ > 0, which means that our equilibrium point is a a node (spirals are not possible because (TrJ) 2 4 − DetJ > 0), and we can understand that An example of the phase portrait for such a system is displayed in Fig. 1, panel b, for the numerical choice of the free parameters as b = −0.1, γ = 0.5, β = −1.2, and α = 1.5; the equilibrium point is denoted with a circle.

Dieterici: interaction term Q 1
The dynamical system to consider is formed by the two Eqs. (A7) and (A11). In this case we get three solutions (x eq , y eq ) which read as: The third must be excluded because it implies H eq = 0. Moreover: for the first and second attractors respectively. We note that the second equilibrium point does not fulfill the physical requirements. In fact For the first attractor we must restrict the values of the model parameters as − 1 5 ≤ b ≤ 0, where the case b = 0 corresponds to a single fluid pure dark energy as a cosmological constant universe, and β < 2(b−1)e −2 . The Jacobian matrix for the dynamical system specified to the first attractor is Fig. 1 An example of the phase portrait for the dynamical system accounting for the interaction term (14), and equation of state (24) in a, (25) in b, (26) in c, for the numerical choice of the free parameters as b = −0.1, γ = 0.5, β = −1.2, and α = 1.5; the equilibrium point is denoted with a circle. We note that the evolution of the system for case (c) stops when the dynamical equations are not anylonger defined, showing the part of the phase plane which can be effectively explored by this system which implies Thus, taking into account the range for the free parameters of the model, and noticing that we get DetJ < 0, which means that our equilibrium point is a saddle. An example of the phase portrait for such a system is displayed in Fig. 1c, for the numerical choice of the free parameters as b = −0.1, γ = 0.5, β = −1.2, and α = 1.5; the equilibrium point is denoted with a circle. We note that in this case the evolution of the system stops when the dynamical equations lose meaning, that is when the argument of the Lambert W function becomes smaller than −1/e.

Redlich-Kwong: interaction term Q 2
The dynamical system to consider is formed by the two Eqs.
(A14) and (A16). In this case we get five solutions (x eq , y eq ) which read as: We notice that the third and fourth equilibria deliver an illdefined H eq , while the second and fifth implies a zero H eq . Thus we can consider only the first point as an appropriate attractor for cosmology, for which we get: (H eq , z eq , q eq , w eq ) which corresponds to a single fluid universe dominated by dark energy in form of a cosmological constant. The model parameters must be further restricted requiring β < −1.
The Jacobian matrix associated to the dynamical system when specified to such equilibrium point reads as: implying Thus, taking into account the range for the free parameters of the model, we get DetJ > 0, which means that our equilibrium point is a node (spirals are not possible because (TrJ) 2 4 − DetJ > 0), and we can understand that it is stable because TrJ < 0. An example of the phase portrait for such a system is displayed in Fig. 2a, for the numerical choice of the free parameters as b = −0.1, γ = 0.5, β = −1.5, and α = 1.5; the equilibrium point is denoted with a circle.
3.5 Modified Berthelot: interaction term Q 2 The dynamical system to consider is formed by the two Eqs. (A14) and (A17). In this case we get four attractors (x eq , y eq ) which correspond to: The first equilibrium point is not physical because it delivers an ill-defined H eq , and we must exclude also the second and fourth solutions because they imply a zero H eq . For the third solution, which is the only one relevant for cosmology, we get: which corresponds to a single fluid universe dominated by dark energy in form of a cosmological constant. The model parameters must be further restricted requiring β < −1. The Jacobian matrix associated to the dynamical system when specified to such equilibrium point reads as: which implies Thus, taking into account the range for the free parameters of the model, we get DetJ > 0, which means that our equilibrium point is a node (spirals are not possible because (TrJ) 2 4 − DetJ > 0), and we can understand that it is stable because TrJ < 0. An example of the phase portrait for such a system is displayed in Fig. 2b, for the numerical choice of the free parameters as b = −0.1, γ = 0.5, β = −1.5, and α = 1.5; the equilibrium point is denoted with a circle.
3.6 Dieterici: interaction term Q 2 The dynamical system to consider is formed by the two Eqs. (A14) and (A18). In this case we get five mathematical equilibria (x eq , y eq ): of which we must exclude the second and fourth because they deliver a zero H eq . For the three cosmologically relevant solutions we get: (H eq , z eq , q eq , w eq ) respectively. For the first attractor, which pictures a darkmatter-free universe, we must restrict − 3e 5 ≤ β < − e 6 ; the particular case β = − e 2 corresponds to a cosmological constant. The second attractor pictures as well a dark-matter-free universe in which dark energy is pictured with a cosmological constant, and for it we must restrict β ≥ − e 2 for having a well-defined Hubble function. The third equilibrium point where the lower limit comes from requiring a negative deceleration parameter and a negative parameter w eq for the dark energy equation of state, while the upper limit comes from requiring a positive dark matter abundance parameter: moreover a negative deceleration parameter requires b < 2γ −1 3 ; then all the other physical requirement are automatically satisfied.
We note that the second attractor constitutes a limit of the first attractor for the particular choice β = − e 2 , which consequently delivers a bifurcation. A bifurcation between the first and third attractor is possible when the condition 2β(1+γ ) (b−γ )e = 1 holds. The three attractors reduce to a unique equilibrium point for the choice (β = −e/2, b = −1). The Jacobian matrices are ill-defined in the first and third equilibria, while in the second it reads: An example of the phase portrait for the dynamical system accounting for the interaction term (15), and equation of state (24) in a, (25) in b, for the numerical choice of the free parameters as b = −0.1, γ = 0.5, β = −1.5, and α = 1.5; the equilibrium point is denoted with a circle. c-g Display the phase space for the system for the same interaction term, but with equation of state (26) for the choices of the free parameters as (γ = 0.1, b = 0.5, β = −0.5), respectively showing a rich behavior and possible bifuractions among the various attractors. Note that in this latter case the evolution equation can be studied as long as the Lambert W function is well defined, thus showing if the system can effectively reach these attractors which implies Being DetJ < 0, this equilibrium point is a saddle. Panels c-g of Fig. 2 display the phase space for this cosmological system for the choices of the free parameters as ( respectively showing a rich behavior and possible bifuractions among the various attractors. Note that in this latter case the evolution equations can be studied as long as the Lambert W function is well defined, thus showing if the system can effectively reach these attractors.
3.7 Redlich-Kwong: interaction term Q 3 The dynamical system to consider is formed by the two Eqs.
(A22) and (A24). In this case we get five mathematical equilibria (x eq , y eq ): of which we can keep only the first since the others deliver an ill-defined (second and third), or zero (fourth and fifth) H eq . For the cosmologically relevant solution we get: The restrictions to apply to the free parameters of the model are − 1 6 ≤ b ≤ 0 (from the values of x eq , z eq , and w eq ), , from the value of H eq . The particular case b = 0 corresponds to the case of a dark-matter-free universe in which dark energy is pictured as a cosmological constant. The Jacobian matrix for the dynamical system specified to this attractor is which implies Since DetJ > 0 this equilibrium point cannot be a saddle, but only a node or a spiral. We note that the value of the parameter α does not affect this analysis, whose outcome consequently depends on the relation between γ and the other two remaining free parameters. Figure 3 depicts the behavior of TrJ, and of (TrJ) 2 4 − DetJ as functions of γ and β for the choice b = −0.1 with the latter being on the top for γ = −1. An example of the phase portrait for such a system is displayed in Fig. 6a- Fig. 3 This figure shows the behavior of TrJ, and of (TrJ) 2 4 −DetJ given by Eqs. (72) and (73), as functions of γ and β for the choice b = −0.1 with the latter being on the top for γ = −1 γ = 1.0, β = −1.12, α = 1.5) respectively; the equilibrium points are denoted with a circle, and are respectively a stable spiral, a stable node, and a unstable node for the previously mentioned choices of the free parameters.
3.8 Modified Berthelot: interaction term Q 3 The dynamical system to consider is formed by the two Eqs. (A22) and (A25). In this case we get three mathematical equilibria (x eq , y eq ): of which we can keep only the first since the others deliver a zero H eq . For the cosmologically relevant solution we get: (H eq , z eq , q eq , w eq ) The restrictions to apply to the free parameters of the model are − 1 6 ≤ b ≤ 0 (from the values of x eq , z eq , and w eq ), and β < − 1 1+b from the value of H eq . The particular case b = 0 corresponds to the case of a dark-matter-free universe in which dark energy is pictured as a cosmological constant. The Jacobian matrix for the dynamical system specified to this attractor is which implies Since DetJ > 0 this equilibrium point cannot be a saddle, but only a node or a spiral. We note that the value of the parameter α does not affect this analysis, whose outcome consequently depends on the relation between γ and the other two remaining free parameters. Figure 4 depicts the behavior of TrJ, and of (TrJ) 2 4 − DetJ as functions of γ and β for the choice b = −0.1 with the latter being on the top for γ = −1. An example of the phase portrait for such a system is displayed in Fig. 6d- 12, α = 1.5) respectively; the equilibrium points are denoted with a circle, and are respectively a stable spiral, a stable node, and a unstable node for the previously mentioned choices of the free parameters.
3.9 Dieterici: interaction term Q 3 The dynamical system to consider is formed by the two Eqs. (A22) and (A26). In this case we get five mathematical equilibria (x eq , y eq ): of which we can keep only the first since the others deliver a zero or an ill-defined H eq . For the cosmologically relevant solution we get: (H eq , z eq , q eq , w eq ) The restrictions to apply to the free parameters of the model are − 1 6 ≤ b ≤ 0 (from the values of x eq , z eq , and w eq ), and − e 2(b+1) < β < − 2 e 2 (1+b) from the value of H eq . The Jacobian matrix for the dynamical system specified to this attractor is which implies Figure 5 shows that DetJ < 0 implying that our equilibrium point is a saddle. An example of the phase portrait for such a system is displayed in Fig. 6g, for the numerical choices 3.10 Summary Table 1 summarizes the nine cases given by the combinations of equation of state for the dark energy and the type of interaction between dark matter and dark energy. Each cell reports the values assumed by the physical quantities (x eq , y eq , H eq , z eq , q eq , w eq ) at equilibrium, the type of the attractor, and the restrictions on the numerical values of the free parameters of the model following from the requirements 0 ≤ x eq ≤ 1, 0 ≤ z eq ≤ 1, − 6 5 = −1.2 ≤ w eq = y eq x eq < 0, q eq = 1 2 (1 + 3y eq ) < 0, well defined real and non-zero H eq (for satisfying the Friedmann equation), and the general requirements −1 ≤ b ≤ 1, −1 ≤ γ ≤ 1, α > 0, β < 0.
To summarize, our analysis about linearly interacting dark matter -nonideal dark energy fluids in the linearized regime of the dynamical evolution allows us to establish some recurrent physical patterns among the different modelings for dark energy.

The strength of interaction (that is the deviation from an
ideal gas behavior p ∝ ρ) inside the dark energy fluid quantified by the parameter α does not affect the stability nature of the late-time attractors. 2. Similarly, α enters only the value of the Hubble function H eq , and not the other physical quantities at the equilibrium like the deceleration parameter, the matter abundance, and the effective parameter of the dark energy equation of state. 3. In the light of the two remarks above, the numerical value of α is not restricted by any late-time cosmological requirement.
4. When an interacting term linear in the dark matter is assumed, the final states are the same as the limiting b → 0 case when the interaction is linear in the dark energy (i.e. that in each table the first row reduces to the second for b = 0). 5. Fixing the type of dark energy -dark matter interaction and varying the equation of state for dark energy brings to final states which differ only for the value of the Hubble constant, that is the age of the universe, with the same other cosmological quantities. 6. The limit of a diverging α → ∞ delivers a zero cosmological constant for all the combinations eos -interaction type. 7. Only the Dieterici equation of state supports bifurcations between the late-time attractors. 8. The late-time attractors that we have found supporting a negative deceleration parameter correspond to a de Sitter universe. 9. Our analysis has identified as well some equilibria characterized by the vanishing of the Hubble function. These points correspond to a Minkowski universe. In fact, from (4) we understand that a zero Hubble function necesserely implies absence of energy density for both dark energy and dark matter. Then, from the equations of state considered in this paper (18)- (20) we obtain that a zero energy density for dark energy implies a zero dark energy pressure. Therefore, from the conservation Eqs. (6)- (7), taking into account our modeling of the interaction term (14)-(16), we can check explicitly that this equilibrium point is a Minkowski universe. Interestingly, this would not be the case when dark energy is modeled according to the Modified Chaplygin Gas [62] with p = A + Bρ α , or adopting the Generalized Chaplygin Gas [14] with p = βρ α (in this latter case with negative α). In fact, a zero dark energy density delivers a non-zero dark energy pressure, which in turn provides a non-zeroḢ through (5), and a non-zero Einstein curvature tensor contrary to the case of a Minkowski universe. The dynamical variables (H ,Ḣ ) seem more suited for investigating the stability of this latter equilibrium point rather than (x, y) [63] (Table 2).

Singularities
The type of singularities that can arise at time t s in a cosmological model according to [64,65] can be classified into five cases: This is a weak singularity [75,76]; where a s , ρ s and p s are some finite constant values. First of all we will explain why certain singularities cannot arise simply looking at the equations of state for dark energy we have adopted. Then, we will classify the remaining types of singularity which can arise in our model looking at the leading terms in the Friedmann evolution equations, essentially following the same line of thinking of [77,78] (Table 3).

Type V singularity
The equations of state for the dark energy modeling we considered, namely the Redlich-Kwong (18), the modified Berthelot (19), and the Dieterici (20) exhibit the following behaviors: (e) (f) (g) (b) (c) (d) Fig. 6 An example of the phase portrait for the dynamical system accounting for the interaction term (16), and equation of state (24) in a-c, (25) in d-f, and (26) Stable spiral, or stable node, or unstable node respectively. Therefore, all the models considered in this papar are type V-singularity-free because they cannot admit at the same time a zero energy density and a diverging effective equation of state parameter. We stress that this result relies only on the particular equations of state we have considered, and not on the interactions terms. Moreover, this is a qualitatively similar behavior than the one that can be obtained for a cosmological modeling which adopts the Chaplygin gas p = − A ρ n as dark energy for which the pressure cannot vanish in the correspondence of a vanishing energy density.

Type I and type III singularity
From (84)-(86) we understand that a type I singularity cannot happen when we model dark energy in terms of the modified Berthelot and of the Dieterici equations of state because it is not possible to have a diverging pressure in correspondence of a diverging energy density, which is qualitatively the same behavior of the Chaplygin gas for positive n. Again, we stress that this observation fully relies only on the particular functional form of the equations of state considered and not on the choice of the interaction terms. The same conclusion holds for the type III singularity.
For the Redlich-Kwong scenario, in proximity of a possible type I or III singularity (if it arises) we can approximate Table 3 This table summarizes the values of the parameters which allows for a type I, or II, or III, or IV, or V singularity, respectively. The simbol # means that that type of singularity cannot occur within our range of interest for the parameters Interaction type EoS

Redlich-Kwong: interaction term Q 1
Using (87) we can approximate the evolution Eqs. (4), and (A5)-(A6) nearby the singularity as The second equation can be integrated as where C 1 is a constant of integration. Therefore, a type I singularity may be possible if where we used the physical restriction −1 ≤ γ ≤ 1, while a type III singularity may be possible if We note that in the absence of the interaction term, i.e. for b = 0 = γ , a type I singularity may still be possible for a phantom dark energy with β < −1, while a type III singularity cannot be realized. Then the evolution eq. for the dark matter density provides where C 2 is another constant of integration. C 1 and C 2 can be determined by imposing ρ de (a 0 ) = ρ de 0 and ρ dm (a 0 ) = ρ dm 0 respectively. In the proximity of a type I singularity we can neglect the first term in the evolution of the energy density for the dark matter respect to the contribution of the dark energy, and write the evolution of the scale factor as where C 3 is a constant of integration which can be determined requiring a(t 0 ) = a 0 . Therefore, we can have a type I singularity, which requires a diverging scale factor at a finite time t = t s if the denominator of the fraction within the square bracket is vanishing and simultaneously the exponent is positive, or if the numerator of the fraction within the square bracket is vanishing and simultaneously the exponent is negative. In the former case we get the following requirement which must be exluded in the light of (92), while in the latter case Therefore, a type I singularity can still arise in the latter case for a phantom fluid with β < −1 even without interactions. We note that the strength of interactions within the dark energy fluid quantified by the parameter α has not played any role in this analysis. If a type III singularity arises, we can write the evolution for the scale factor as and we can neglect the second term in its rhs in light of (93) and using that the exponential is dominating over the polynomial. Therefore Observing that we understand that we can get a type III singularity when (93) applies. Again the parameter α has not played any role because in the high-energy regime the Redlich-Kwong fluid behaves effectively like an ideal fluid.

Redlich-Kwong: interaction term Q 2
The evolution for the energy density of dark matter in this case reads as where C 1 is some constant of the integration. The evolution for dark energy is where C 2 is another constant of integration. A type I singularity may happen if The former condition cannot happen for −1 γ 1 and −1 b 1. The latter implies − 6 5 β < −1, that is that dark energy must be a quintessence fluid. Moreover, we ignore the case in which the denominator is zero because the exponential term dominates over the polynomial term. Nearby the singularity we can approximatė which can be integrated for the scale factor as where C 3 is another constant of integration. The scale factor in this case is always regular at a finite time for − 6 5 β < −1. Therefore a type I singularity cannot occur.
A type III singularity can happen when The first condition implies that β= b−γ 1+γ , and therefore requires b < γ for being realized. The latter condition cannot be realized since −3(b+1) γ +1 < 0. Thus, nearby the singularity we can approximatė which can be integrated for the scale factor as where C 3 is a constant of integration. Therefore a type III singularity, characterized by a(t) → a s , cannot be realized because 1+γ 3(1+b) > 0, as already discussed, and taking into the account vanishing factor in the denominator which follows from (106).

Redlich-Kwong: interaction term Q 3
In this limit, we can approximate the energy conservation Eqs. (A20) and (A21) witḣ which deliver the following time evolution for the dark energy density: where C 1 and C 2 are two arbitrary constants of integration.
Since neither u nor v can diverge, a type III singularity, characterized by lim a→a s ρ de = ∞, cannot occur in this framework. For investigating the possible occurrence of a type I singularity, characterized by lim a(t)→∞ ρ(t) = ∞, we observe that we must require u + v > 0. Since u + v will dominate over u − v, it is enough to study the sign of the former. Explicitly we must require which after some algebraic manipulations can be recast as In particular, this requires b =-1, and more im general the condition − 6 5 β < 1 will impose − 1 where in the last step we imposed β < 0. Applying the parameterization of [77] to the energy density rather than to the Hubble function we write where in the last step we kept only the leading term. A diverging energy density in the neighborhood of t = t s requires n < 0. By inverting for the scale factor we get and since u + v > 0 (under the conditions we have derived previously), we can get a(t) → ∞ for t → t s . Therefore, a type I singularity can occur when the following conditions are met simultaneously: β < − 1 1+b , − 1 6 β 1 and b < γ .

Type II singularity
We may have a type II singularity in correspondence of the values of the energy densities which make the denominator of the equations of state vanish. We observe that: in which we used α > 0 as from [44]. Therefore, a type II singularity cannot arise in the models involving the Redlich-Kwong and the modified Berthelot equations of state, but it may arise in a framework supported by the Dieterici equation of state. Following [77], nearby a singularity, we can parameterize the Hubble function as In a type II singularity H (t) is finite because ρ(t) is finite and we use (4). Therefore, n 0. Integratinġ we obtain the evolution of the scale factor as where C 1 is some constant of the integration and t s is the time at the singularity. For n 0, a(t) is finite about t = t s , as we need in a type II singularity. Therefore, we just need to check explicitly whether we have a finite energy density nearby the singularity for each type of interaction. First of all, by expanding the Dieterici EOS about ρ = 2 α we get in which the leading term is the first one.

Dieterici: interaction term Q 1
Keeping only the leading term in the pressure and substituting it into the energy conservation equations we geṫ From (126) we get where the double sign corresponds to approaching ρ s from the left and from the right respectively. Therefore, for avoiding a divergence in the energy density we must require that α = 0 and γ = 1. In this type of singularity α is playing an important role (contrary to other types of singularity). In fact for an ideal fluid with p βρ, we would get p → ∞, and a type II singularity would not be realized. This result is compatible with ρ s = 2 α ⇒ α = 0 (for having a finite energy density), but we get one more piece of information such that γ = 1. By plugging this ρ de intoρ dm , and integrating over the time, we understand that these two conditions α = 0 and γ = 1 are enough for guaranteeing a finite ρ dm as well.

Dieterici: interaction term Q 2
By repeating the same steps, the dark energy density, when the second interaction is considered, can be obtained as where the double sign must be interpreted in the same way in the previous paragraph. Similarly, as in previous case we must require α = 0. The energy density for dark matter is To have a finite energy density we must required that γ = −1.

Dieterici: interaction term Q 3
For the case of the third interaction term we obtain where the double sign must be interpreted in the same way in the previous paragraph. Therefore, we must require α = 0, that is that the dark energy must be pictured as a nonideal fluid. The dark matter energy density for this case reads as which can be finite when α =0 and γ = −1, (taking into account that β = 0 for dark energy fluid).

Type IV singularity
To find out if there is such a singularity or not, it is only necessary to check whether the second derivative of Hubble function is ill-defined, since, in this case, all the higher order derivatives will be ill-defined automatically. From the Friedman equations we can get the first derivative for Hubble function aṡ Taking the second derivative of Hubble function and using energy conservation equations geẗ In Type IV singularity this term is regular by definition.
If a type IV singularity occurs, we must have a divergent pressure for dark energy | p de | → ∞. By looking carefully at the second termṗ de which reads aṡ We can see that it can be divergent when ∂ p ∂ρ → ∞ orρ → ∞. Note that theρ depends on the type of the interaction. For each equation of state we calculate • For Redlich-Kowng equation of state: Therefore, a type IV singularity may happen if ρ = . However in this case we get ρ < 0. So there is no type IV singularity.
• For Modified Betherlot equation of state: This quantity is always regular for α > 0. Thus a type IV singularity does not happen in this case. • For Dieterici equation of state: In this case we may have a type IV singularity if ρ=2/α. However in this case we get | p de | → ∞. Thus there is no type IV singularity. Therefore, we can have type IV singularity when → ⎧ ⎨ ⎩ γ = 1 for 1st interaction γ = −1 for 2nd interaction no possibilities for 3rd interaction.

Conclusions
The theoretical construction of a cosmological model is an interplay between the following three independent aspects: a gravitational theory, the symmetry group of the manifold, a modeling for the matter content. It is important to clarify the specific role of each three of them on driving the evolution of the Universe. For example, the effects of a Chaplygin gas can be accounted for also by an appropriate f (R) theory [79], or by an effective inhomogeneous cosmological model [80], and similarly the effects of dark matter can be mimicked by a non-local gravitational theory [81]. These claims triggered a line of research trying to break the possible degeneracy between physically and conceptually different frameworks [82][83][84]. In particular, we are interested in understanding how this interplay affects two basic aspects of the cosmological model under investigation: the existence of late-time equilibrium solutions characterized by a negative deceleration parameter, and the types of singularities which can arise. As the starting point, in this manuscript, we explicitly checked the existence of accelerating attractors in some models based on general relativity, the flat Friedmann metric, and some non-ideal dark energy fluids (already proposed in the literature) linearly interacting with dark matter. In our models a phase transition between decelerating and accelerating universes is possible (see Appendix B). We showed that the uniqueness of the accelerating late-time state depends on the modeling of the dark energy because in one case bifurcations are allowed. In our class of models, all the late-time attractors which support a negative deceleration parameter correspond to a de Sitter universe. Then, we explained that certain types of singularity cannot arise just looking at the specific equations of state we considered, without the need of any information about the theory of gravity and the metric element we adopted. Moreover, we showed that the deviations from an ideal fluid for the dark energy do not affect the existence of other types of singularity. Similarly, we showed that this quantity affects the value of the Hubble function at equilibrium (which is the age of the universe), but not of the other cosmological parameters. Furthermore the request of a negative deceleration parameter has allowed us to establish a set of constraints among the free parameters of the models under investigation. In the forthcoming works we consider appropriate to investigate whether the classification of the other singularities we exhibited remains true when we consider dark energy interacting with scalar fields carrying a non-zero pressure, or when we introduce high-energy (branes) or long-distance corrections to general relativity for extending [85], possibly accounting for the role of a massive graviton [86,87], bouncing cosmology [88,89], in quintom cosmology [90], in mimetic gravity [91], or when we abandon the Copernican principle (e.g. by adopting a Bianchi or an inhomogeneous model). Also, we will check whether our claim that the non-ideality of the fluid does not affect the values of the cosmological parameters at equilibrium, but for the age of the Universe, should be re-examined. Our results about the classification of the singularities depending on certain relationships between the parameters of the model are also important in light of the studies of their quantum corrections (accounted for by the Wheeler-de Wit equation), as done for example in the case of dark energy models based on the non-ideal Shan-Chen fluid in which the cases of singularities characterized by a finite or diverging energy density where compared and contrasted [92]. and in general they depend on the values of the free parameters α 1 ,…,α n . It is possible that, when more than one equilibrium is found, a particular choice of the values of the free parameters can bring to the identification of two or more equilibria into a single point: in this case we say that we have a bifurcation. In general, the system is expected to tend to a certain equilibrium point at late times (thus we can speak also of attractors), whose value depend on the initial conditions chosen (if more than one equilibrium is allowed). Hartman-Grobman theorem provides the classification of the stability of equilibrium points in the linearized regime [96]. Let be the Jacobian matrix of the system evaluated at the equilibrium point of interest. Then, according to the nature of its eigenvalues, or equivalently of its determinant DetJ and of its trace TrJ, the stability of the equilibrium point can be classified as follows: • DetJ < 0 ⇔ the eigenvalues are real and of opposite sign ⇔ the phase portrait is a saddle (which is always unstable); • If 0 < DetJ < (TrJ) 2 4 ⇔ the eigenvalues are real, distinct, and of the same sign ⇔ the phase portrait is a node, stable if TrJ < 0, unstable if TrJ > 0; • If 0 < (TrJ) 2 4 < DetJ ⇔ the eigenvalues are neither real nor purely imaginary ⇔ and the phase portrait is a spiral, stable if TrJ < 0, unstable if TrJ > 0.
Essentially this classification describes whether the system returns back to the equilibrium point or depart from it after a small perturbation.
In the following we derive the dynamical autonomous system for each interaction model that we consider: 5.1 Q 1 = 3Hbρ de + γρ de The choice of the interaction term (14) implies the following laws for the conservation of the dark energy and dark matter: Consequently the dynamical system to consider becomes: The evolution equation for y in our explicit cases (24), (25), and (26) reads as respectively.

Q
The choice of the interaction term (15) implies the following laws for the conservation of the dark energy and dark matter: Consequently the dynamical system to consider becomes: where the latter in our explicit cases (24), (25), and (26), is given by: respectively.

Q
The choice of the interaction term (16) implies the following laws for the conservation of the dark energy and dark matter: Consequently the dynamical system to consider becomes: The evolution equation for y in our explicit cases (24), (25), and (26) reads as respectively.

Appendix B: Phase transitions
In this appendix we investigate the phase transitions between early and late-time epochs of the universe in each cosmological model which we have considered in our analysis. This can be done by looking at the evolution of the deceleration parameter respect to the e-folding number. The deceleration parameter is in which the evolution of y is computed from the solution of the dynamical system. In Fig. 7 we display the results for our nine cases.