Discriminating interacting dark energy models using Statefinder diagnostic

,


I. INTRODUCTION
Since the advent of modern cosmology and the discoveries of the late twentieth century, our view of the Universe has been firmly shifted.In this regard, current observational evidence coming from precise measurements of Supernovae Ia (SnIa) [1] [2], the large scale structure (LSS) from the Sloan Digital Sky Survey [3] and the cosmic microwave background (CMB) anisotropies [4] indicates that our Universe is flat or almost flat, expanding with an accelerated rate of expansion, and dominated by dark energy (DE) and cold dark matter (DM).The accelerating expansion of the present Universe is attributed to the DE, which is an exotic component with negative pressure, such as the cosmological constant (CC) Λ [5][6][7], which can be modeled as a perfect fluid with energy density ρ Λ = Λ/8πG and negative pressure p Λ = −ρ Λ .The Lambda-cold-dark-matter (ΛCDM) model has proven to be successful at explaining most observations of our universe [8].Despite this success, the ΛCDM model is plagued of theoretical issues.The first, and perhaps most significant problem, called the cosmological constant problem [9,10], is the catastrophic disagreement between the value of the energy density of Λ, ρ Λ , calculated from observation, and the vacuum energy density, ρ vac , predicted by quantum field theory (QFT), which is about 10 120 larger than the observed value.Since the theoretically expected value of the cosmological constant is much larger than the observed one, one may hope that Λ vanishes in some way.In this case, an alternative explanation for the current accelerated expansion of the universe is required.Such explanations can be classified in modified matter models and modified gravity models.For a review of DE models, see Refs.[11][12][13].The second significant problem with Λ is the so-called coincidence problem [14,15].The coincidence parameter r is defined as the ratio between DM and DE densities, r ≡ ρc ρ de , which at the present time has a value of r0 ∼ O (1).In most models of the universe, this situation is indeed highly coincidental, as a very specific set of initial conditions is required to yield the correct relative energy densities in the present epoch.
In addition to the already mentioned theoretical/conceptual issues, some anomalies and tensions have been found between cosmological and astrophysical data.In particular, the tension in estimation of the Hubble constant H 0 from different astronomical/cosmological observations assuming the standard ΛCDM model of the Universe has remained a big issue in the current cosmology [16].A discrepancy between 4σ and 6σ [17] in the measures of H 0 has been observed by comparing early measurements from Planck Collaboration [18] and late observations with supernovae Type Ia [19].In addition to the H 0 tension, a second tension has emerged during the recent years: the S 8 tension [20].In ΛCDM model, the S 8 parameter quantifies the amplitude of matter fluctuations in the late universe, and is defined as (Ω matter /0.3) 0.5 σ 8 , where σ 8 is the standard deviation of the density fluctuation in an 8 h −1 Mpc radius sphere.The measurement derived from low-redshift probes [21][22][23] is systematically 2 − 3 σ lower than that indirectly measured by the Planck-2018 CMB data [18].Theoretically, tensions between measurements in the early and the late universe could represent the signature of new physics beyond the ΛCDM model [24].On this subject, interacting DE models have been reexaminated in the light of latest cosmological observations, since a non-vanishing coupling between the dark sector of the universe can eventually alleviate the emerging cosmological tensions of ΛCDM model [25][26][27][28][29][30][31][32].Remarkably, the assumption that DE and DM do not evolve separately but interact with each other non-gravitationally, which allows for an energy transfer between DE and DM, was considered previously in order to solve the cosmic coincidence problem .For extensive reviews on interacting DE models, see Refs.[62,63].The coupling between DE and DM is modeled through a rate of energy exchange between both components.In the absence of a fundamental theory, this rate of energy exchange can not be obtained from first principles.Therefore, a phenomenological approach for the coupling is used [64].For a given model, it is possible to choose the parameters in order the model be consistent with data (e.g., SnIa, CMB, BAO and H(z)) that constraint the expansion history at background level [65][66][67][68][69][70].As several DE models predict very similar expansion histories, mostly of them are still in agreement with the available observational data.What is more, different mathematical model based on the same observational data could produce equivalent results [71].Thus, we need a criterion to discriminate between DE models.They should agree with current observations and predict, at least, very similar expansion histories.The above-mentioned problem was one of the main reasons to introduce alternative mechanisms to identify models which could be quite similar than other.On this subject, the Statefinder diagnostic was introduced in 2003 by [72,73].Such approach is a tool that makes possible to differentiate between the very distinct and competing cosmological scenarios involving DE at background level.The Statefinder diagnostics is performed taking advantage of higher order of derivatives on the scale factor a(t), in such a way that it introduces the so-called Statefinder parameters, r, s, which are corrections to the Hubble rate H(z) and deceleration parameter q(z) of higher order in a(t).Thus, by computing the correction behind H(z) and q(z) we could distinguish between models apparently equivalent.This diagnostic proposal introduces new geometrical dimensionless parameters that characterize the properties of DE regardless of the model, because they depend on the observable Hubble parameter and its derivatives [74].Statefinder diagnostic have been applied to several models, e.g., some authors have applied the Statefinder diagnostic to discriminate holographic DE (HED) models from the ΛCDM model [75], to analyze Barrow holographic DE [76], interacting DE models [77][78][79][80], among others [81][82][83][84][85].
Considering that an interaction between DE and DM offers an interesting framework to study phenomenology beyond to ΛCDM model, the main goal of the present work is to analyse several interacting DE, recently proposed in the literature, via Statefinder diagnostic.In doing so, for each interacting DE model we compute the Statefinder parameters as functions of the redshift, studying their high and low-redshift limits.We also plot the evolution trajectory on the s − r and q − r planes in order to determine its deviation from ΛCDM.In addition, according to their behaviour on the s−r and q −r planes, we discriminate between the several interacting DE studied, breaking the degeneracies of the models at background level.Basically, we perform the Statefinder diagnostic for 17 several forms of the energy transfer rate Q between DE and DM, belonging to the following categories: i) linear models in energy densities of DE and DM, ii) non-linear models, iii) models with a change of direction of energy transfer between DE and DM, iv) models involving derivatives of the energy densities, v) parametrized interactions through a function of the coincidence parameter r, and finally we also consider two kinds of models with a self-interaction between DM, without DE.These models have been already studied in the literature and constrained with observational data available at that time, see e.g.[25, 26, 30, 54-56, 60, 62, 63, 86-92].
Our work is organized as follows: after this introduction, we present a brief review of the ΛCDM model as well as the basic equations which describe an interaction between DE and DM in Section II.In the third Section, we present the Statefinder parameters which will be applied for the several models, while in section IV the background dynamics of each particular model is studied.In particular, analytical solutions for the corresponding Hubble rate as functions of the redshift H(z) are presented.In the fifth section we discuss the Statefinder analysis for the several models, and in Section VI we make the comparison between models in light of the Statefinder analysis.Finally we summarize our findings and present our conclusions in Section VII.We adopt the mostly positive metric signature, (−, +, +, +), and we work in natural units where c = ℏ = 1.

A. ΛCDM model
We consider a FLRW Universe where a(t) is the scale factor at cosmic time t and K is the Gaussian curvature of the space-time.Setting κ 2 = 8πG, with G being the Newton's constant, the scale factor satisfies the Friedmann equations [93] where H = ȧ/a is the Hubble rate and ρ A and p A denote the energy density and pressure of each individual fluid component, respectively, and they are related by an en equation-of-state (EoS) parameter w A = p A /ρ A .The curvature term can be brought to the right-hand side of Eq.( 2) by defining ρ K = −3K/(8πGa 2 ).On the other hand, the energy density associated to Λ can be written as ρ Λ = Λ/8πG, while its EoS becomes w Λ = −1.In addition to the curvature and cosmological constant terms, the ΛCDM model also includes other energy density components: baryons (b), cold DM (c), and radiation (r), characterized by the EoS parameters w b = w c = 0 and w r = 1/3, respectively.In ΛCDM it is assumed that the different matter components do not have any other interaction, therefore the energy conservation equation for each component holds Then, by integrating Eq.( 7) for baryons, cold DM, and radiation, the Friedmann equation for ΛCDM (2) can be written as where we have used the fact that each energy density is usually expressed in terms the dimensionless density parameter, defined as Ω A = ρ A /ρ cr with ρ cr = 3H 2 /(8πG) being the critical density, and the sub-index "0" denotes the current value of any given quantity.The Hubble constant is usually expressed as H 0 = 100h km s −1 Mpc −1 , where the parameter h denotes the observational uncertainty.In Eq.( 5), Ω m0 denotes the current contribution of non-relativistic matter (baryons+cold DM).By convention, the dimensionless density parameter associated with the curvature term is Ω K = −K/(aH) 2 [94,95].The cosmological parameters of the ΛCDM model, derived from the CMB temperature fluctuations measured by Planck Collaboration with the addition of external data sets are given by [96] where the spatial flatness (K = 0) has been assumed, while the contribution of radiation becomes Ω r0 ∼ 10 −5 .By using the relation between the scalar factor and the redshift a = (1 + z) −1 , the Friedmann equation ( 5) for ΛCDM with zero curvature may be written as where E(z) ≡ H(z)/H 0 is defined as the dimensionless Hubble rate.At z = 0 (today), the following constraint is satisfied If we set Λ = 0 and consider a more general DE (de) component as the responsible of the accelerated expansion, having energy density ρ de and EoS w de , a combination of CMB, SN, and BAO measurements, assuming a flat Universe, found w de = −1.03± 0.03 [18], consistent with the cosmological constant case w Λ = −1.

B. Interaction between DE and DM
Now, we are going to consider the case where an interaction between a dynamical DE component (with and EoS parameter w de ) and DM is allowed.Despite the additional postulated interaction, the total energy density of the dark sector is conserved, and the energy momentum tensor of the dark sector,T µ ν = T µ (de)ν + T µ (c)ν , satisfies therefore, if an interchange of energy between DE and DM is introduced, according to [62,66].
where the F ν is the four-vector of interaction between dark components.We can project the Eqs.( 9) parallel to the four-velocity u µ , and we obtain and respect part orthogonal to the velocity, we must use the projector In addition, the general expression for the energy-momentum tensor tensor describing perfect fluid is [97], [95] Now, we use Eqs.( 12), ( 10) and ( 11) to obtain the following Euler equations In the context of a flat FLRW universe, u µ = (1, 0, 0, 0) in the comoving coordinates, thus We use the notation u µ F µ = Q.Now, we write the Eqs.( 13) and ( 14) in the form Where Q is the rate of energy exchange between DE and DM.Here, Q > 0 reflects that the energy flows from DE to DM, while Q < 0 the opposite.In a spatially flat FLRW universe filled with DE, DM, baryons, and radiation, we can split the equation ( 7) for baryons and radiation, yielding the following set of equations ρr + 4Hρ r = 0, In phenomenological construction of interacting DE models, it is assumed that Q can be modeled as a certain function of energy densities and the Hubble rate [62,63].Several expressions for the rate of energy exchange Q have been studied in the literature: models where Q is a linear function (widely studied) of densities of DM or DE [26,30,54,86,[89][90][91], nonlinear interactions [54,60,62,91], proportional to deceleration parameter [87], and others combinations, where Q is constructed from derivatives of the densities of DM or DE, or varying couplings with dependence on cosmic time or redshift [86], among others.The conservation equations for baryons and radiation, Eqs.( 21) and (22), respectively, are solved by separation of variables, yielding However, the solutions of the equations ( 19) and ( 20) depend on the mathematical structure of Q, therefore, it is not always possible to obtain an analytical solution for E(z) = H(z)/H 0 .Then, we focus on interacting DE models that allow analytical solutions for E(z) in order to apply the Statefinder analysis.

III. STATEFINDER DIAGNOSTIC
In Ref. [73], the authors introduced the Statefinder parameters by expanding in Taylor series the scale factor a(t).We are interested in small values of |t − t 0 |.The scale factor a(t) is expanded as follows or in the form [98] ( where ) is nth derivative of the scale factor with respect to cosmic time.A n is used to define the Statefinder parameters q, r and s.Recalling that q is defined above as the deceleration parameter, it can be written in terms of A 2 as We rewrite Eq. (3) in terms of the deceleration parameter as follows where w i = p i /ρ i , Ω i = 8πGρ i /3H 2 0 and we have assumed that Ω b + Ω c + Ω de = 1.The parameter r is the next (after the Hubble rate H and the deceleration parameter q) member of the set of kinematic parameters that describe the expansion of the Universe [99], which is defined as or in terms of DE density where the parameter r was first introduced in [100], and in [101] as well as the jerk 'j'.On the other hand, the parameter s is defined as a combination of r and deceleration parameter, and it should not to be confused with the snap (the fourth time derivative) [102] s substituting q and r, It is observed that s does not explicitly depend on the DE density.The reason of why s is chosen in such a way is due to the fact that the features chosen for the description of DE may be geometrical, if they come directly from the metric of space-time, but also physical, if they depend on the characteristics of the fields that represent DE.Physical qualities vary on models, while geometrical characteristics are universal [99].
For our analysis, it will be useful to express these parameters in terms of the dimensionless Hubble rate and the redshift z.The cosmic times derivatives are written as redshift derivatives according to where E 2 (z) = H 2 (z)/H 2 0 as been used.In this way, the set Statefinder parameters becomes Therefore, if the role of DE is played by a cosmological constant, i.e w de = −1, then the value of r remains at r = 1 throughout the matter-dominated epoch and at all future times (for z ≲ 10 4 ).In the s − r plane, the fixed point {s, r} = {0, 1} in a Universe containing a cosmological constant and non-relativistic matter corresponds to ΛCDM.
According to [73,103], a first limit case (Ω m = Ω c + Ω b = 1, Ω de = 0), without radiation, corresponds to the standard cold dark matter (SCDM), in which the Universe presents a decelerated power-law expansion, according to a(t) ∝ t 2/3 .In this case, the fixed point is found to be {s, r} = {1, 1}.The other limit case (Ω m = 0, Ω de = 1) corresponds to a dS Universe, which expands at constant rate, a(t) ∝ exp Λ 3 t.The following plots correspond to the case for ΛCDM model.In Figs.(1a), (1b), (1c) we show the plots of Statefinder parameters against redshift q(z), r(z) and s(z), and in Figs.(1d), (1e) are shown the parametric plots in s − r and q − r planes.In addition, the fixed point {q, r} = {−1, 1} in the q − r plane corresponds to the asymptotic dS solution [74].By comparing the evolution trajectories on the r − q and r − s planes between the corresponding to several DE and ΛCDM models, it is possible to find some similarities and analyze the deviation and compatibility with ΛCDM model.

IV. INTERACTING DE MODELS: ANALYTIC SOLUTIONS FOR H(z)
In this Section, we present several mathematical expressions for the coupling Q between DE and DM as well as the corresponding solutions for the dimensionless Hubble rate E = H(z)/H 0 , which will be used in order to analyze each model by means the Statefinder diagnostic.First, we present models 1, 2, where Q is linear with respect to the energy densities of DM and DE, respectively.Model 3 is a linear combination of DE and DM energy densities.Models 4, 5, and 6 are non-linear.Models 7, 8, and 9 involve a coupling proportional to decelerating parameter q.Model 10 considers a linear combination of the derivatives (respect to ln(a 3 )) of the energy densities DM and DE.Furthermore, models 11, 12, 13, 14, and 15 are interactions that can be parameterized as a function of the ratio r between density DM and DE.Finally, as an "extra bonus", models 16 and 17 are two kinds of special models without DE, where exists a self-interaction between DM.It should be noted that all the investigated models are dimensionally consistent with respect to the expression for Q.In this paper, all the expressions for Q we focused on may written in general terms as Q = 3HγR(ρ c , ρ de ), where γ is a dimensionless constant and R is a real function having units of energy density.Moreover, each Q is linear in the Hubble rate H, since this facilitates the consistency of units and the solving of differential equations by replacing Q in Eqs.(19) and (20).On the other hand, we will assume w = const.where it corresponds.

A. Linear interactions
Going further the simplest interacting DE models where Q is constant, linear models have brought a lot of interest in their study.Accordingly, we are going to study three linear models, Q proportional to the DM density, the second FIG.1: Statefinder parameters plots q(z), r(z), s(z), q − r and s − r for ΛCDM.
proportional to the DE, and a general case being a linear combination of DE and DM energy densities.The models with the rate of energy transfer Q ∝ ρ c and Q ∝ ρ de have been widely studied in the literature.We use the expressions presented in Refs.[90], [59], where Qis proportional to the constant 3γ, since this simplifies the expression for E 2 (z).

MODEL Q1 = 3γHρc
For the first case to be presented, the rate Q is proportional to the DM energy density where γ is a dimensionless constant and gives us the strength of the coupling.If we set K = 0, the Friedmann equation ( 2) becomes Assuming that w de is a constant, the authors in Refs.[57,90], provide the analytical solutions for the densities ρ c and ρ de by inserting Eq.( 36) into Eqs.(19) and (20), thus If we substitute Eqs. ( 23), (24), the expression for ρ c , ρ de and a = (1 + z) −1 into (37), and dividing by the critical density, the expression for the dimensionless Hubble rate E(z) = H(z)/H 0 is expressed as Besides, the function E(z) for w de = −γ is a special case to consider.We calculate the limit for that purpose.When w → −γ the expression for the dimensionless Hubble rate becomes 2. MODEL 2. Q2 = 3γHρ de The second interaction within this category depends linearly on the DE energy density In the same way as the first interaction, according to Refs.[57], [90], [59], if Eq.( 36) is replaced into Eq.(19), the solution for ρ de is found to be Accordingly, the expression for the dimensionless Hubble rate yields where Ω m = Ω c + Ω b , Ω r = 2.469 × 10 −5 h −2 (1 + 0.2271N ef f ) and N ef f = 3.046.If γ < 0, this implies that the energy transfer is from DM to DE, and for MODELS 1 and 2 if we set γ = 0, there is no interaction and both models reduced to wCDM, while for w de = −1, wCDM reduces to Eq. ( 6) corresponding to ΛCDM.
Analogously to the first case, if w de → −γ, the expression for the dimensionless Hubble (44) takes the form Our third model is a linear combination of DE and DM energy densities which has been studied in Refs.[54,57,86,88,89,104].The coupling parameters λ c and λ de denote the strengths of the interaction between the components of the dark sector and their signs imply the direction of energy transfer.It is relevant to mention that in [88], the analytic solution found for this interaction assumed that the coupling parameters are very small (i.e.|λ c | ≪ 1 and |λ de | ≪ 1 ), and the authors neglect the product λ c λ de in the solution, however, we will consider this product to be not null when using the expressions in the subsequent analysis.In Ref. [89], the conservation equations (19,(20) were unified and the authors derived the following master equation when the coupling ( 46) is considered which is a second order differential equation, where ρ T = ρ c +ρ de and primes denote derivatives respect to the variable x = ln(a 3 ).In order to solve the above equation, for w de = const., the solution of Eq. ( 47) according to Refs.[54], [89] is found to be where ρ 1 , ρ 2 are integration constants, and m ± are considering (λ c + w de + λ de ) 2 − 4λ c λ de > 0 to obtain m + and m − real and distinct.From Eq. ( 2) and neglecting the contribution of radiation, the Hubble rate takes the following form The explicit analytic solutions for DM and DE energy densities are respectively.Furthermore, in terms of the density parameters for the equivalent two fluids Ω 1 and Ω 2 , the density parameters for DM (Ω c0 ) and DE (Ω de,0 ) are respectively written as where Ω i = 8πGρ i /3H 2 0 .Using the aforementioned equations ( 53) and ( 54), the parameters Ω 1 and Ω 2 may be expressed.Taking into account Ω c0 + Ω de,0 = Ω 1 + Ω 2 , and from Friedmann equation ( 2) at the present time t = t 0 .Then, it is obtained For a flat Universe, i.e.Ω K = 0, we have that If we write Eq. ( 50) in terms of the normalized densities and replace Ω c , Ω de and a = (1 + z) −1 , the dimensionless Hubble rate is expressed as follows and expanding it, we get It is straightforward to check that if we set λ c = λ de = 0 and w de = −1, there is no interaction and the model reduces to ΛCDM.

B. Non-linear interactions
The next three non-linear interactions we have already been studied in [56,57,62].First, we recall that the coincidence parameter r is defined as the ratio between DM and DE energy densities If ρ = ρ c + ρ de , the individual energy densities are written as This allows us to obtain a differential equation for r, which makes it easier to find its solution, as we will show later in the parameterized interactions.
1. MODEL 4. Q4 = 3Hγ ρcρ de ρ For the first non-linear interaction where ρ = ρ c + ρ de .According to [56], the following analytic solutions were found r = r0 a 3(w de +γ) , Now, we construct the expression for the dimensionless Hubble rate from the Friedmann equation ΛCDM is recovered for γ = 0 and w de = −1.For z ≫ 1 (matter-dominated epoch), we obtain the expected evolution of density ρ ∝ a −3 .In the asymptotic case, if w de → −γ, the above expression (67) becomes and when γ = 1, we have The second non-linear term to be studied has the following dependence on the energy densities where ρ = ρ c + ρ de .In this case, the analytical solutions for the ratio r and total energy density ρ are r = r0 then, the expression for the dimensionless Hubble rate is (73) For a ≪ 1 (the high-redshift limit), the ratio r becomes a constant, r → |w de | /γ.In the opposite limit (a ≫ 1), ρ ∝ a −3 , as in the ΛCDM case.Similarly for MODEL 4 when γ = 0 and w de = −1, ΛCDM is recovered.A special case corresponds for w de → −1 and γ → 1.The expression for the dimensionless Hubble rate (73) exhibits a convergent behavior to The third interaction term belonging to this category is where ρ = ρ c + ρ de .For w de < 0, i.e. w de = − |w de | the solutions are while the dimensionless Hubble rate yields In this case, the expression ( 77) is always defined for any values of γ and w de .Additionally, the ratio r scales as a −3|w de | for a ≪ 1.For w de = −1, this coincides with the scaling of its ΛCDM counterpart.In the far-future limit, however, we obtain r ⇒ γ/ |w de |, (a ≫ 1), i.e., the energy-density ratio remains finite, while it tends to zero in the ΛCDM model.The energy density scales as a −3 for a ≪ 1, i.e., we recover an early matter-dominated period.In the limit a ≫ 1 one has From this last Eq., we can see that this generally does not correspond to a dS phase.
C. Models with a change of direction of energy transfer In Ref. [105], the authors investigated the interaction regardless the explicit form of Q, by using the most current observational data available at that time.The authors divided the whole redshift range z into a few and concluded that the interaction term δ(z) = Q/(3H) was a constant in each redshift interval.They realized that δ(z) is likely to cross the line where the interaction does not occur (δ = 0), therefore, the sign of the interaction Q must change in the approximate redshift range of 0.45 ≲ z ≲ 0.9.This result poses a problem since most of the interactions studied in the literature are of linear type as in MODELS 1 and 2, which are always positive or negative and, therefore, cannot give the possibility of changing their signs.Taking into account the work [105], the author proposed a new type of interaction in [55] and [87], where the deceleration parameter q was included in the coupling Q, thus allowing that the interaction Q changes sign when the Universe goes from deceleration (q > 0) to acceleration (q < 0).In this Section, we will investigate MODELS 7, 8, and 9, previously studied in [87], and later by [62,92].In these cases, Q is proportional to the deceleration parameter, which produces a change of direction of energy transfer.For simplicity, in these works, the authors set w de = −1.
1. MODEL 7. Q7 = q(α ρc + 3βHρc) The first interaction within this category is a linear combination of energy density of DM and its time derivative For this kind of coupling and using Eq.( 20) and ( 3), a transcendental differential equation for H(a) is obtained which cannot be solved analytically for α ̸ = 0. Nevertheless, in the particular case of α = 0, the analytic solution for the dimensionless Hubble rate is found to be When β = 0, it yields corresponding to ΛCDM, where the contribution coming from baryons and radiation components is neglected.We note that if β → −1 or β → − 2 3 , the expression (81) becomes undefined.
2. MODEL 8. Q8 = q(α ρtot + 3βHρtot) In this second case, the coupling involves the sum of DM and DE densities, i.e ρ = ρ c + ρ de and its time derivative Similarly to the previous case, we arrive at a transcendental differential equation for H(a) for α = 0, the interaction reduces to Q = 3βqHρ tot , and the corresponding solution for the dimensionless Hubble rate yields where r 1 ≡ 4 + β(4 + 9β) and If β = 0, E(z) reduces to which corresponds to ΛCDM.
3. MODEL 9. Q9 = q(α ρde + 3βHρ de ) Now, we are going to present a model where Q is a linear combination of DE energy density and its time derivative In the same way as two previous cases, we arrive at a transcendental differential equation for H(a) We solve the equation for α = 0 and obtain the following dimensionless Hubble rate where r 2 ≡ |2 − β| and If we set β = 0 ΛCDM model is recovered.If β → 2 3 , the expression (90) becomes D. Interaction involving derivatives of the energy densities In Ref. [59], the authors studied an interaction where Q is constructed to be a lineal combination of DE and DM energy densities but involving derivatives (derivatives with respect to ln(a 3 ), where a is the scale factor) as where ′ ≡ d dln(a/a0) 3 = d 3Hdt .This model was previously studied in [54] and later in [106].Assuming that w de is constant, in Ref. [59] the authors provide an analytical solution for the total energy density ρ = ρ de + ρ c and hence the analytical expression for E 2 (z).Firstly, the solution for ρ is given by where and We assume that (1 − α) 2 w 2 de − 4αw de > 0 for β ± to be real and distinct.If de − 4αw de = 0.In that case, F ± is undefined.Thus, the dimensionless Hubble rate can be written as It is easy to verify that if α = 0 there is no interaction.We obtain from Eq. ( 96), β + = −1, and β − = −(1 + w de ), therefore F + = −Ω x and F − = Ω c , therefore, the model reduces to wCDM.

E. Parameterized interactions
In the following Subsection, we present some interactions which are derived from a real function f (r) of the coincidence parameter r is introduced, which is the ratio between the energy densities of DM and DE, and it is given by Eq.( 60).The main motivation for studying the interactions parametrically is to solve the cosmic coincidence problem.By solving the differential equations one has an analytical solution for the ratio, allowing us to see the evolution of the ratio as a function of the scale factor or redshift.Some of the interactions examined in the preceding models can be obtained by appropriately selecting the parametric function f (r).In Ref. [60], five models were investigated by the authors.In this reference is set w de = −1 and the sign of Q is changed, hence, the equations ( 19) and ( 20) have the following form where Q is written as Q = 3HγR(ρ c , ρ de ), whit γ being a dimensionless constant and R is a real function which has units of energy density.By taking the time derivative of the ratio r between the energy density of DM and DE as Eq. ( 60), we obtain the expression Using the equations ( 99), ( 100) and ( 101) together, we get where We are interested in the analytical solution for (102), and we can obtain the expression for ρ c = ρ c (a), ρ de = ρ de (a) when r = r(a).
The following models correspond to different choices of f (r).
The first parametric form is chosen to be which is equivalent to the interaction of MODEL 4, given by Eq.( 62) Besides, the coincidence parameter has the following solution while the energy density of DM yields By using the relation a = (1 + z) −1 , we have the expression for the dimensionless rate Hubble as follows It is important to mention that when we analyze the models numerically in next Section, the sign of γ is the opposite of that in the models with equivalent interactions since the sign for Q has been chosen with the opposite sign.
In the same way as we did in MODEL 4, we investigate the case where In such a case, the dimensionless Hubble rate (108) takes the form which is equivalent to the equation (69).
2. MODEL 12. f (r) = 1 r By choosing the following explicit form for f (r) results in a coupling Q which is equivalent to the interaction of MODEL 6 given by Eq.( 74) and the solution for r(a) is Besides, the analytical solution for the DE density is given as while the dimensionless Hubble rate is expressed as 3. MODEL 13. f (r) = r An equivalent interaction as MODEL 5, given by Eq.( 70) can be obtained if f (r) depends linearly on the coincidence parameter, i.e.
For this parametrization, the coincidence parameter has the solution and the DM energy density results in In this way, the dimensionless Hubble rate is written as For a parameterization of the type the corresponding rate Q becomes equivalent to MODEL 2, which is given by Eq.( 42) and the coincidence parameter may be written as follows Therefore, the dimensionless Hubble rate becomes The interaction of the MODEL 1 can be reproduced if f (r) has the following explicit dependence on r For this last parametrization, the solution for the ratio r is found to be r(a and accordingly In this model, the dimensionless Hubble rate yields We observe that E 2 (z) reduces to ΛCDM when γ = 0.

F. Self-interaction between DM
In Ref. [107], the authors instead of considering a DE component, they argue that a self-interacting DM produces the accelerated expansion.They investigate background cosmic evolution by assuming a direct interaction between two DM particle types.One species is denoted by a m subscript, while the other one is denoted by a x subscript.It is important to point out that baryons and radiation are not considered in this work as self-interacting models for small redshift will be studied.By defining a Langrangian density, conservation equations for each type of DM particles are derived, which are similar to the equations ( 19) and ( 20), thus In that work, two specific models are studied: a symmetric model and an asymmetric model with respect to the Q m and Q x interaction rates.In the first model, a simple linear interaction between the both DM components is proposed, and it is similar to those investigated between DE and DM interacting models, with the benefit of knowing the mathematical behavior through the analytical solution of the density equations.As a result, the interaction rates are defined as follows thus, solving Eqs. ( 129) and (130), and assuming p m = p x = 0, it was obtained the following dimensionless Hubble rate where Ω i = ρ 0 i /3H 2 0 .As it can be seen, the equation for E 2 (z) is symmetrical for both types of DM components.It is crucial to note that when α, β ≪ 1 (but non null) implies that E 2 (z) reduces to This indicates that having a small (but non-zero) interaction between two dust-like components, the model we obtain is very similar to a model with no interaction containing dust and a cosmological constant.Linear interaction functions may be defined in an asymmetric way, such that which implies that the solution for the corresponding Hubble rate yields In the limit z → −1 for the expression E 2 (z), we get clearly the term between the square brackets vanishes when z → −1.It is concluded that, regardless of the value of β, in the limit when z → −1, ΛCDM is recovered.

V. STATEFINDER ANALYSIS OF INTERACTING DE MODELS
Previously, we have shown how to obtain the dimensionless Hubble rate E(z) for the several interacting models.Now, in the present Section we will compute the Statefinder parameters q(z), r(z), and s(z) for all the interacting DE models presented in previous section.We present the equations for q(z), r(z), and s(z) for each model and will not be shown due to that, generally, they are quite large.We will also obtain the present and asymptotic values of them as it is shown in the Table (I).
We will analyze the models for three important epochs, namely, early epoch (z ≫ 1), present epoch (z = 0), and future epoch (z → −1).As we focus on the evolution for late times, particularly in the transition from the matterdominated era to the present epoch, we neglect the radiation component in the numerical calculation of E(z) and the Statefinder {q, r, s}.In the subsequent figures, we plot the behavior of the evolutionary trajectories in the q − r and s − r planes, in order to compare the performance of the q − r, s − r parametric graphs for each model with ΛCDM.The plots corresponding to (1d) and (1e) will overlay with the curves evaluated with the fitted values of the normalized densities and constants for the respective model.With the values of the Statefinder parameters that we have calculated, we will be able to establish the deviations of this model from ΛCDM [103].On this subject, we will say that an interacting model is compatible with ΛCDM when for the matter-dominated era, the present time, and in the future, the Statefinder parameters are exactly or very close to those predicted by ΛCDM.When this is only satisfied for one o two epochs, we will say that it is partially compatible and if the Statefinder parameters deviate significantly from those predicted by ΛCDM at all epochs or if there exists a large deviation in the matter-dominated era, we will say that the model is incompatible.In order to draw the plots of H(z) with error bars, we use the values from Ref. [108] with 0 ≤ z ≤ 2.36.
A. Linear models We extract the values of the normalized densities from Ref. [59] in order to evaluate this model numerically.We have employed the values h = 0.669, Ω c = 0.301, Ω b = 0.049, γ = 0.071, w de = −1.03,corresponding to SNIa+H(z)+BAO+f gas +CMB, Table II, column C. As we can see in the graph of Fig. (2a), the number of measurements for high redshift is decreasing, also we can see that for redshift higher than 0.5 the uncertainty of the measurements increases.In general the H(z) curves obtained for each model fit accurately to the points with error bars.Of these, we will only show a few for each model category since they are quite similar.In Fig. (2b) it is shown the parametric curve of MODEL 1 in the (q, r) plane.For early times the trajectory starts in the lower right corner of the figure, in the region bounded by 0 < q < 1 and 0 < r.When z ≫ 1, the parameters {q, r, s} → {0.42, 0.79, 0.93}, hence, the values for q, r, and s deviate significantly in the standard matter-dominated era from ΛCDM, since for ΛCDM, the values for q, r and s are given by q = 0.5, r = 1 and s = 0, respectively.However, the value obtained for s is close to 1, which corresponds to SCDM [103], [73].We have evaluated the model in the past at z = 0.8 as a benchmark.For the subsequent parametric plots of all models, we will follow the same procedure; we will use different redshift values associated to the past.We plot that representing the interaction between DE and DM (dashed line), as well as, the Hubble rate for ΛCDM model (solid line), along with the data from Ref. [108].We have used H0 = 66.9 km s −1 Mpc −1 , from [59].For this model, Figures (b) and (c) show the q − r and s − r, planes, respectively.The arrows on the curves show the direction of evolution.
point, in order to better observe the temporal evolution of the model.At the present time, z = 0, we obtain {q, r, s} = {−0.50,0.99, 0.003}.It is evident that the curve has a linear behavior and evolves in the negative direction of q and positive of r.We note that close to the present time, the curve crosses the red straight line r = 1 corresponding to ΛCDM.The value of q deviates significantly from the value of q 0 ≃ −0.55, associated to ΛCDM.The values of r and s are very close from the values of r 0 = 1 and s 0 = 0, obtained for ΛCDM.At late times, the curve evolves and moves away from the dS fixed point in the future (q = −1; r = 1), as we have calculated in the limit z → −1, we have {q, r, s} → {−1.0, 1.1, −0.03}.Thus, we conclude that the model does not converge to ΛCDM in the future.In Fig. (2c) we see the (s, r) plane and we note that the plot starts in the region bounded by 0 < s < 1 and r > 0.
The trajectory is non-linear and proceeds in the negative direction of s and positive of r.We corroborate that around the present time, it crosses the fixed point (0, 1) corresponding to ΛCDM.According to our analysis, we will say that this model is partially compatible with ΛCDM.

MODEL Q2 = 3γHρ de
In order to evaluate this model numerically, we used the values of the normalized densities from Ref. [59].We will use SNIa+H(z)+BAO+f gas +CMB, Table I, column C. h = 0.671, Ω c = 0.300, Ω b = 0.048, γ = 0.07, w de = −1.05.When z ≫ 1, the parameters {q, r, s} → {0.5, 1, 0.02}, and we see that the value of q and r correspond to the standard matter-dominated era, and s is very close to the value s = 0 for ΛCDM.We observe in the Figs.(3a) and (3b) the asymptotic behaviour of r(z) and q(z) for high redshift.For MODEL 2 in the (q, r) plane, we show in Fig. (3c) that the trajectory is linear and starts in the upper right corner, in the region bounded by 0 < q < 1 and 0 < r < 1.Moreover, the curve evolves in the negative direction of q and negative of r.At the present time, z = 0, we obtain {q, r, s} = {−0.53,0.94, 0.02}, which are very close values for the Statefinder parameters for ΛCDM.When z → −1 the parameters {q, r, s} → {−0.97, 0.91, 0.02}.Thus, it is observed a significant deviation from the dS fixed point in the future.In the (s, r) plane, Fig. (3d) shows that the curve starts to the right side of the point corresponding to ΛCDM, in the region bounded by 0.01 < s < 0.03 and 0 < r < 1.The trajectory is linear in behaviour and proceeds in the negative direction of r, while holding fixed s = 0.02.It is clearly observed that the behaviour of the curve in the s − r plane is very similar to those of Quintessence models [72], [109].It is observed that s remains constant, while r decreases asymptotically up to r = 0.91.We corroborate the information provided by the (q, r) plane; the trajectory of model deviates from ΛCDM.It is concluded that this model is partially compatible with ΛCDM.
3. MODEL 3. Q3 = 3λcHρc + 3λ de Hρ de For this model we extract the values of the normalized densities and constants from Ref. [89].Table III.Ω k = 0, Ω m0 = 0.285, λ c = 0.039, λ de = −0.024,w de = −0.987.When we plot the curves on r − s and r − q planes, the baryons are considerate by separate (we choose for the model Ω b0 ≈ 0.04) and the model has the same behaviour.In [89] the baryons are not considered separately in the tables.Also, radiation is not considered in their work but we consider it in the equations, however, we set it to zero when numerically evaluating the model.For z ≫ 1, the parameters q, r and s behaves as {q, r, s} → {0.44, 0.84, 0.96}.In this case q, and r deviate significantly in the standard matter-dominated era for ΛCDM, but the value s is close to 1, which corresponds to SCDM.Fig. (4a) shows the trajectory of MODEL 3 on the (q, r) plane.This curve starts in the lower right corner of the figure, in the region bounded by 0 < q < 1 and r > 0. The trajectory has linear behaviour and proceeds in the negative direction of q and positive of r.At the present time z = 0, we obtain {q, r, s} = {−0.56,0.99, 0.005}, which are very close values to, those expected for ΛCDM.
After the present time, the trajectory crosses the red straight line r = 1 corresponding to ΛCDM, then it evolves very close to the dS fixed point in the future, which is consistent with the computation in the limit z → −1.In this case the parameters {q, r, s} → {−1, 1.1, −0.012}.Fig. (4b) shows the (s, r) plane, where the curve starts in the region bounded by 0 < s < 1 and r > 0. The trajectory is non-linear and evolves in the negative direction of s and positive of r.We observe that after the current time, it crosses the fixed point (0, 1) corresponding to ΛCDM.This model is partially compatible with ΛCDM, and it behaves similar to MODEL 1, which can be seen in the two parametric plots.This suggests that in the linear combination, the DM contribution is larger than the coming from DE.
B. Non-linear interaction 1. MODEL 4. Q4 = 3Hγ ρcρ de ρ We observe in the H(z) plot that for this model the curve moves away from the H(z) values with the error bars for redshift greater than 2. The values we will use to evaluate this model are obtained from the Ref. [60], which we will cite below to investigate models where interactions are parameterized.For this model, we have used the values Ω m0 = 0.321, γ = −0.06,w de = −1, using SNe Ia + H 0 + CC + BAO, from Table 3 of [60].However, in order to evaluate the non-linear MODELS 4, 5, and 6 numerically, we will use the value of γ with the positive sign, since in [60], the sign of Q is changed when defining the conservation equations.Furthermore, the contribution coming from baryons is not considered separately in Table 3 of values in same work.When we evaluate our model numerically, we neglect baryons and assign all matter to DM, but to evaluate the model for ΛCDM we will separate baryons.When we consider baryons separately and evaluate our model numerically, the behaviour of the curves does not change significantly.As it can be seen from Fig. (5b), that for MODEL 4 in the (q, r) plane, the curve starts in the region bounded by 0 < q < 1 and 0 < r < 1, in the upper right corner.For z ≫ 1, the parameters q, r, and s behaves as {q, r, s} → {0.5, 1, 0.06}, where q and r are exactly as ΛCDM predicts, and s is very close to zero, in this case we recover the standard matter-dominated era.The trajectory is non-linear in behaviour (like a concave-up parabola) and evolves in the negative direction of q and at r drops to the midpoint of the graph and then rises.Moreover, at the present time, z = 0, we obtain {q, r, s} = {−0.52,0.94, 0.019}.The value of q = −0.52 is close to the value for ΛCDM, q 0 ≃ −0.55, and the values of r and s are also close to 1 and 0, respectively.After, in the limit z → −1, the model gives {q, r, s} → {−1, 1, 0}, and the curve converges to the dS fixed point in the future.In the (s, r) plane, the plot starts in the region bounded by approximately 0.05 < s < 0.07 and 0.9 < r ≤ 1.The trajectory is very similar to those of the (q, r) plane, with a non-linear performance.We check the information given by the (q, r) plane, and it is observed that the model curve converges to ΛCDM in the limit z → −1.
We conclude that this model is compatible with ΛCDM, since in all three epochs the model is very similar to ΛCDM.It is possible that with more updated data the values for the Statefinder parameters {q, r, s} at present may be closer to {−0.55, 1, 0}.In order to evaluate this model numerically, we will use the following values: Ω m0 = 0.320, γ = −0.038(recall that, we will use this value with positive sign), w de = −1, using SNe Ia + H 0 + CC + BAO, from Table 3, Ref. [60].In Fig. (6a), we show the curve for MODEL 5 in the (q, r) plane.For early epochs, the curve starts in the lower right corner of the figure, in the region bounded by 0 < q < 1 and r > 0. For z ≫ 1, the parameters q, r and s yield {q, r, s} → {0.44, 0.84, 0.96}, and all three parameters values deviate significantly from the values expected for the standard matter-dominated era.Only the value for s is close to that of SCDM (s = 1).We observe that the trajectory is linear and evolves in the negative direction of q and positive of r.Before the present time, crosses the red straight line r = 1 corresponding to ΛCDM.At z = 0, we obtain {q, r, s} = {−0.56, 1, −0.019}.The parameter q is very close to the expected value for the current time, and r is exactly the value predicted by ΛCDM, while s is very close to zero.In the future, when we take the limit z → −1, the model gives {q, r, s} → {−1, 1.2, −0.037}.Then, the parameter r deviates significantly from ΛCDM.The parameter q is exactly the value predicted by ΛCDM and s is very close to zero.Fig. (6b) shows the (s, r) plane, and the curve starts in the region bounded by 0 < s < 1 and r > 0. The trajectory is non-linear and evolves in the negative direction of s and positive of r.We observe that before the current time, it crosses the fixed point (0, 1) corresponding to ΛCDM.The curve passes through ΛCDM but then evolves and away from it.We see that this model performs better near the present time since the values at that time are closer to ΛCDM than for other times.Then, we conclude this model is partially compatible with ΛCDM.
3. MODEL 6. Q6 = 3Hγ In order to evaluate this model numerically, we will use the following values Ω m0 = 0.326, γ = −0.08,w de = −1, using SNe Ia + H 0 + CC + BAO, from Table 3, Ref. [60].Fig. (7a) shows the evolution for MODEL 6 in the (q, r) plane.For early epochs, the trajectory starts in the upper right corner of the figure, in the region bounded by 0 < q < 1 and 0 < r ≤ 1.In this epoch, when z ≫ 1 the parameters q, r and s yield {q, r, s} → {0.5, 1, 0}, and the values of q, r and s agree with the matter-dominated standard era.The trajectory behaves much like MODEL 2, where the interaction also depends on the DE energy density.We see that in this case, the trajectory is curved and evolves in the negative q and negative r direction.At the past the curve moves away from the red straight line r = 1 corresponding to ΛCDM, and then moves forward.At z = 0, we obtain {q, r, s} = {−0.51,0.84, 0.05}.The parameter q is very close to the expected value for the present time, while r is far away to the value predicted by ΛCDM, and s is very close to zero.After, the trajectory moves away significantly from the dS fixed point, because in the limit z → −1, the model gives {q, r, s} → {−0.89, 0.69, 0.074}.Fig. (7b) shows the (s, r) plane, where the trajectory starts in the region bounded by 0 < s < 1 and 0 < r ≤ 1 and has non-linear behaviour, evolves in the positive direction of s and negative of r.We observe that it starts at the fixed point (0, 1) corresponding to ΛCDM, and then evolves to current time and diverges away.It is concluded that this model is partially compatible with ΛCDM.C. Models with a change of direction of energy transfer 1. MODEL 7. Q7 = q(α ρc + 3βHρc) For α = 0, it yields Q 7 = 3βqHρ c .We will evaluate this model numerically using the following values Ω c0 = 0.2738, Ω b0 = 0, β = −0.01,w de = −1, of TABLE I from Ref. [87].For MODEL 7 in the (q, r) plane, we observe from Fig. (8b) that the trajectory is non-linear and starts in the upper right corner, in the region bounded by 0 < q < 1 and 0 < r < 1.The curve evolves to the negative direction of q and negative of r.For z ≫ 1, the parameters q, r and s behaves as {q, r, s} → {0.51, 1, 1}, and q and r correspond to the standard matter-dominated era, and s has the value corresponding to SCDM.Before the present time, the trajectory crosses the red straight line r = 1 corresponding to ΛCDM and evolves.At z = 0, we obtain {q, r, s} = {−0.59,0.99, 0.002}, which are close to the values expected for ΛCDM.In the limit z → −1, the model gives {q, r, s} → {−1, 1, 0}, which correspond exactly with to the dS fixed point in the future.Thus, the model converges to ΛCDM.In the (s, r) plane, Fig. (8c) shows that the curve starts in the upper left corner, in the region bounded by 0 < s and r > 1.The trajectory is non-linear, passing through the point (0, 1) corresponding to ΛCDM before the present time.After the present time, the curve generates a twist and closes itself, converging to the point (0, 1).We corroborate the information provided by the (q, r) plane; the trajectory of the model converges to ΛCDM.As a consequence, we conclude that this model is compatible with ΛCDM, because the values of Statefinder parameters are very close to the values predicts for ΛCDM, in general, and in the future, the model evolves towards ΛCDM.
2. MODEL 8. Q8 = q(α ρtot + 3βHρtot) In order to evaluate this model numerically, we use the following values: Ω c0 = 0.2764, Ω b0 = 0, β = −0.0247,w de = −1, Table I from Ref. [87].Fig. (9a) shows the evolution for MODEL 8 in the (q, r) plane.For early epochs, the trajectory starts in the upper right corner of the figure, in the region bounded by 0 < q < 1 and 0 < r ≲ 1.1.In this epoch, when z ≫ 1 the Statefinder parameters q, r and s behave as {q, r, s} → {0.52, 1.1, 1}, and the values of q and r agree with the matter-dominated standard era, while s deviates significantly from the value s = 0 for ΛCDM.The value s = 1 corresponds to SCDM, being very similar to MODEL 7 in this epoch.We see that in this case, the trajectory is linear and evolves in the negative q and negative r direction.At the past, the curve moves away from the red straight line r = 1 corresponding to ΛCDM, and evolves.At z = 0, we obtain {q, r, s} = {−0.59,0.93, 0.02}, which are close to the values expected for ΛCDM.After, the trajectory moves away significantly from the dS fixed point, because in the limit z → −1, we obtain {q, r, s} → {−0.96, 0.89, 0.024}.9b) shows the (s, r) plane, where the trajectory starts in the region bounded by 0 < s < 1 and 0 < r ≲ 1.1 and has non-linear behaviour, evolves in the positive direction of s and negative of r.We observe that it starts at the fixed point (0, 1) corresponding to ΛCDM, and then evolves to the current time and diverges away.It is concluded that this model is partially compatible with ΛCDM.The values used for to evaluate this model numerically are: Ω c0 = 0.2717, Ω b0 = 0, β = 0.0136, w = −1, from TABLE I from Ref. [87].As in this Ref., the author used w de = −1, the "Λ" subscript that we use in the graphs is replacing in the equations by the subscript "de".We see in Fig. (10a) the curve for MODEL 9 in the (q, r) plane.The trajectory starts in the lower right corner of the figure, in the region bounded by 0 < q < 1 and r > 0. For z ≫ 1, the parameters q, r and s behaves as {q, r, s} → {0.5, 1.03, 0.007}, and all three parameter values are very close from the values expected for the standard matter-dominated era.We observe that the trajectory is non-linear and evolves in the negative direction of q and positive of r.Before the present time (z = 0), it crosses the red straight line r = 1 corresponding to ΛCDM.At the present, we obtain {q, r, s} = {−0.59,1.02, −0.008}.The parameters are close to the expected value predicted by ΛCDM for the present time.In the future, z → −1, the model gives {q, r, s} → {−1.02, 1.06, −0.014}.Then, the trajectory moves away significantly from the dS fixed point.Fig. (10b) shows the (s, r) plane, and the curve starts in the region bounded by 0 < s < 1 and r > 0. The trajectory is non-linear and evolves in the negative direction of s and positive of r.We observe that before present time, it crosses the fixed point (0, 1).The curve passes through ΛCDM, but evolves and away from it.We say that this model is partially compatible with ΛCDM.

D. Interaction involving derivatives
In order to evaluate this model numerically, we will use the values: h = 0.70, Ω c = 0.37, Ω b = 0.045, γ = 0.071, w de = −1.4,α = −0.15,extracted from Ref. [59], SNIa+H(z)+BAO+f gas , Table IV, Column B. We have preferred to use the values in column B, because in column C, although the CMB measurements are included, we have α = 0.0±0.1,so there is almost no interaction.In Fig. (11b) we see the curve for MODEL 10 in the (q, r) plane.For early times, the trajectory starts in the lower right corner of the figure, in the region bounded by 0 < q < 1 and r > 0. For z ≫ 1, the parameters q, r and s behave as {q, r, s} → {0.38, 0.68, 0.86}.All the three parameter values deviate significantly from the values expected for the standard matter-dominated era.The value of s = 0.86 is close to SCDM (r = 1, s = 1), but lower than the other  models that evolve in this way.We observe that the trajectory is linear and evolves in the negative direction of q and positive of r.Before the present time (z = 0) it crosses the red straight line r = 1 which corresponds to ΛCDM.At z = 0, we obtain {q, r, s} = {−0.73,2.35, −0.37}, and then all the three parameters values deviate significantly from the values expected for ΛCDM.In the future, for z → −1, it yields {q, r, s} → {−1.7, 4.1, −0.47}.Then, all the parameters deviate significantly from ΛCDM.We have used the values Ω m0 = 0.321, γ = −0.06,w de = −1, using SNe Ia + H 0 + CC + BAO, from Table 3 [60].This values are the same that for MODEL 4 to which this parameterization corresponds.show the parametric curves in the q − r and s − r planes, respectively.
the same, see Table (I).On the other hand, we have included the plot of the ratio between the DM and DE densities against z.We have corroborated what the authors of the reference [60] claim, that for γ < 0, the cosmic coincidence problem may be considerably alleviated.At early times, DM dominates over DE, if z ≫ 1 then r(z) → ∞.If z → −1 then r(z) → 0, thus, DE dominates over DM in the far future.For negative value of the coupling, the ratio r(z) yields the order of 1 earlier than in the ΛCDM model (γ = 0).

MODEL 1f (r) = r
This model is equivalent to MODEL 5, so we will use the following values: Ω m0 = 0.320, γ = −0.038,w de = −1, using SNe Ia + H 0 + CC + BAO, from Table 3, Ref. [60].In Fig. (14a), we see the curve for MODEL 13 in the (q, r) plane.The trajectory starts in the lower right corner of the figure, in the region bounded by 0 < q < 1 and 0 < r < 1.We calculate the parameters q, r and s at z ≫ 1, then {q, r, s} → {0.44, 0.84, 0.96}, and then all three parameter values deviate significantly from the values expected for the standard matter-dominated era.The value of s = 0.96 is very close to SCDM.We observe that the trajectory is non-linear and evolves in the negative direction of q and positive of r.At the present z = 0, we obtain {q, r, s} = {−0.52,0.98, 0.006}.The parameters are close to the values predicted by ΛCDM for the present time.After the present time it approaches to the red straight line r = 1 corresponding to ΛCDM.We verify this, since in the future, in the limit z → −1, the model gives {q, r, s} → {−1, 1, 0}.Then, the trajectory converges to the fixed dS point.Fig. (14b) shows the (s, r) plane, the curve starts in the region bounded by 0 < s < 1 and r > 0. The trajectory is non-linear and evolves in the negative direction of s and positive of r.We notice that at late times, it converges to the fixed point (0, 1) associated to ΛCDM.Hence, this model becomes partially compatible with ΛCDM.In this case, we have used the values: Ω m0 = 0.320, γ = −0.037,w de = −1, considering SNe Ia + H 0 + CC + BAO, from Table 3, Ref. [60].In the high redshift limit, z ≫ 1, the asymptotic values of q, r and s behave as {q, r, s} → {0.5, 1, 0.037}.For the present epoch z = 0, we obtain {q, r, s} = {−0.52,0.89, 0.037}.In the limit z → −1, the model gives {q, r, s} → {−0.94, 0.84, 0.037}.This model is equivalent to MODEL 2, and its behaviour in the (q, r) and (s, r) plane becomes the same, but as for the parameterized model we have used updated data.We can compare Figs.(3c) and (3d) with Figs.(15a) and (15b), respectively.We notice that the trajectory in the (s, r) plane evolves in the negative direction of r while holding fixed s = 0.037, but in MODEL 2 the trajectory gets lower by s = 0.02.This model is equivalent to MODEL 1, but, in this case, we have used the values Ω m0 = 0.309, γ = −0.0019,w de = −1, considering SNe Ia + H 0 + CC + BAO, from Table 3, Ref. [60].For early times, i.e. z ≫ 1, the parameters q, r and s behave as{q, r, s} → {0.5, 0.99, 1}.The value of s = 1 corresponds to SCDM.We recover the standard matter-dominated era due the values obtained for r and q.
At the present time, we get {q, r, s} = {−0.54,0.991, 0.001}.In the limit z → −1, the model gives {q, r, s} → {−1, 1, 0}.We compare Figs.(2b) and (2c) with Figs.(16a) and (16b), respectively.It is important to point out that for MODEL 1, the trajectories in the q − r and s − r planes do not converge to the point (−1, 1) and (0, 1), respectively.However, the trajectories for MODEL 15 in both planes do converge to these points.This model has similar behavior to MODEL 1, but for the parameterized model we have used updated data and w = −1.This may explain the convergence in the future to ΛCDM, since in MODEL 1 we have used w ̸ = −1.In this way, we say that this model is compatible with ΛCDM.Given the symmetry of MODEL 16, it is evaluated using the values extracted from Ref. [107], Ω m = Ω x = 0.11±0.03,and α = β = 0.25 ± 0.15.Furthermore, if we evaluate numerically the Symmetric and Asymmetric models using Ω m = 0.27, and setting α = β ≈ 0.001 (we cannot use α and β null), we get that the model is similar to ΛCDM.In Fig. (17a), we show the dimensionless Hubble rate against z, and we represent the self-interaction (dashed line) overlap to the Hubble rate for the ΛCDM model (solid line).We have calculated using the above values that for z < 0.56 both functions are almost identical, but they start to differ from z > 0.56.In the same reference, the authors show that this happens for z = 0.5.We find that for z ≫ 1, the parameters {q, r, s} → {0.88, 2.4, 1.2}, observing that all three values deviate significantly from the values expected for the standard matter-dominated era.Recalling that in the limiting case (Ω c + Ω b = 1, Ω de = 0), neglecting radiation, the value s predicted by the SCDM model is non-zero (r = 1, s = 1).In this case, s = 1.2 is far way from the value s = 1.From Fig. (17b), we observe that the trajectory In the (q, r) plane is linear and starts in the upper right corner, in the region bounded by 0 < q < 1 and r > 1.The curve evolves in the negative direction of q and negative of r.At the present, z = 0, one finds {q, r, s} = {−0.67,1.25, −0.071}, and also that q and r are not close values for the Statefinder parameters expected for ΛCDM, {q 0 , r 0 , s 0 } = {−0.55, 1, 0}.When z → −1, the parameters {q, r, s} → {−1, 1, 0}.Thus, the parameters present an asymptotic behaviour toward the dS fixed point in the future.In the (s, r) plane, Fig. (17c) shows that the curve starts in the upper right corner of the figure, in the region bounded by s < 0 and r > 1.The trajectory is non-linear in behavior and proceeds in the negative direction of r, but in the positive direction of s.We corroborate the information provided by the (q, r) plane; the trajectory of the model converges to ΛCDM.We conclude that this model is incompatible with ΛCDM, since, although it converges to ΛCDM, the present values and in the epoch that matter dominates, deviate significantly from ΛCDM.MODEL 17 is evaluated numerically using the values: Ω m = 0.14 ± 0.09, Ω x = 0.146 ± 0.085, α = 0.09 ± 0.17 and β = −0.01 ± 0.17.In Fig. (18a), we show the dimensionless Hubble rate, like the previous model, we plot the self-interaction (dashed line), which overlaps with the plot of Hubble rate for the ΛCDM model (solid line).Unlike the previous model, the H(z) curves of the asymmetric model and the reduced ΛCDM curve fit better for a larger redshift interval.We find that for z ≫ 1, the parameters {q, r, s} → {0.51, 1.02, 1}.It is observed that the values of q and r are very close to the expected values for the standard matter-dominated era.Additionally, the value of s = 1 corresponds to SCDM (r = 1, s = 1).In the (q, r) plane from Fig. (18b), we note that the trajectory is linear and starts in the upper right corner, in the region bounded by 0 < q < 1 and r < 1.The curve evolves in the negative direction of q and negative of r.At the current time, z = 0, it is found that {q, r, s} = {−0.6,1.01, −0.002}, and according, r and s are close values for the expected Statefinder parameters for ΛCDM, but the value of q(0) deviates significantly.When z → −1 the parameters {q, r, s} → {−1, 1, 0}.Thus, the parameters exhibit asymptotic behaviour towards the dS fixed point in the future, as also does the symmetric model.In the (s, r) plane, Fig. (18c) shows that the curve starts in the upper left corner of the figure, in the region bounded by s < 0 and r > 1.The trajectory has a non-linear behavior and advances in the negative direction of r, and in the positive direction of s.We can compare with the information provided by the (q, r) plane, the trajectory of the model converges to ΛCDM.We conclude that this model is compatible with ΛCDM.Although the value of q(0) deviates from ΛCDM, the values q, r, and s (as a whole) are close to ΛCDM in the all epochs.In contrast with the Symmetric model, we note that this model is able to reproduce a standard matter-dominated era.

VI. COMPARISON BETWEEN MODELS
Now, we are able to compare between the models by looking at the trajectories from where they start.We consider Figs.(19a) and (19b) for comparison.In accordance to [103] and [110], in the case of models 7, 16, and 17, it is observed that the pair (s, r) starts on the left-hand side of the fixed point (0, 1) corresponding to ΛCDM, which is characteristic of the hybrid expansion law (HEL), Chaplygin gas and Galileon models, such that s < 0 and r > 1 [110,111].We point out that this behaviour is different from the case of the Quintessence model, for which the trajectory in the (s, r) plane is observed to start in the region 0 < s < 1 and r < 1.On the other hand, we can see that the trajectory of MODEL 7 in Fig. (8b), is very similar to the HEL model, which the trajectory in the (q, r) plane starts in the region bounded by 0 < q < 1 and r > 1.The plot for model 10 in the (s, r) plane is very similar to that of models 5 and 9.For the latter model, the deviations are larger than in the similar models 5 and 9.As we established above, for the limiting case (Ω m = 1, Ω de = 0, and neglecting radiation), the values of the Statefinder parameters are r = 1, s = 1, corresponding to SCDM.It is observed in the Table (I) that when z ≫ 1, s becomes equal to 1 or very close to 1 obtained for models 1, 3, 5, 7, 8, 13, 15 and 17.We can conjecture and give an explanation for these values by assuming that it may be since in these models the Q function depends mainly on DM and because the impact of DE is smaller on the interaction, since the analysis takes place in the matter-dominated era.In addition, we can compare the MODEL 7 with a Statefinder analysis performed for running vacuum models (RVM).If we take the plots in r − s plane, we see that Fig. (8c) is similar to Fig. 1 of [80].As noted above, after the present time, the curve r − s generates a twist and closes on itself, converging to the point (0, 1).The similarity is that in our MODEL 7 and model I of [80], since they have a Q rate linearly depending on dark matter density.In contrast, for the DE energy density dependence (model II in [80]), we can see that the curve does not converge to the point (0, 1).In model II the asymptotic fixed point is a scaling solution, where the scaling solution is when the energy density of DE depends on an inverse power of the scale factor a(t).

VII. CONCLUSIONS
The nature of the components of the dark sector of our Universe is still unknown.Several models have been proposed in order to explain the accelerating expansion of the Universe.In this context, interacting DE models have been proposed to solve or alleviate the emerging cosmological tensions of ΛCDM model as well as the cosmic coincidence problem.In this paper, we have used the Statefinder diagnostic with the purpose to discriminate between several interacting DE models.In particular.we have investigated 17 interacting models between DE and DM, which have been already studied in the recent literature and constrained with astronomical and cosmological data.Specifically, the investigated models belong to the following categories: following categories: i) linear models in energy densities of DE and DM, ii) the current values of the (s, j) or (q, j) pair, while the dark dots on the curves reflect the matter-dominated phases of the models.Figure taken from [110].
Q q(z ≫ 1) r(z ≫ 1) s(z ≫ 1) q(0) r(0) s(0) q(−1) r(−1) s(−   non-linear models, iii) models with a change of direction of energy transfer between DE and DM, iv) models involving derivatives of the energy densities, v) parametrized interactions through a function of the coincidence parameter r, and finally we also consider two kinds of models with a self-interaction between DM, without DE.The models we have chosen allow us to solve the conservation differential equations and to find an analytical expression for the dimensionless Hubble rate E(z) = H(z)/H 0 .Thus, by using E(z) we have computed the expressions for the Statefinder parameters q(z), r(z), and s(z).Then, these parameters were evaluated for the matter-dominated era (z ≫ 1), at the present epoch (z = 0), and in the future (z → −1).These values are shown in Table (I), as a summary.
It was necessary to plot in the q − r and s − r planes the respective parametric curves for each case.The trajectory of each model was compared with the points associated to ΛCDM, i.e., (0, 1) in the s − r plane and the dS fixed point (−1, 1) in the q − r plane.Additionally, we found that the "anomalous" behavior of the parameter s for z ≫ 1, corresponds to the point (1, 1), associated to the SCDM model.Differences and similarities were found between models of the same category (linear, non-linear, etc.), depending on the deviations obtained with respect to ΛCDM.Comparing all of these models, some of them converge to ΛCDM at late times, such as models 4, 7, 11, 13, 15 and the self-interacting models 16 and 17.In that sense, the worst performing model in the future is MODEL 10, since its deviations are very large with respect to the others.
In the current epoch, the models 5, 3 and 15 show which are values very close to those predicted by ΛCDM, i.e., {q 0 , r 0 , s 0 } = {−0.55, 1, 0}.But, to a lesser extent, models 13, 7, 9 and 17 are close to ΛCDM.On the other hand, the worst behavior corresponds to MODEL 10, and the symmetric model, since these deviate significantly for the current epoch, especially the parameter q.
In models 2, 4, 6, 9, 11, 12 and 14 the values found for the three Statefinder parameters behave very similar to ΛCDM in the matter-dominated era, i.e., {q(z ≫ 1), r(z ≫ 1), s(z ≫ 1)} = {0.5, 1, 0}.At this same epoch, models 1, 3, 5, 7, 8, 13, 15, and the asymmetric model predict values of s close to 1, which we established is related to SCDM.We inferred that this was due to the larger impact of DM in the interactions, given the form of the Q rate and the small contribution of the DE energy density.Furthermore, at this epoch, we see that models 16, 10, 1, 3, 5, and 13 deviate significantly from the q and r values predicted by ΛCDM.
The models parameterized and their equivalents have similar behaviours, except 5-13 and 1-15.Actually, we found that models 5 and 1 do not converge to ΛCDM in the future, but 13 and 15 do instead.As we said, this may be due to the fact that the parameterized models use more updated data and w = −1.This leads us to conclude that the models analyzed by Statefinder analysis cannot be definitively discarded, as more updated data can lead to better model performance.
As we have shown, we can compare our models with other different DE models extensively studied in the literature, such as HEL, Chaplygin gas, Galileon models, Quintessence , and DGP from [110].This can be done by considering the deviations from ΛCDM and the shape of the trajectories in the (s, r) plane, as we look at the origins of the curves.
For future research, we would like to complement this work and to corroborate the convergence of the models by analyzing dynamical systems and observing the stability of the fixed points, as has been done for other DE models [80,112].This work opens up the possibility to further investigation on others interacting DE models that have a more complex expression for Q, for example, where the EoS parameter w or strengths of the couplings depend on redshift or time.We hope to be able to address this point in a future work.

FIG. 2 :
FIG. 2: In Figure (a) it is shown the Hubble rate H(z) against z, for the model with the coupling function Q1 = 3γHρc,

FIG. 3 :
FIG. 3: In this figure we show the plots of the Statefinder parameters for the coupling function Q2 = 3γHρ de .Figs.(a) and (b) show the evolution of r(z) and q(z), respectively.We represent the interacting model DE-DM (black dashed line), which overlaps with ΛCDM (blue solid line).For this model, the figures (c) and (d) show the q − r and s − r planes, respectively.

FIG. 5 :
FIG. 5: Fig.(a) shows the Hubble rate H(z) as a function of redshift z for the model Q4 = 3Hγ ρcρ de ρ .We have used H0 = 69.44 km s −1 Mpc −1 , from Ref. [60].In Figures (b) and (c) we show the parametric curves in the q − r and s − r planes, respectively.

FIG. 6 : 2 cρ.
FIG. 6: Left panel: In this figure, we have plotted the parametric curve in the q − r plane.It shows the evolution for the model with the coupling function Q5 = 3Hγ ρ 2 c ρ .Right panel: The curve in the s − r plane it is shown.

FIG. 7 :
FIG. 7: Left panel: This figure shows the evolution in the q − r plane for the model with the coupling function Q6 = 3Hγ ρ 2 de ρ .Right panel: The graph shows the curve for this model in the s − r plane.

FIG. 8 :FIG. 9 :
FIG. 8: Figure (a) shows the dimensionless Hubble rate H(z)/H0 as a function of redshift z for the model with the coupling function Q7 = 3βqHρc.We represent the interaction with a dashed line, and also the Hubble rate of the ΛCDM model with the solid line.For this model, Figures (b) and (c) show the curves in the q − and s − r planes, respectively.

Fig. (
Fig.(9b) shows the (s, r) plane, where the trajectory starts in the region bounded by 0 < s < 1 and 0 < r ≲ 1.1 and has non-linear behaviour, evolves in the positive direction of s and negative of r.We observe that it starts at the fixed point (0, 1) corresponding to ΛCDM, and then evolves to the current time and diverges away.It is concluded that this model is partially compatible with ΛCDM.

3. MODEL 9 .FIG. 10 :
FIG. 10: Figures (a) and (b) show the evolution for the model with the coupling function Q9 = 3βqHρ de .Left panel: Figure shows the q − r parametric curve.Right panel: Figure shows the s − r parametric curve.
Fig.(11c)shows the s − r plane, and the curve starts in the region bounded by 0 < s < 1 and r > 0. The trajectory is non-linear and evolves in the negative direction of s and positive of r.The plot in the s − r plane shows that in the first place, the curve passes through ΛCDM.After, it passes through the point associated with the present time, and then it deviates significantly in the future.We say that this model is incompatible with ΛCDM, since the Statefinder values deviate significantly from the expected values in all three epochs.E. Parameterized interactions 1. MODEL 11. f (r) = 1

f
FIG. 12: In Figure (a) it is shown the Hubble rate H(z) as a function of redshift z for the model with function f (r) = 1.We have used H0 = 69.44 km s −1 Mpc −1 , from [60].For the same model, Figure (b) shows the evolution of the ratio r, where we have represented the interacting model (black dashed line), which overlaps with ΛCDM (blue solid line).Figures (c) and (d)show the parametric curves in the q − r and s − r planes, respectively.

FIG. 13 :f
FIG. 13: In the figures we show the evolution for the model with function f (r) = 1 r .Left panel: Figure shows the q − r parametric curve.Right panel: Figure shows the s − r parametric curve.

FIG. 15 :
FIG. 15:For the model with function f (r) = 1 + 1 r , the parametric curves in the q − r and s − r planes are shown.

ffFIG. 16 :
FIG. 16:We show the parametric curves in the q − r and s − r planes for the model with parametric function f (r) = 1 + r.

3 3 3
FIG. 17: In Figure (a), it is shown the dimensionless Hubble rate H(z)/H0 as a function of redshift z for the Symmetric model, Qm = −3Hαρm, Qx = −3Hβρx.It is represented the self-interaction (dashed line), as also, the Hubble rate of ΛCDM model (solid line).For this model, Figures (b) and (c) show the parametric curves in the q − r and s − r planes, respectively.

Q m = - 3 3 3 FIG. 18 :
FIG. 18: In Figure (a), it is shown the dimensionless Hubble rate H(z)/H0 as a function of redshift z for the Asymmetric model, Qm = −3Hαρm, Qx = −3Hβ(ρm + ρx), and we represent the self-interaction (dashed line), as also, the Hubble rate of ΛCDM model (solid line).For this model, Figures (b) and (c) show the parametric curves in the q − r and s − r planes, respectively.

FIG. 19 :
FIG. 19: The parametric curves in the q − j plane are shown in Figure (a).The dS state q = −1 is shown by the vertical purple line.The plots in the s − j plane are displayed in Figure (b).The HEL model is shown by the green curve in both panels.The ΛCDM, DGP, Chaplygin gas, and Galileon models are represented by the Red, Cyan, Magenta, and Blue curves, respectively.The ΛCDM point (0, 1) is where horizontal and vertical dashed lines connect.The dots on the curves representthe current values of the (s, j) or (q, j) pair, while the dark dots on the curves reflect the matter-dominated phases of the models.Figure taken from[110].

TABLE I :
Summary of present and asymptotic values of Statefinder parameters for all the interacting models studied in the present research.