We revise the cosmological interaction between dark energy and dark matter. More precisely we focus on models that support compartmentalization or coexistence in the dark sector of the universe. Within the framework of a homogeneous and isotropic spatially flat FriedmannLematreRobertsonWalker geometry we analyse the asymptotic behaviour of the physical parameters for two interacting models where dark energy and dark matter have constant equations of state parameters in the presence of dark radiation when dark energy is described by a quintessence scalar field. For each model we determine the asymptotic solutions and attempt to understand how the interaction affects the cosmological evolution and history.
['General Relativity and Quantum Cosmology']
Compartmentalization and Coexistence in the Dark Sector of the Universe
{'Compartmentalization and Coexistence in the Dark Sector of the Universe': "1, 2, 3, 4, ∗ 1, † 4, ‡ 2, 5  \nCentre for Space Research, North-West University, Potchefstroom 2520, South Africa  \nAndronikos Paliathanasis, Kevin Duffy, Amlan Halder, and Amare Abebe 1 Institute of Systems Science, Durban University of Technology, Durban 4000, South Africa 2  \n3 Departamento de Matem'aticas, Universidad Cat'olica del Norte,  \nAvda. Angamos 0610, Casilla 1280 Antofagasta, Chile  \n4 School of Technology, Woxsen University, Hyderabad 502345, Telangana, India 5 National Institute for Theoretical and Computational Sciences (NITheCS), South Africa  \nWe revise the cosmological interaction between dark energy and dark matter. More precisely, we focus on models that support compartmentalization or coexistence in the dark sector of the universe. Within the framework of a homogeneous and isotropic, spatially flat Friedmann-Lemaˆıtre-Robertson-Walker geometry, we analyse the asymptotic behaviour of the physical parameters for two interacting models, where dark energy and dark matter have constant equations of state parameters, in the presence of dark radiation, when dark energy is described by a quintessence scalar field. For each model, we determine the asymptotic solutions and attempt to understand how the interaction affects the cosmological evolution and history.  \nPACS numbers: 98.80.-k, 95.35.+d, 95.36.+x", '1. INTRODUCTION': 'The nature of dark energy and dark matter is a mystery. These exotic fluids do not interact with light, and we observe only the gravitational effects they are responsible for; they form the so-called dark sector of the universe. Dark matter has been introduced to explain the velocity profile of galaxies [1]. On the other hand, dark energy [2] has been  \nintroduced to explain cosmic acceleration [3, 4]. In contrast to dark matter, which exerts gravitational effects, dark energy exerts repulsive forces that lead to anti-gravitational effects. In gravitational physics, matter sources are described by the formalism of fluid dynamics.  \nOn large scales, where the universe is isotropic and homogeneous, as described by the Friedmann-Lemaˆıtre-Robertson-Walker (FLRW) geometry, dark matter is described by a pressureless fluid source, also known as dust, while dark energy is described by a perfect fluid with a negative equation of state parameter [2].  \nThe cosmological constant provides the simplest mechanism for describing dark energy [3]. Despite the fact that the cosmological constant introduces the fewest degrees of freedom and is in agreement with the majority of cosmological data, it suffers from various problems. For a recent discussion on the challenges of the cosmological constant, we refer the reader to [5]. In order to overcome the puzzles of the cosmological constant, cosmologists have dealt with the introduction of other dark energy models with a dynamical equation for the state parameter. The origin of the new dynamical variables is mainly attributed to the existence of scalar fields [6-13], or to modifications of the gravitational action integral [14-22].  \nInteraction in the dark sector of the universe [23], that is, energy transfer between dark matter and dark energy, has been proposed as a potential mechanism to explain the cosmic coincidence problem [24-27]. Indeed, the existence of this interaction leads to new dynamical behaviour of the physical variables. Because of an interaction it is possible the nature of the dark energy fluid can change. Although the interaction between dark energy and dark matter is mainly a phenomenological consideration, there are various theoretical models that predict an interaction, such as a Weyl integrable spacetime [28, 29], the Chameleon mechanism [30, 31], a conformal transformation [32, 33], Chiral scalar fields [34-37] and others. The analysis of the cosmological observations with the interacting models has shown that the existence of a nonzero energy transfer term between these two fluids can be a potential solutions for the cosmological tensions, that is, the H 0 tensions and the σ 8 tension [38-45]. There is a plethora of interacting models which have been proposed in the literature, see for instance [46-62] and references therein.  \nIn biological systems, understanding dynamic interactions and predicting future states can be quite complex due to the multitude of interacting components and variables. Compartmental models serve as a powerful tool in simplifying these complexities by segregating a biological system into interconnected, homogeneous compartments representing different  \nstates or classes of elements within the system. Compartmental models can provide detailed quantitative analyses by defining transfer rates between compartments, which is crucial for understanding system dynamics and predicting future states, especially useful in epidemiology. These models also allow for robust parameter estimation and can be adapted or generalized for various studies. Additionally, they integrate data from multiple sources, enhancing analysis accuracy and robustness across fields like ecology and epidemiology. A similar approach is used in this study, where dark matter and dark energy are considered as compartments of variable composition. An additional term in the energy conservation equations represents the energy exchange between these compartments. Analyses of the differential equations of the system are used to explore the long-term behavior of the universe, using techniques to simplify the equations.  \nIn this study, we employ dynamical system analysis to investigate the evolution of the physical parameters for some dark energy and dark matter interacting models which describe compartmentalization or coexistence in the dark sector of the universe. The analysis of phase space is a powerful method for studying the asymptotic behaviour of a given cosmological model, as well as for the reconstruction of the provided cosmological history [63-72]. The method has been widely applied in the analysis of background equations in model cosmology, as well as in the analysis of the asymptotic behaviour of cosmological perturbations [73, 74]. The structure of the paper is as follows.  \nIn Section 2, we introduce the cosmological model under our consideration, which is that of an isotropic and homogeneous spatially flat FLRW geometry. Moreover, the continuity equation for the cosmological fluid is discussed, and the concept of interaction is introduced. In Section 3, we assume that the cosmological fluid consists of dark energy and dark matter. These two fluids are assumed to have constant equation-of-state parameters, but they have a nonzero interaction term. We consider two interacting functions which can describe compartmentalization or coexistence between the two fluids in the dark sector of the universe. We employ the H -normalization approach and perform a detailed analysis of the solution space for the field equations. In particular, we calculate the stationary points and determine their stability properties, as well as the physical properties of the asymptotic solutions.  \nIn order to make our analysis more practical, in Section 4, we introduce a radiation term. The existence of radiation leads to a new behaviour in the evolution of the dynamics. Moreover, in Section 5, we replace the minimally coupled radiation fluid with that of dark  \nradiation. This is an exotic matter source which interacts with the other two fluids of the dark sector of the universe. In Section 6, we assume that the dark energy fluid is described by a scalar field minimally coupled to gravity. In this framework, dark energy has a dynamical equation of state parameter. We focus only on the case where dark matter exists and investigate the solution space. Because of the new dynamical variables introduced by the scalar field, the evolution of the physical solutions differs from those found before. Finally, in Section 7, we summarize the results and draw our conclusions.', '2. THE COSMOLOGICAL MODEL': "According to the cosmological principle, the geometry which describes the universe on large scales is the spatially flat FLRW metric with line element  \nds 2 = -dt 2 + a 2 ( t ) ( dx 2 + dy 2 + dz 2 ) . (1)  \nFunction a ( t ) describes the radius of a three-dimensional space, and the expansion is given by the Hubble function H = ˙ a a , ˙ a = da dt . The line element (1) is isotropic and homogeneous where from the Einstein tensor  \nG µν ≡ R µν -R 2 g µν , (2)  \nonly the diagonal terms survive  \nG µν = diag ( 3 H 2 , -2 ˙ H -3 H 2 , -2 ˙ H -3 H 2 , -2 ˙ H -3 H 2 ) . (3)  \nThe gravitational field equations of Einstein's General Relativity are  \nG µν = T µν , (4)  \nwhere T µν is the energy momentum tensor which attributes the matter source in the universe.  \nFor the FLRW geometry, and for the Einstein tensor (3) it is obvious that T µν inherits the properties of G µν ; consequently, the cosmological fluid is describe by an isotropic and homogeneous perfect fluid, i.e.  \nT µν = ( ρ + p ) u µ u ν + pg µν , (5)  \nwhere u µ is the comoving observer, u µ = δ µ t , ρ ( t ) is the energy density and p ( t ) is the pressure component. Hence, the field equations are  \n3 H 2 = ρ , (6)  \n2 ˙ H +3 H 2 = -p . (7)  \nThe Bianchi identity gives G µν ; ν = 0, from where we calculate the equation of motion for the cosmological fluid  \nT µν ; ν = 0 , (8)  \nor for the FLRW geometry  \n˙ ρ +3 H ( ρ + p ) = 0 . (9)  \nThe effective equation of state parameter for the cosmological fluid is defined as w eff = p ρ . With the application of the field equations the latter parameter is expressed as  \nw eff = -1 -2 3 ˙ H H 2 . (10)", '2.1. The Dark Sector of the Universe': "Using the latest cosmological observations, the universe is dominated by two fluid sources known as dark matter and dark energy. These two fluids do not absorb or radiate light and are detected only through the gravitational phenomena related to them. In the period of galaxy formation, dark matter dominated the universe, while today dark energy dominates.  \nThe energy momentum tensor is defined as  \nT µν = T m µν + T DE µν . (11)  \nDark matter is the 'missing' matter source in galaxy formation, whose gravitational effect can explain the velocity profile for the stars in galaxies. At cosmological scales the pressure component is neglected, that is,  \nT m µν = ρ m u µ u ν . (12)  \nOn the other hand, dark energy is responsible for the late-time rapid expansion of the universe. It has a negative pressure component, which is related to anti-gravitational forces. The energy-momentum tensor associated with dark energy is  \nT DE µν = ( ρ d + p d ) u µ u ν + p d g µν , (13)  \nfrom where we can define the equation of state parameter w d = p d ρ d . In the following, we consider w d < 0. By replacing (11) in (8) we find  \n( T m µν g νκ + T DE µν g νκ ) ; κ = 0 , (14)  \nwhich, equivalently, can be expressed as  \n( T m µν g νκ ) ; κ = Q , ( T DE µν g νκ ) ; κ = -Q. (15)  \nHere, the interaction function Q ( t ) describes the energy transfer between the two fluids. Thus from (15) in terms of coordinates it follows  \n˙ ρ m +3 Hρ m = Q, (16)  \n˙ ρ d +3 H ( ρ d + p d ) = -Q. (17)  \nWhen the fluids are minimally coupled to each other, Q ( t ) = 0, and from (15) we find the two conservation equations  \n˙ ρ m +3 Hρ m = 0 , (18)  \n˙ ρ d +3 H ( ρ d + p d ) = 0 . (19)  \n̸  \nNevertheless, it is possible for the two fluids to interact and the function Q ( t ) to be nonzero, i.e. Q ( t ) = 0. In this case there exists energy transfer between the two fluids. For Q ( t ) > 0, there is mass transfer from dark energy to dark matter, while for Q ( t ) < 0 dark matter is converted into dark energy.  \nA nonzero interaction component in the dark sector of the universe is not neglected by the cosmological data. There are studies which support the existence of mass transfer [75, 76]. Although there are theoretical models where dark energy and dark matter are two different contributions in a universe of a generic fluid source.", '2.2. Cosmographic parameters': 'There are generally two main approaches to understanding the dynamical history of our universe. The first is a model-independent approach, which focuses on kinematical quantities derived directly from the space-time metric. The second approach involves fixing either the matter content of the universe or the underlying gravitational framework to explore  \ns ( t ) = 1 a d 4 a dt 4 [ 1 a da dt ]  \nits dynamics. Kinematical quantities offer greater flexibility and novelty compared to the model-dependent approach, as they rely solely on geometrical aspects, making the analysis more robust. The study of kinematical quantities in a homogeneous and isotropic universe is known as cosmography [86].  \nUnder the background of the homogeneous and isotropic universe characterized by the FLRW universe, one can define the kinematical quantities as [86, 87]  \nH ( t ) = 1 a da dt , q ( t ) = -1 a d 2 a dt 2 [ 1 a da dt ] -2 j ( t ) = 1 a d 3 a dt 3 [ 1 a da dt ] -3 , -4 ,  \nwhere H , q , j , s are respectively called the Hubble, deceleration, jerk, snap parameters and except for the Hubble parameter, the last three parameters are dimensionless. With these terminology, the expansion scale factor a ( t ) can be written in a Taylor expansion as follows  \na ( t ) = a 0 [ 1 + H 0 ( t -t 0 ) -1 2 q 2 0 ( t -t 0 ) 2 + 1 3! j 0 ( t -t 0 ) 3 + 1 4! s 0 ( t -t 0 ) 4 + O [( t -t 0 ) 5 ] ] ,  \nin which any sub-index attached to any quantity refers to its present value. Parameters q, j and s in terms of the Hubble functions are defined  \nq = -1 -˙ H H 2 (20)  \nj = H H 3 -3 q -2 (21) ...  \ns = H H 4 +4 j +3 q ( q +4) + 6 , (22)  \nequivalently, the jerk and snap parameters are expressed in terms of the deceleration parameter as  \nj = q +2 q 2 -˙ q H , (23)  \ns = -q (1 + 2 q ) (5 + 3 q ) + (5 + 6 q ) ˙ q H -q H 2 (24)  \n,  \nThus, from the evolution of q, j, s , one can understand the expansion of the universe (accelerating or decelerating), the rate of acceleration ( j ) and its next derivative. Although jerk and snap parameters are asymptotic near today in the following we derive their dynamical evolution for the cosmological models of our consideration.', '3. DARK MATTER - DARK ENERGY': "In this piece of study, we follow a phenomenological approach and consider an interacting function such that the two fluids to be described by the simplest model of 'competitive fluids', a Lotka-Volterra type system, or a compartmental model such as the SIR model.", '3.1. The Compartmental Interaction': "From the interacting function [77, 78]  \nQ A ( t ) = Q 0 ( t ) ρ m ρ d , (25)  \nintroduced in [62], equations (15) become  \n˙ ρ m +3 Hρ m = Q 0 ( t ) ρ m ρ d , (26)  \n˙ ρ d +3 H ( ρ d + p d ) = -Q 0 ( t ) ρ m ρ d , (27)  \nwhere Q 0 ( t ) has dimensions of H -1 . The latter system can be expressed as  \n˙ ρ m +3 Hρ m = Q 0 ( t ) ρ m ρ d , (28)  \n˙ ρ d +3(1 + w d ) Hρ d = -Q 0 ( t ) ρ m ρ d . (29)  \nThis interaction model does not suffer from early instabilities [78].  \nAlthough the aforementioned system has the form of the Lotka-Volterra system with non-constant coefficients, the Hubble function H ( t ) is related to the energy densities ρ m , ρ d , from the Friedmann equation (6), that is, 3 H 2 = ρ m + ρ d . Consequently, the dynamical system (28), (29) is a type of compartmental model. These models are mainly applied in epidemiology and are known as SIR models.  \nWe introduce Q 0 ( t ) = α H , where α is a constant, and the dimensionless variables  \nΩ m = ρ m 3 H 2 , Ω DE = ρ d 3 H 2 , τ = ln a. (30)  \nRecall that the dimensionless variables have the constraint { Ω m , Ω DE } ∈ [0 , 1].  \nHence, equations (28), (29) read  \nd Ω m dτ = 3( α + w d ) Ω m Ω DE , (31)  \nd Ω DE dτ = 3Ω DE ( w d (Ω DE -1) -α Ω m ) . (32)  \nWe remark that in terms of the dimensionless variables, where the Hubble function has been eliminated by the dynamical system, it is clear that this model is a compartmental model.  \nFurthermore, from the Friedmann's equation (6) it follows  \nΩ m +Ω DE = 1 . (33)  \nWith the use of the latter constraint we can reduce by one the dimension of the dynamical system (31), (32).  \nHence the final equation is  \nd Ω DE dτ = 3( α + w d ) Ω DE (Ω DE -1) . (34)  \n̸  \na. Stationary points We assume that w d is a constant, and then we calculate the stationary points. For α + w d = 0, there exist two stationary points, point P 1 with Ω DE ( P 1 ) = 0 and P 2 with Ω DE ( P 2 ) = 1.  \nThe point P 1 describes an asymptotic solution where dark matter dominates in the universe, and the effective cosmological fluid has the equation of state parameter w eff ( P 1 ) = 0. The stationary point is an attractor for ( α + w d ) > 0.  \nOn the other hand, point P 2 corresponds to an asymptotic solution where the dark energy dominates, that is, w eff ( P 2 ) = w d . Last but not least, P 2 is an attractor when ( α + w d ) < 0.", '3.1.2. Analytic solution': 'The nonlinear differential equation (34) can be solved analytically. Indeed, the solution in terms of the scale factor τ = ln a is given as follows:  \nΩ DE ( a ) = [ 1 + Ω 0 a 3( α + w d ) ] -1 , (35)  \nfrom which we construct the Hubble function  \nH ( a ) = ˆ H 0 [ 1 + Ω 0 a 3( α + w d ) ] w d 2 ( w d + α ) [ Ω 0 a 3( α + w d ) ( α + w d ) ] -1+ w d 2 ( α + w d ) . (36)  \nHere Ω 0 and ˆ H 0 are integration constants. If w d = -1, the latter analytic expression for the scale factor becomes  \nH ( a ) = ˆ H 0 ( 1 + Ω 0 a 3( α -1) ) -1 2( α -1) . (37)  \nFor α = 0, we recover the ΛCDM universe, that is,  \nH ( a ) = H 0 ( Ω 0 Λ +Ω 0 m a -3 ) 1 2 , (38)  \nwhere we have replaced ˆ H 0 = H 0 (Ω 0 Λ ) 1 2 , and Ω m = Ω 0 Λ Ω 0 .', '3.2. The Lotka-Volterra Interaction': 'The Lotka-Volterra model describes coexistence. In order to recover such behaviour between the dark energy and the dark matter we introduce the interaction [80]  \nQ B ( t ) = Q 0 ( t ) ρ m ρ d + Q 1 ( t ) ρ d , (39)  \nwhere Q 0 ( t ) has the dimensions of H -1 and Q 1 ( t ) has the dimensions of H .  \nFor this interaction the continuous equations for the two fluids reads  \n˙ ρ m +3 Hρ m = Q 0 ( t ) ρ m ρ d + Q 1 ( t ) ρ d . (40)  \n˙ ρ d +3(1 + w d ) Hρ d = -Q 0 ( t ) ρ m ρ d -Q 1 ( t ) ρ d . (41)  \nWe proceed with an asymptotic analysis of the field equations for this interacting model, where we assume that Q 0 ( t ) = α H and Q 1 ( t ) = βH , with α , β constant parameters.  \n3.2.1. Asymptotic analysis  \nWe apply the same algorithm as before and introduce dimensionless variables (30). The field equations with dimensionless variables form the dynamical system  \nd Ω m dτ = 3( α + w d ) Ω m Ω DE + β Ω DE , (42)  \nd Ω DE dτ = 3Ω DE ( w d (Ω DE -1) -α Ω m ) -β Ω DE , (43)  \nwith constraint equation (33).  \nApplying the constraint (33) results in the single differential equation  \nd Ω DE dτ = (3 ( α + w d ) (Ω DE -1) -β ) Ω DE . (44)  \na. Stationary points The stationary points of equation (44) are the two points ¯ P 1 with Ω DE ( ¯ P 1 ) = 0 and ¯ P 2 with Ω DE ( ¯ P 2 ) = 1 + β 3( α + w d ) .  \nThe point ¯ P 1 describes a universe dominated by dark matter, that is, Ω m ( ¯ P 1 ) = 1. On the other hand, point ¯ P 2 describes a universe of coexistence between the two fluids. The solution is physically accepted when 0 ≤ Ω DE ( ¯ P 2 ) ≤ 1, that is, -1 ≤ β 3( α + w d ) ≤ 0. The equation of state parameter for the effective fluid in solution ¯ P 2 is derived w eff ( ¯ P 2 ) = w d ( 1 + β 3( α + w d ) ) .  \nAs far as the stability of the stationary points is concerned, point ¯ P 1 is attractor when 3 ( α + w d ) + β > 0, while ¯ P 2 is an attractor when -3 ( α + w d ) < β ≤ 0.  \nWe obtain similar results for the more general model of the Lotka-Volterra family, Q ( t ) = Q 0 ( t ) ρ m ρ d + Q 1 ( t ) ρ m + Q 2 ( t ) ρ d . Specifically, for this model we find two stationary points that describe coexistence. When Q 1 ( t ), or Q 2 ( t ) becomes zero, one of the coexistence points reduced to that of P 1 , or P 2 , for the studied before for the compartmental interaction.  \nFurthermore, when α = 0 , β = 0, the ΛCDM universe is recovered for ˆ H 0 = H 0 Ω Λ , and Ω m = Ω Λ Ω 0 .', '3.2.2. Analytic solution': 'The analytic solution of equation (44) in terms of the scale factor is  \nΩ DE ( a ) = ( 1 + β 3 ( α + w d ) ) ( 1 + Ω 0 a 3( α + w d )+ β ) -1 , (45)  \nand the analytic expression for the Hubble function is  \nH ( a ) = ˆ H 0 ( 1 + Ω 0 a 3( α + w d )+ β ) w d 2 ( α + w d ) a -3 ( 1+ w d ( 1+ β 3 ( α + w d ) )) . (46)  \nThus, for the case where w d = -1, the latter solution becomes  \nH ( a ) = ˆ H 0 ( 1 + Ω 0 a 3( α -1)+ β ) -1 2( α -1) a β ( α -1) . (47)  \nThe results in Eqs. (36) and (46) show there is a large wiggle room to address the issue of the Hubble tension (and the σ 8 tension). While one needs to do a proper constraint analysis  \non the parameter space of such interactions using latest and cosmological upcoming data and this is within the scope of our upcoming investigation - it is evident that this solution provides a richer range of possible values of H 0 than that provided by ΛCDM, and hence has the potential to reduce the Hubble tension. In fact, the interaction model presented in Eqs. (25) and (39) can be thought of as more general cases of the interaction model presented in [79] in which the authors show such an interaction model has an observational support in reducing both the Hubble and σ 8 tensions.', '4. DARK ENERGY - DARK MATTER - RADIATION': 'Although the contribution of the radiation in the total cosmological fluid can be neglected nowadays, radiation played an important role in the early stages of the universe. In order to present a more realistic model, we consider the energy tensor for the cosmological model to depend on the three fluids,  \nT µν = T m µν + T DE µν + T r µν (48)  \nwhere T r µν describes the radiation component defined as  \nT r µν = ( ρ r + p r ) u gm u ν + p r g µν , p r = 1 3 ρ r . (49)  \nThe conservation law (14) is modified as  \n( T m µν g νκ + T DE µν g νκ + T r µν g νκ ) ; κ = 0 , (50)  \nwhere since radiation does not interact with the rest of the dark components of the universe, it follows  \n( T m µν g νκ + T DE µν g νκ ) ; κ = 0 , ( T r µν g νκ ) ; κ = 0 , (51)  \nConsequently, in a FLRW background the continuous equation for the radiation fluid reads  \n˙ ρ r + 4 3 Hρ r = 0 . (52)', '4.1. The Compartmental Interaction': "Let us consider the compartmental interaction Q ( t ) as given in expression (25) with Q 0 ( t ) = α H . We introduce the new dimensionless variable Ω r = ρ r 3 H 2 which describes the energy density for the radiation fluid, and Ω r ∈ [0 , 1].  \nThe cosmological field equations in terms of the dimensionless variables (30) are  \nd Ω m = 3( α + w d ) Ω m Ω DE +Ω m Ω r ,  \ndτ (53) d Ω DE dτ = 3Ω DE ( w d (Ω DE -1) -α Ω m ) + Ω DE Ω r , (54)  \nwhere for the first Friedmann's equation we derive the constraint  \nΩ r = 1 -Ω m -Ω DE . (55)  \nWe proceed with the analysis of the stationary points for the two-dimensional system (53), (54). For each stationary point, we discuss the physical properties of the asymptotic solution and its stability properties. Indeed, the equation of state parameter for the effective fluid in terms of the dimensionless variables is  \nw eff = w d Ω DE + 1 3 Ω r . (56)  \na. Stationary points The dynamical system (53), (54) admits four stationary points A = (Ω m ( A ) , Ω DE ( A )), with coordinates; A 1 = (1 , 0); A 2 = (0 , 1); A 3 = (0 , 0) and A 4 = ( 1 -3 w d 3 α , -1 3 α ) .  \nThe stationary points A 1 and A 2 correspond to universes without a radiation component, that is, Ω r ( A 1 , 2 ) = 0. That is, points A 1 , A 2 have the same physical properties with that of P 1 and P 2 found before. We linearize the dynamical system (53), (54) around the stationary points and we calculate the eigenvalues e 1 ( A 1 ) = -1, e 2 ( A 1 ) = -3 ( w d + α ) and e 1 ( A 2 ) = -1 + 3 w d , e 2 ( A 2 ) = 3 ( w d + α ). Thus, the point A 1 is an attractor for α > 0 and -α < w d < 0, while A 2 is an attractor for { α < 0 , w d < 0 } , or { α > 0 , w d < -α } .  \nConcerning the asymptotic solution at point A 3 , it describes a universe dominated by the radiation fluid, that is, Ω r ( A 3 ) = 0. The eigenvalues around the stationary point A 3 are e 1 ( A 3 ) = 1, e 2 ( A 3 ) = 1 -3 w d , which means that point A 3 is always a source.  \nFinally, point A 4 describes a universe where all the fluid contributes in the cosmological evolution. However, because we have assumed w d < 0, it follows that the point is not physically accepted.  \nTwo-dimensional phase-space portraits for the dynamical system (53), (54) are presented in Fig. 1. The portraits are for w d = -1 and for two values of the parameter α , where  \npoint A 1 or A 2 are attractors. We remark that for w d = -1 point A 2 describes the de Sitter universe.  \nFurthermore, in Fig. 2 we present the qualitative evolution for the energy density parameters Ω m , Ω DE and Ω r , for values of the free parameters such that the future attractor is the point A 2 and describes the de Sitter solution. Finally, the qualitative evolution of the effective equation of state parameter w eff and of the cosmographic parameters are given in Fig. 3.  \n<!-- image -->  \nFIG. 1: Compartmental: Phase-space portraits for the dynamical system (53), (54) for w d = -1 and α = 1 10 (left fig.) and α = 2 (right fig.). The physical accepted stationary points are marked with red dots. In the left fig. point A 2 is an attractor, while in the right fig. the attractor is the stationary point A 2 .  \n<!-- image -->", '4.2. The Lotka-Volterra Interaction': 'In this Section, we study the asymptotic behaviour of the field equations for the interaction model (39), which leads to a modified Lotka-Volterra system.  \nFIG. 2: Compartmental: Qualitative evolution of the energy density parameters Ω m , Ω DE and Ω r , for values of the free parameters such that the future attractor of the dynamical system (53), (54) is the point A 2 and describe the de Sitter solution. Plots are for different initial conditions for the variables Ω m ( a ) , Ω DE ( a ) and Ω r ( a ).  \n<!-- image -->', '4.2.1. Asymptotic analysis': 'For interaction (39) the field equations in terms of the dimensionless variables { Ω m , Ω DE , Ω r , τ } are expressed as follows  \nd Ω m = 3( α + w d ) Ω m Ω DE + β Ω DE +Ω m Ω r ,  \ndτ (57) d Ω DE dτ = 3Ω DE ( w d (Ω DE -1) -α Ω m ) -β Ω DE +Ω DE Ω r , (58)  \nwith constraint equation (55).  \nMoreover, the equation of state parameter for the effective cosmological fluid is given by expression (56). In the following section, we investigate the stationary points of this two-dimensional system.  \na. Stationary points The dynamical system (57), (58) admits the stationary points ¯ A = ( Ω m ( ¯ A ) , Ω DE ( ¯ A )) with coordinates ¯ A 1 = (1 , 0) ; ¯ A 2 = ( -β 3( α + w d ) , 1 + β 3( α + w d ) ) ; ¯ A 3 = (0 , 0) and ¯ A 4 = ( 1 -3 w d + β α , -1 3 α + β 3 α (1 -3 w d ) ) .  \n̸  \nThe points ¯ A 1 , ¯ A 2 describe solutions without radiation, i.e. Ω r ( ¯ A 1 , 2 ) = 0. Hence, asymptotic solutions have similar properties as those described by points ¯ P 1 and ¯ P 2 respectively. Furthermore, the point ¯ A 3 describes a universe dominated by the radiation fluid. Finally, point ¯ A 4 is physically accepted for w d = 1 -β 3 , with β > 1 and α = 0. However, in this case, the point ¯ A 4 is identical to ¯ A 3 and always describes the radiation solution.  \nThe eigenvalues of the two-dimensional linearized system (57), (58) around the stationary points ¯ A 1 and ¯ A 2 are e 1 ( ¯ A 1 ) = -1, e 2 ( ¯ A 1 ) = -(3 ( w d + α ) + β ); e 1 ( ¯ A 2 ) = -1 + w d ( 3 + β α + w d ) , e 2 ( ¯ A 2 ) = 3( w d + α ) + β . Hence, point ¯ A 1 is an attractor for α > -β 3 , -3 α + β 3 < w d < 0; while point ¯ A 2 is an attractor when a < 0 and β < -3 ( α + w d ) or α > 0 with constraints { w d < -α, β < -3 ( w d + α ) } or { -α < w d < 0 , -(1+ w d )( α + w d ) w d < β < -3 ( α + w d ) } .  \nThe radiation solution described by point ¯ A 3 is always unstable, since the corresponding eigenvalues are calculated as e 1 ( ¯ A 3 ) = 1, e 2 ( ¯ A 3 ) = 1 -3 w d -β . Finally, point ¯ A 4 provides the eigenvalues e 1 ( ¯ A 4 ) = αβ + √ D ( α,β,w d ) 2 α (1 -3 w d ) , e 2 ( ¯ A 4 ) = αβ -√ D ( α,β,w d ) 2 α (1 -3 w d ) with D ( α, β, w d ) = α ( 4 (1 -3 w d ) 2 (1 -3 w d + β ) ((1 -3 w d ) ( α + w d ) -βw d ) + aβ 2 ) . It follows that when the point is physically accepted, it is always unstable.  \nAt this point, we remark that it is possible for point ¯ A 4 to exist in the dynamics for other values of the free parameters { α, w d , β } . Thus, in order to avoid nonphysical solutions in the model parameters w d and β are constrained such that point ¯ A 4 is identical to ¯ A 3 . However, in this case where w d = 1 -β 3 , β > 1, there exists a unique attractor described by point ¯ A 2 , with α < 0. In Fig. 4 we present phase-space portraits for the two-dimensional dynamical system (57), (58) for w d = 1 -β 3 , β > 1 and α < 0, where ¯ A 2 is the unique attractor.  \nThe qualitative evolution for the energy density parameters Ω m , Ω DE , Ω r and of the effective equation of state parameter w eff are presented in Figs. 5 and 6.  \nFurthermore, in Fig. 2 we present the qualitative evolution of the energy density parameters Ω m , Ω DE and Ω r for values of the free parameters such that the future attractor is the point A 2 and describes the de Sitter solution. Finally, the qualitative evolution of the effective equation of state parameter w eff and of the cosmographic parameters are given in Fig. 3.', '5. DARK ENERGY - DARK MATTER - DARK RADIATION': "Dark radiation is an exotic form of radiation which has been proposed to interact only with dark matter [81, 82]. Such an interaction scenario has recently gained attention as a potential candidate for resolving the Hubble tension [83-85]. We will, in this section, study the dynamical systems implications of the interaction, but we will return to a full-fledged tensions constraining analysis vis-'a-vis cosmological data in a subsequent work.  \nWe introduce the energy momentum tensor  \nT R µν = ( ρ R + p R ) u µ u ν + p R g µν , p R = 1 3 ρ R , (59)  \nwhich we assume describes the dark radiation. The total energy momentum of the universe is defined as  \nT µν = T m µν + T DE µν + T R µν , (60)  \nand using the Bianchi identity we can now define the following system  \n( T m µν g νκ ) ; κ = Q -V , ( T DE µν g νκ ) ; κ = -Q , ( T R µν g νκ ) ; κ = V, (61)  \nwhere the new variable V describes the interaction between dark matter and dark radiation.", '5.1. The EMR Model': 'We consider the compartmental interaction with the parameter Q defined by expression (25), and we define the parameter V as V = V 0 ( t ) ρ m .  \nHence, in the spatially flat FLRW background, equations (61) read as  \n˙ ρ m +3 Hρ m = Q 0 ( t ) ρ m ρ d -V 0 ( t ) ρ m . (62)  \n˙ ρ d +3(1 + w d ) Hρ d = -Q 0 ( t ) ρ m ρ d , (63)  \n˙ ρ R + 4 3 Hρ R = V 0 ( t ) ρ m . (64)  \nThis dynamical system defines an analogue of a SIR model for cosmology. We shall call it the dark EMR model. In an SIR model, the interaction is between an infected population and a recovered population. In this analogue, the interaction is between dark matter and dark radiation.  \nWe proceed with the analysis of the asymptotics.', '5.1.1. Asymptotic analysis': 'We introduce a new dimensionless variable Ω R = ρ R 3 H 2 , and we write the cosmological field equations in an equivalent form  \nd Ω m dτ = (3( α + w d ) Ω DE +Ω R -γ ) Ω m , (65)  \nd Ω DE dτ = 3Ω DE ( w d (Ω DE -1) -α Ω m ) + Ω DE Ω R , (66)  \nwith constraint  \nΩ R = 1 -Ω m -Ω DE . (67)  \nFinally the equation of state parameter for the total cosmological fluid is defined as  \nw eff = w d Ω DE + 1 3 Ω R . (68)  \nWe proceed with an analysis of the asymptotic behaviour.  \na. Stationary points The stationary points B = (Ω m ( B ) , Ω DE ( B )) of the twodimensional system (65), (66) are B 1 = (1 -γ, 0) , B 2 = (1 , 0), B 3 = (0 , 0) and B 4 = ( (1 -3 w d ) 3 α ( 1 -γ 3( α + w d ) ) , -1 3 α + γ (1+3 α ) 9 α ( α + w d ) ) .  \nThe asymptotic solution at point B 1 describes a universe with dark matter and dark radiation, i.e. Ω R ( B 1 ) = γ . The point is physically accepted for 0 < γ < 1. This means that there is energy transfer from dark matter to the dark radiation fluid. Points B 2 and B 3 have the same physical properties with that of points A 2 and A 3 , respectively. Indeed, the point B 2 describes a universe where only dark energy contributes to the cosmological fluid, while at the point B 3 the asymptotic solution describes a universe dominated by the radiation fluid.  \nFinally, point B 4 is an extension of the previous point A 4 , where the three fluids contribute in the universe Ω R ( B 4 ) = 1 -γ -3 w d 3 α and w eff ( B 4 ) = 1 3 ( 1 -(1 -3 w d ) 3( α + w d ) γ ) . If we require that point B 1 , which describes the matter era, exist always, that is, 0 < γ < 1, then point B 4 is physical accepted when { 0 < α < 1 3 , -α < w d < 0 , 3( α + w d ) 1+3 α ≤ γ ≤ 3 ( α + w d ) } , or α > 1 3 and { -α < w d < (1 -3 α ) 3 , 3( α + w d ) 1+3 α ≤ γ ≤ 3 ( α + w d ) } or { (1 -3 α ) 3 < w d < 0 , 3( α + w d ) 1+3 α ≤ γ ≤ 1 } For w d = -1, the latter conditions becomes { 1 < α < 4 3 , 3( α -1) 1+3 α ≤ γ ≤ 3 ( α -1) } , or { α > 4 3 , 3( α -1) 1+3 α ≤ γ ≤ 1 } .  \nAs far as the stability properties of the stationary points are concerned, point B 1 is an attractor for α > 3 w d -γ 3( γ -1) while B 2 is an attractor for α < 1 3 ( γ -3 w d ). Furthermore, the point B 3 always describes an unstable solution. The eigenvalues of the linearized system (65), (66) around the stationary point B 4 are e 1 ( B 4 ) = (3 w d -1) αγ + √ ¯ D ( α,w d ,γ ) 6 α ( α + w d ) , e 2 ( B 4 ) = (3 w d -1) αγ -√ ¯ D ( α,w d ,γ ) in which  \n6 α ( α + w d )  \n¯ D ( a, w b , γ ) = (3 w d -1) α ( 12 ( α + w d ) 2 ((2 + 3 α ) γ -3 ( α + w d )) -α (5 + 12 α ) γ 2 ) . (69)  \nWhen ¯ D ( α, w d , γ ) < 0, then (3 w d -1) αγ 6 α ( α + w d ) and the stationary point can be a stable sink, otherwise for ¯ D ( α, w d , γ ) > 0 there are values of the free parameters where the point can be an  \nattractor. Last but not least, for the value of the free variables where point B 4 is physical accepted, the unique attractor is point B 4 .  \nIn Fig. 7 we present the region space of the free parameters ( a, w b , γ ), where the asymptotic solution described by point B 4 is stable. The phase-space portrait of the dynamical system (65), (66) are given in 8. Moreover, the qualitative evolution of the physical parameters { Ω m , Ω DE , Ω R } , w eff and of the cosmographic parameters are presented in 9 and 10 respectively.', '6. QUINTESSENCE AS DARK ENERGY': 'In the previous sections, we have assumed that the dark energy has a constant equation of state parameter. However, this assumption suffers from various problems. The simplest gravitational model which describes the dynamical behaviour of the dark energy fluid is the scalar field.  \nQuintessence is a scalar field minimally coupled to gravity, where in the case of homogeneous spacetimes the corresponding fluid is described by a perfect fluid with energy density ρ ϕ and pressure component p ϕ defined as  \nρ ϕ = 1 2 ˙ ϕ 2 + V ( ϕ ) , (70)  \np ϕ = 1 2 ˙ ϕ 2 -V ( ϕ ) , (71)  \nwhere ϕ = ϕ ( t ) is the scalar field and V ( ϕ ) is the potential function which defines the scalar field mass.  \nThe equation of state parameter for the scalar field is calculated as  \nw ϕ = 1 2 ˙ ϕ 2 -V ( ϕ ) 1 2 ˙ ϕ 2 + V ( ϕ ) , (72)  \nfrom where we observe that | w ϕ | ≤ 1, and w ϕ →-1, when the scalar field potential V ( ϕ ) dominates.  \nWhen there is no interaction between the scalar field and the dark matter, equation (19) provides the equation of motion for the scalar field, which is a second-order differential equation, that is,  \n˙ ϕ ( ¨ ϕ +3 H ˙ ϕ + V ,ϕ ) = 0 . (73)  \nIn the following sections, we investigate the asymptotic dynamics for the two models under consideration, that is, the compartmental interaction model (25) and the LotkaVolterra interaction model (39). For simplicity of our calculations, we omit the radiation term.', '6.1. The Compartmental Interaction': 'Consider the interaction in (25) with Q 0 ( t ) = α ˙ ϕ √ 6 H 2 ; which results in the following system  \n˙ ρ m +3 Hρ m = α ˙ ϕ √ 6 ρ m H 2 ( 1 2 ˙ ϕ 2 + V ( ϕ ) ) , (74)  \n¨ ϕ +3 H ˙ ϕ + V ,ϕ = -α ρ m √ 6 H 2 ( 1 2 ˙ ϕ 2 + V ( ϕ ) ) . (75)  \n6.1.1. Asymptotic analysis  \nWe employ the H -normalization approach, thus, in terms of new dimensionless variables  \nx = ˙ ϕ √ 6 H , y = √ V ( ϕ ) √ 3 H , Ω m = ρ m 3 H 2 , λ = V ,ϕ V , τ = ln a (76)  \nthe field equations read as  \ndx = 1 [ 3 x ( x 2 -y 2 -1 ) -6 λy 2 -3 α ( x 2 + y 2 ) Ω m ] , (77)  \ndy dτ = 1 2 y [ 3 + √ 6 λx +3 ( x 2 -y 2 ) ] , (78)  \ndτ 2 √  \ndλ dτ = √ 6 xλ 2 (Γ ( λ ) -1) , Γ( λ ) = V ,ϕϕ V V 2 ,ϕ , (79)  \nwith the constraint equation  \nΩ m = 1 -x 2 -y 2 . (80)  \nFrom the constraint equation it follows that 1 -x 2 -y 2 ≥ 0, which means that x, y take values in the two unitary disk. Furthermore, by definition, y > 0. Recall that the energy density for the scalar field Ω ϕ is defined as Ω ϕ = x 2 + y 2 .  \nLast, but not least, in terms of the new variables the effective equation of state parameter is defined as  \nw eff ( x, y, λ ) = x 2 -y 2 . (81)  \nFor simplicity of our calculations and in order to keep the number of dimensions and the free parameters of the dynamical system low, we assume that the scalar field potential has a constant potential function V ( ϕ ) = V 0 . This means that the mass of the scalar field as given by equation (75) depends on the energy density of dark matter.  \nFurthermore, the parameter λ is always zero. Thus the dynamical evolution of the physical variables for this model is given by the two-dimensional dynamical system  \ndx dτ = 1 2 [ 3 x ( x 2 -y 2 -1 ) -3 α ( x 2 + y 2 ) ( 1 -x 2 -y 2 )] , (82)  \ndy dτ = 1 2 y [ 3 + 3 ( x 2 -y 2 )] . (83)  \na. Stationary points The stationary points of the dynamical system (82), (83) are defined on the two-dimensional plane T = ( x ( T ) , y ( T )), with coordinates T 1 = (0 , 0), T ± 2 = ( ± 1 , 0), T 3 = ( -1 α , 0 ) , T 4 = (0 , 1) and T 5 = ( x 4 ( α ) , √ 1 + x 2 4 ( α ) ) .  \nThe point T 1 with Ω m ( T 1 ) = 1 describes the matter-dominated era. The linearized system around the stationary points provides the eigenvalues e 1 ( T 1 ) = -3 2 and e 2 ( T 1 ) = 3 2 , from which we infer that the matter solutions are always unstable and T 1 is a saddle point.  \nOn the other hand, T ± 2 describe solutions where the universe is dominated by the kinetic term of the scalar field, Ω m ( T ± 2 ) = 0, and the effective cosmological fluid is that of stiff matter, i.e. w eff ( T ± 2 ) = 1. The corresponding eigenvalues are e 1 ( T ± 2 ) = 3 and e 2 ( T ± 2 ) = 3 (1 ± α ). Thus, asymptotic solutions are always unstable. T + 2 is a saddle point for α < -1, otherwise it is a source, while T -2 is a saddle point for α > 1, otherwise it is also a source.  \nThe point T 3 describes an asymptotic solution where the kinetic term of the scalar field and the dark matter contribute to the cosmological fluid, i.e. Ω m ( T 3 ) = 1 -1 α 2 . The point is physically accepted for α 2 > 1, while when α 2 = 1 it reduces to the points T ± 2 . For the cosmological fluid we calculate w eff ( T 3 ) = 1 α 2 , which means that it can not describe acceleration. Nevertheless, for α = 1 √ 3 , point T 3 describes the radiation solution. The eigenvalues are e 1 ( T 3 ) = 3 ( 1+ α 2 ) 2 α 2 , e 2 ( T 3 ) = -3 ( 1 -α 2 ) 2 α 2 , from where we infer that the point when it is physically accepted is a saddle point, otherwise it is a source.  \nThe stationary point T 4 is an attractor, e 1 ( T 4 ) = -3, e 2 ( T 4 ) = -3, and the asymptotic solution is that of the de Sitter universe, w eff ( T 4 ) = -1, Ω m ( T 4 ) = 0.  \nFinally, for point T 5 , parameter x 4 is given by the equation α = 1 x 4 ( 1+2 x 2 4 ) . However, y ( T 5 ) = √ 1 + x 2 4 ( α ) > 1, that is, the point is not physically accepted. The  \neigenvalues of the linearized system are e 1 ( T 5 ) = 3 { -1+ √ 5+4 x 2 4 [ 9+4 x 2 4 ( 5+3 x 2 4 )] } 2 ( 1+2 x 2 4 ) , e 1 ( T 5 ) = 3 { -1 -√ 5+4 x 2 4 [ 9+4 x 2 4 ( 5+3 x 2 4 )] } 2 ( 1+2 x 2 4 ) , thus T 5 is always a saddle point.  \nIn Fig. 11 we present phase-space diagrams for the dynamical system (82), (83), where we observe that T 4 is the unique attractor. Moreover, in Figs. 12 and 13 we present the qualitative evolution of the physical parameters Ω m , Ω ϕ , w eff and of the cosmographic parameters for various sets of initial conditions.  \nAt this point, it is important to mention that for α 2 ≤ 1, the physically accepted stationary points are only those of the quintessence model without interaction. Therefore, the physically accepted stationary points that describe non-zero interaction exist only when α 2 > 1.', '6.2. The Lotka-Volterra Interaction': 'For the Lotka-Volterra interaction (39) with Q 0 ( t ) = α ˙ ϕ √ 6 H 2 and Q 1 ( t ) = β ˙ ϕ 2 H , in terms of the dimensionless variables the field equations are expressed as follows  \ndx dτ = 1 2 [ 3 x ( x 2 -y 2 -1 ) -√ 6 λy 2 -3 ( x 2 + y 2 ) ( α Ω m +2 βx ) ] , (84)  \ndλ dτ = √ 6 xλ 2 (Γ ( λ ) -1) , Γ( λ ) = V ,ϕϕ V V 2 ,ϕ . (86)  \ndy dτ = 1 2 y [ 3 + √ 6 λx +3 ( x 2 -y 2 ) ] , (85)  \nwith algebraic constraint (80).  \nHence, by applying the constraint (80) and considering a constant potential function, that is, λ = 0, we end up with the two-dimensional system  \ndx dτ = 1 2 [ 3 x ( x 2 -y 2 -1 ) -√ 6 λy 2 -3 ( x 2 + y 2 ) ( α ( 1 -x 2 -y 2 ) +2 βx ) ] , (87)  \ndy dτ = 1 2 y [ 3 + √ 6 λx +3 ( x 2 -y 2 ) ] . (88)  \na. Stationary points The stationary points ¯ T = ( x ( ¯ T ) , y ( ¯ T )) of the two dimensional dynamical system (87), (88), have the following coordinates, ¯ T 1 = (0 , 0), ¯ T 2 = (0 , 1), ¯ T 3 = ( x 3 ( α, β ) , 0) and ¯ T 4 = ( x 4 ( α, β ) , √ 1 + ( x 4 ( α, β )) 2 ) . We proceed with a discussion of the physical properties for the asymptotic solutions at the stationary points.  \nThe asymptotic solution at point ¯ T 1 describes a universe dominated by the matter era, similar to point T 1 . The eigenvalues are calculated e 1 ( ¯ T 1 ) = -3 2 , e 2 ( ¯ T 2 ) = 3 2 , which means that ¯ T 1 is always a saddle point.  \nMoreover, point ¯ T 2 describes the de Sitter solution as point T 4 . The eigenvalues are e 1 ( ¯ T 2 ) = -3, e 2 ( ¯ T 2 ) = -3 (1 + β ). Hence, ¯ T 2 is an attractor when β > -1.  \nFor the stationary point ¯ T 3 , parameter x 3 ( α, β ) is given by the algebraic equation β = (1+ αx 3 ) ( x 2 3 -1 ) 2 x 2 , where for β = 0, the points T ± 2 and T 3 are recovered. The asymptotic solution at the point ¯ T 3 cannot describe acceleration, i.e., w eff ( ¯ T 3 ) = ( x 3 ) 2 ≥ 0. We remark that ¯ T 3 can describe one or three real points. We calculate the eigenvalues e 1 ( ¯ T 3 ) = 3 2 (1 + x 2 3 ) and e 2 ( ¯ T 3 ) = 3 2 (2 + αx 3 (1 + x 2 3 )). Hence, point ¯ T 3 always describes an unstable solution.  \nFinally, the point ¯ T 4 is not physically accepted because y ( ¯ T 4 ) > 0. Function x 4 ( α, β ) is given by the algebraic equation β = αx 4 -1 1+2 x 2 4 . The eigenvalues of the linearized system near to the stationary point are e 1 ( ¯ T 4 ) = 3 2 ( 1+2 x 2 4 ) ( -1 + √ 1 + 4 x 4 (1 + x 2 4 ) (1 + 2 x 2 4 ) ( α +4 x 4 (1 + αx 4 (1 + x 2 4 ))) ) and e 2 ( ¯ T 4 ) = 3 2 ( 1+2 x 2 4 ) ( -1 -√ 1 + 4 x 4 (1 + x 2 4 ) (1 + 2 x 2 4 ) ( α +4 x 4 (1 + αx 4 (1 + x 2 4 ))) ) . From these eigenvalues it is clear that ¯ T 4 can be an attractor. However, this will be unphysical and the free parameters α , β should be constrained such that ¯ T 4 is a source or saddle point.  \nIn Fig. 14 we present the contour plot that relates the free parameters α , β and x 4 ; we also present the region where the real parts of the eigenvalues are negative. We observe that the model is physically accepted approximately for β > -1. Phase-space portraits of the two-dimensional dynamical system (87), (88) are given in Fig. 15. Finally, the evolution of the physical parameters Ω m , Ω ϕ , w eff and of the cosmographic parameters are given in Figs. 16 and 17. The values of the free parameters have been selected such that the unique attractor is the de Sitter point ¯ T 2 .', '7. CONCLUSIONS': "Cosmological models with interaction in the dark sector are crucial for describing the universe. Due to the interaction, new degrees of freedom are introduced into the gravitational model, which drive the dynamics and help explain the observations.  \nIn this study, we examined gravitational models that describe compartmentalization and coexistence in the dark sector of the universe. Drawing inspiration from population dynamics  \nmodels in biology, we analogized dark energy and dark matter as two species evolving in the 'universe.' Furthermore, dark radiation was introduced as the third 'species.'  \nTo study the interaction models, we employed dynamical analysis and a phase-space investigation. We used the H-normalization approach and calculated the asymptotic solutions and their stability properties in terms of dimensionless variables. Each asymptotic solution corresponds to a specific stationary point of the cosmological field equations, representing a distinct era of cosmic evolution.  \nWe considered dark matter to be described by a pressureless fluid source, while for dark energy we assumed two cases. In the first case, dark energy is described by an ideal gas with a constant equation of state parameter. In the second scenario, we considered dark energy as described by a quintessence scalar field. Due to the potential function, the dynamics of the second cosmic scenario leads to a more complex cosmological history. From this analysis, we can provide constraints on the dynamical variables of the interacting model and discuss the cosmological viability of the models.  \nFrom this research, the interactions between dark matter and dark energy could significantly influence the dynamics of the universe. The study demonstrates that these interactions could potentially address cosmological tensions such as H 0 and σ 8 tensions, suggesting that the interaction terms might be key to solving some persistent problems in cosmology.  \nFurthermore, the models analysed show varied asymptotic behaviours, indicating the complexity and sensitivity of cosmic evolution to the specific characteristics of dark sector interactions. We conclude that understanding these interactions not only provides insights into the nature of dark matter and dark energy but also helps in predicting the future dynamics of the universe.  \nThe analysis and results underscore the importance of considering different interaction models and their implications on cosmological scales. In future work, we plan to investigate the extent to which such interacting models can help reduce cosmological tensions. This research lays a foundational framework for such investigations, highlighting the critical role of theoretical models in cosmology.  \nData Availability Statements: Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.", 'Acknowledgments': "AP thanks the support of VRIDT through Resoluci'on VRIDT No. 096/2022 and Resoluci'on VRIDT No. 098/2022. KD was funded by the National Research Foundation of South Africa, Grant number 131604. AA acknowledges that this work is based on research supported in part by the National Research Foundation of South Africa (Grant Numbers 151059). Part of this work was supported by Proyecto Fondecyt Regular 2024, Folio 1240514, Etapa 2024  \n- [1] D. W. Sciama, Modern Cosmology and the Dark Matter Problem, Cambridge University Press, Cambridge (1994)\n- [2] L. Amendola and S. Tsujikawa, Dark Energy: Theory and Observations, Cambridge University Press, Cambridge (2013)\n- [3] A.G. Riess et al., Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant, Astron. J. 116, 1009 (1998)\n- [4] M. Kowalski et al., The Three-Dimensional Power Spectrum of Galaxies from the Sloan Digital Sky Survey, Astrophys. J. 686, 749 (2008)\n- [5] L. Perivolaropoulos and F. Skara, Challenges for ΛCDM: An update, New Astronomy Reviews 95, 101659 (2022)\n- [6] A. Nicolis, R. Rattazzi and E. Trincherini, Phys. Rev. D 79, 064036 (2009)\n- [7] C. Deffayet, G. Esposito-Farese and A. Vikman, Phys. Rev. D 79, 084003 (2009)\n- [8] G.W. Horndeski, Int. J. Ther. Phys. 10, 363 (1974)\n- [9] L. Arturo Ure˜na-L'opez, J. Phys. Conf. Ser. 761, 012076 (2016)\n- [10] B. Ratra and P.J.E Peebles, Phys. Rev. D 37, 3406 (1988)\n- [11] A. Linde, D-Term Inflation, Phys. Rev. D 49, 748 (1994)\n- [12] E.O. Kahya and B. Pourhassa, The universe dominated by the extended Chaplygin gas, Astroph. Space Sci. 353, 677 (2014)\n- [13] O. Bertolami and V. Duvvuri, Chaplygin inspired Inflation, Phys. Lett. B 640, 121 (2006)\n- [14] C. Brans and R.H. Dicke, Phys. Rev. 124, 195 (1961)\n- [15] J. O'Hanlon, Phys. Rev. Lett. 29 137 (1972)\n- [16] K.-i. Maeda and N. Ohta, Phys. Lett. B 597, 400 (2004)\n- [17] R. Ferraro and F. Fiorini, Phys. Rev. D 75, 084031 (2007)\n- [18] R. Ferraro and F. Fiorini, Phys. Rev. D 78, 124019 (2008)\n- [19] A. Paliathanasis, J.D. Barrow and P.G.L. Leach, Cosmological solutions of f(T) gravity, Phys. Rev. D 94, 023525 (2016)\n- [20] S. Nojiri and S.D. Odintsov, Phys. Lett. B 631, 1 (2005)\n- [21] T. Clifton, Spherically symmetric solutions to fourth-order theories of gravity, Class. Quant. Grav. 23, 7445 (2006)\n- [22] L. Heisenberg, Review on f(Q) gravity, Physics Reports 1066, 1 (2024)\n- [23] L. Amendola, Phys. Rev. D 62, 043511 (2000)\n- [24] R.-G. Cai and A. Wang, Cosmology with interaction between phantom dark energy and dark matter and the coincidence problem, JCAP 03, 002 (2005)\n- [25] L. Amendola, G.C. Campos and R. Rosenfeld, Consequences of dark matter-dark energy interaction on cosmological parameters derived from type Ia supernova data, Phys. Rev. D 75, 083506 (2007)\n- [26] J. F. Jesus, A. A. Escobal, D. Benndorf and S. H. Pereira, Can dark matter-dark energy interaction alleviate the cosmic coincidence problem?, Eur. Phys. J. C 82, 273 (2022)\n- [27] N. Cruz, S. Lepe and F. Pena, Dark energy interacting with two fluids, Phys. Lett. B 663, 338 (2008)\n- [28] J.M. Salim and S.L. Saut'u, Gravitational theory in Weyl Integrable Spacetime, Class. Quantum Grav. 13, 353 (1996)\n- [29] J.M. Salim and S.L. Sautu, Gravitational collapse in Weyl Integrable Spacetimes, Class. Quantum Grav. 16, 3281 (1999)\n- [30] J. Khoury and A. Wetlman,Phys. Rev. Lett. 93, 171104 (2004)\n- [31] L. Amendola, Coupled Quintessence, Phys. Rev. D 62, 043511 (2000)\n- [32] V. Faraoni, Cosmology in Scalar-Tensor Gravity, Kluwer Academic Publishers, Dordrecht, (2004)\n- [33] G. Leon and A. Paliathanasis, Stability of a modified Jordan-Brans-Dicke theory in the dilatonic frame, Gen. Rel. Grav. 52, 71 (2020)\n- [34] Y.-F. Cai, E.N. Saridakis, M.R. Setare and J.-Q. Xia, Phys. Rept. 493, 1 (2010)\n- [35] A. Paliathanasis, Class. Quantum Grav. 37, 195014 (2020)  \n- [36] S.V. Chervon, Quantum Matter 2, 71 (2013)\n- [37] A. R. Brown, Phys. Rev. Lett. 121, 251601 (2018)\n- [38] S. Kumar and R. C. Nunes, Phys. Rev. D 96, no. 10, 103511 (2017)\n- [39] E. Di Valentino, A. Melchiorri and O. Mena, Phys. Rev. D 96, 043503 (2017)\n- [40] S. Pan, W. Yang and A. Paliathanasis, MNRAS 493, 3114 (2020)\n- [41] S. Kumar, R. C. Nunes and S. K. Yadav, Eur. Phys. J. C 79, no. 7, 576 (2019)\n- [42] A. Pourtsidou and T. Tram, Phys. Rev. D 94, no. 4, 043518 (2016)\n- [43] R. An, C. Feng and B. Wang, JCAP 1802, no. 02, 038 (2018)\n- [44] M. Lucca, Dark energy-dark matter interactions as a solution to the S 8 tension, Phys. Dark Univ. 34, 100899 (2021)\n- [45] V. Poulin, J.L. Bernal, E. Kovetz and M. Kamionkowski, The Sigma-8 Tension is a Drag, Phys. Rev. D 107, 123538 (2023)\n- [46] J. Gleyzes, D. Langlois, M. Mancarella and F. Vernizzi, JCAP 1508, 054 (2015)\n- [47] C. G. Boehmer, N. Tamanini and M. Wright, Phys. Rev. D 91, no. 12, 123002 (2015)\n- [48] C. G. Boehmer, N. Tamanini and M. Wright, Phys. Rev. D 91, no. 12, 123003 (2015)\n- [49] S. Pan, G. S. Sharov and W. Yang, Phys. Rev. D 101, no.10, 103533 (2020)\n- [50] G. D'Amico, T. Hamill and Nemanja Kaloper, Phys. Rev. D 94, 103526 (2016)\n- [51] S. Pan, J. de Haro, W. Yang and J. Amor'os, Phys. Rev. D 101, no.12, 123506 (2020)\n- [52] D. Pav'on and B. Wang, Gen. Rel. Grav. 41, 1 (2009)\n- [53] L. P. Chimento, Phys. Rev. D 81, 043525 (2010)\n- [54] W. Yang, S. Pan and J. D. Barrow, Phys. Rev. D 97, no. 4, 043529 (2018)\n- [55] M. Thorsrud, D. F. Mota and S. Hervik, JHEP 1210, 066 (2012)\n- [56] S. Pan, S. Bhattacharya and S. Chakraborty, Mon. Not. Roy. Astron. Soc. 452, 3038 (2015)\n- [57] S. Pan and G. S. Sharov, Mon. Not. Roy. Astron. Soc. 472, 4736 (2017)\n- [58] G. Caldera-Cabral, R. Maartens, L.A. Ure˜na-L'opez, Phys. Rev. D 79, 063518 (2009)\n- [59] A. Paliathanasis, S. Pan and W. Yang, Int. J. Mod. Phys. D 28, 1950161 (2019)\n- [60] M.A. van der Westhuizen and A. Abebe, Interacting dark energy: clarifying the cosmological implications and viability conditions, JCAP 01, 048 (2024)\n- [61] E. Gonz'alez, C. Maldonado, N. S. Mite and R. Salinas, WIMP dark matter in bulk viscous non-standard cosmologies (2024) [arXiv:2409.03083 [hep-ph]].\n- [62] Y.C. Ong, Universe 9, 437 (2023)  \n- [63] L. Amendola, D. Polarski and S. Tsujikawa, IJMPD 16, 1555 (2007)\n- [64] G. Leon and E.N. Saridakis, JCAP 1504, 031 (2015)\n- [65] G. Leon and E.N. Saridakis, JCAP 1504, 031 (2015)\n- [66] T. Gonzales, G. Leon and I. Quiros, Class. Quantum Grav. 23, 3165 (2006)\n- [67] A. Giacomini, S. Jamal, G. Leon, A. Paliathanasis and J. Saveedra, Phys. Rev. D 95, 124060 (2017)\n- [68] G. Acquaviva and Nihan Katirci, Dynamical analysis of logarithmic energy-momentum squared gravity, Phys. Dark Univ. 38, 101128 (2024)\n- [69] S. Mishra and S. Chakraborty, EPJC 79, 328 (2019)\n- [70] L. Chen, Dynamical analysis of loop quantum R 2 cosmology, Phys. Rev. D 99, 06425 (2018)\n- [71] A. Chattrerjee, Dynamical analysis of coupled curvature-matter scenario in viable f(R) dark energy models at de Sitter phase, Class. Quantum Grav. 41, 095007 (2024)\n- [72] G. Acquaviva and A. Beesham, Dynamical analysis of a first order theory of bulk viscosity, Class. Quantum Grav. 35, 195011 (2018)\n- [73] A. Algo, C. Uggla and J. Wainwright, Perturbations of the Lambda-CDM model in a dynamical systems perspective, JCAP 09, 045 (2019)\n- [74] S. Basilakos, G. Leon, G. Papagiannopoulos and E.N. Saridakis, Dynamical system analysis at background and perturbation levels: Quintessence in severe disadvantage comparing to ΛCDM, Phys. Rev. D 100, 043524 (2014)\n- [75] S. Pan, W. Yang, E. Di Valentino, D. F. Mota and J. Silk, IWDM: the fate of an interacting non-cold dark matter - vacuum scenario, JCAP 07 (2023), 064\n- [76] D. Benisty, S. Pan, D. Staicova, E. Di Valentino and R. C. Nunes, Late-Time constraints on Interacting Dark Energy: Analysis independent of H 0 , r d and M B , Astron. Astrophys. 688 (2024), A156\n- [77] G. Mangano, G. Miele and V. Pettorino, Coupled quintessence and the coincidence problem Mod. Phys. Lett. A 18, 831 (2003)\n- [78] N.A. Koshelev, On the growth of perturbations in interacting dark energy and dark matter fluids, Gen. Rel. Grav. 43, 1309 (2011)\n- [79] Di Valentino, E., Melchiorri, A., Mena, O., & Vagnozzi, S. (2020). Interacting dark energy in the early 2020s: A promising solution to the H 0 and cosmic shear tensions. Physics of the Dark Universe, 30, 100666.  \n- [80] Y.C. Ong, An Effective Sign Switching Dark Energy: Lotka-Volterra Model of Two Interacting Fluids, Universe 9, 437 (2023)\n- [81] I. Masina, Dark matter and dark radiation from evaporating primordial black holes, Eur. J. Phys. Plus 135, 552 (2020)\n- [82] K.V. Berghaus, P.W. Graham, D.E. Kaplan, G.D. Moore, S. Rajendran, Phys. Rev. D 104, 083520 (201)\n- [83] Di Valentino, E., Mena, O., Pan, S., Visinelli, L., Yang, W., Melchiorri, A., & Silk, J. (2021). In the realm of the Hubble tension'a review of solutions. Classical and Quantum Gravity, 38(15), 153001\n- [84] Bagherian, H., Joseph, M., Schmaltz, M., & Sivarajan, E. N. (2024). Stepping into the Forest: Confronting Interacting Radiation Models for the Hubble Tension with Lymanα Data. arXiv preprint arXiv:2405.17554.\n- [85] Gariazzo, S., & Mena, O. (2024). On the dark radiation role in the Hubble constant tension. In The Hubble Constant Tension (pp. 457-480). Singapore: Springer Nature Singapore.\n- [86] S. Weinberg, Gravitation and cosmology: Principles and applications of the general theory of relativity, (Wiley, New York, 1972).\n- [87] M. Visser, Jerk, snap, and the cosmological equation of state, Class. Quant. Grav. 21, 26032616 (2004).  \n<!-- image -->  \n<!-- image -->  \nα  \n=  \n1  \n/  \n10,  \nw  \nd  \n=-  \n1,  \nΩ  \n0  \nDE  \n=  \n0.01  \n<!-- image -->  \nFIG. 3: Compartmental: Qualitative evolution of w eff , of the deceleration paramter q and of the statefinders j and s for values of the free parameters such that the future attractor for the dynamical system (53), (54) is the point A 2 and the initial conditions are of that of Fig. 2.  \n<!-- image -->  \n<!-- image -->  \nFIG. 4: Lotka-Volterra: Phase-space portraits for the dynamical system (57), (58) for w d = 1 -β 3 , α = -3 with β = 4 (left fig.) and β = 3 (right fig.). The physical accepted stationary points are marked with red dots. The unique attactor in the two figures is point ¯  \n<!-- image -->  \nA 2 which describes a universe with a coexistence of dark matter and dark energy.  \nFIG. 5: Lotka-Volterra: Qualitative evolution of the energy density parameters Ω m , Ω DE and Ω r for the dynamical system (57), (58) for values of the free parameters such that the future attractor is point ¯ A 2 . Plots are for different initial conditions for the variables Ω m ( a ) , Ω DE ( a ) and Ω r ( a ).  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nα  \n=-  \n3,  \nw  \nd  \n=-  \n1,  \nβ  \n=  \n4,  \nΩ  \n0  \nDE  \n=  \n0.01  \n<!-- image -->  \n0  \n0  \nFIG. 6: Lotka-Volterra: Qualitative evolution of w eff , of the deceleration paramter q and of the statefinders j and s for values of the free parameters such that the future attractor of the dynamical system (57), (58) is point ¯ A 2 and initial conditions are that of Fig. 5.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFIG. 7: EMR-I: Region space where the asymptotic solution described by point B 4 is stable.  \n<!-- image -->  \n<!-- image -->  \nFIG. 8: EMR-I: Phase-space portraits for the dynamical system (65), (66) for values of the free parameters where B 4 is physically accepted. The dynamical space admits four stationary points with B 4 the unique attractor, and B 1 a saddle point.  \n<!-- image -->  \nFIG. 9: EMR-I: Qualitative evolution of the energy density parameters Ω m , Ω DE and Ω r for the dynamical system (65), (66) for values of the free parameters such that B 4 is physically accepted. Plots are for different initial conditions for the variables Ω m ( a ) , Ω DE ( a ) and Ω r ( a ).  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nα  \n=  \n1.06,  \nγ  \n=  \n0.2,  \nw  \nd  \n=-  \n1,  \nΩ  \n0  \nDE  \n=  \n10  \n-  \n5  \n<!-- image -->  \nFIG. 10: EMR-I: Qualitative evolution w eff , of the deceleration paramter q and of the statefinders j and s for values of the free parameters such that the future attractor of the dynamical system (65), (66) is the point B 4 and initial conditions are that of Fig. 9.  \n<!-- image -->  \n<!-- image -->  \nFIG. 11: Quintessence (Compartmental): Phase-space portraits for the dynamical system (82), (83) for α = -2 (left fig.) and α = 1 2 (right fig.) The stationary points are marked with red dots. We remark that the unique attractor is point T 4 . However, the structure of the phase-space is different when point T 3 exist.  \n<!-- image -->  \nFIG. 12: Quintessence (Compartmental): Qualitative evolution of the energy density parameters Ω m and Ω ϕ the dynamical system (82), (83) for α = -2, and different sets of initial conditions.  \n<!-- image -->  \n4  \nFIG. 13: Quintessence (Compartmental): Qualitative evolution w eff , of the deceleration paramter q and of the statefinders j and s as they are given by the solution of the dynamical system (82), (83) and initial conditions of that of Fig. 12. The de Sitter universe described by point T 4 is the unique future attractor.  \n<!-- image -->  \nα  \n=-  \n2,  \nλ  \n=  \n0,  \ny  \n0  \n-  \n4  \n<!-- image -->  \n4  \n<!-- image -->  \nα  \n=-  \n2,  \nλ  \n=  \n0,  \ny  \n-  \n4  \n<!-- image -->  \n0  \n=  \n10  \n=  \n10  \n<!-- image -->  \nFIG. 14: Quintessence (Lotka-Volterra): Left fig.: Contour plot which relates the values x 4 , α and β . Right fig.: Region when point ¯ T 4 is an attractor.  \n<!-- image -->  \n<!-- image -->  \nFIG. 15: Quintessence (Lotka-Volterra): Phase-space portraits for the dynamical system (87), (88) for α = -2, β = -1 4 (left fig.) and β = 5 (right fig.) The stationary points are marked with red dots. The unique attractor is the de Sitter point T 2 .  \n<!-- image -->  \nFIG. 16: Quintessence (Lotka-Volterral): Qualitative evolution of the energy density parameters Ω m and Ω ϕ the dynamical system (87), (88) for α = -2, β = -1 4 and different sets of initial conditions.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nα  \n=-  \n2,  \nλ  \n=  \n0,  \nβ  \n=-  \n0.25,  \ny  \n0  \n=  \n10  \n-  \n4  \n<!-- image -->  \n4  \nFIG. 17: Quintessence (Lotka-Volterra): Qualitative evolution w eff , of the deceleration paramter q and of the statefinders j and s as they are given by the solution of the dynamical system (82), (83) and initial conditions of that of Fig. 16. The de Sitter universe described by point T 2 is the unique future attractor.  \n<!-- image -->"}