Dynamics and statefinder analysis of a class of sign-changeable interacting dark energy scenarios

We revise the dynamical properties of a class of cosmological models where the dark sector interacts through an interacting term that changes sign during evolution. In particular, we obtain the critical points and we investigate the existence and stability conditions for cosmological solutions, describing radiation, matter and dark energy dominated eras. We find that all the studied models admit a stable critical point corresponding to an accelerated phase. We use background data to find the best fit parameters for one of the studied models, resulting an interacting parameter with a definite sign within $1\sigma$ confidence level, consistent with the results of the dynamical system analysis. We also compute the statefinder parameters and plot the $r-q$ and $r-s$ planes, where we observe different trajectories when we vary the interaction parameter for a specific model and when we vary the interacting scenario. We can in this sense distinguish among models, including $\Lambda$CDM.


Introduction
All the available data in cosmology seems to indicate that our current universe consist of 69% of dark energy and 26% of dark matter [1], but to the present understanding of physics, both components have unknown nature [2,3]. Usually we consider the standard cosmological scenario ΛCDM, where the dark components evolve independently. A different approach is to consider that both components interact [4,5]. Interacting scenarios allow to alleviate the coincidence problem [4] and also are useful in addressing the tension in the H 0 parameter [6]. Among interacting scenarios it is common to consider those models in which the interaction sign remains constant during evolution but there are some less explored scenarios where the interaction sign can change [7,8,9]. In cosmology, dynamical system methods are a very useful technique to investigate the asymptotic behavior of a given cosmological model, specially when the analytical solutions are unknown. In Ref. [10] the authors review the dynamics of a plethora of dark energy models applying dynamical system analysis, an updated version of this work is presented in Ref. [11]. Both works include interacting scenarios. The authors of Ref. [12] introduce the cosmological pair r and s, dubbed statefinder parameters, depending on derivatives of the scale factor up to third order, where the evolution of these parameters a e-mail: fabiola.arevalo@umayor.cl b e-mail: acidm@ubiobio.cl allows to differentiate among dark energy models. For instance, in Ref. [13] the authors use the statefinder parameters to characterize interacting quintessence models of dark energy. In Refs. [14,15] the authors apply dynamical systems methods and perform a statefinder analysis to interacting dark energy models. Even though the dynamical system analysis shows a similar asymptotic behavior for these models, the statefinder diagnosis allows to distinguish among models, when the evolution of the statefinder parameters is studied. The goal of this work is to revise the dynamics of a class of sign-changeable scenarios which do not have analytical solutions, then it is interesting to characterize the asymptotic behavior of these models through dynamical system analysis and study its performance in fitting the known asymptotic behavior of our universe. The studied models are then compared through the evolution of the statefinder parameters to distinguish its behavior from ΛCDM. This work is organized as follows. In section 2 we describe interacting scenarios and the class of signchangeable interactions to study. Section 3 is devoted to the dynamical system analysis of each of the models described in section 2. In section 4 we use background data to perform the observational contrast for one of the studied models. In section 5 we compute the statefinder parameters and we study the features of the trajectories in the r − q and r − s planes. Finally, in section 6 we present our final remarks on this work.

Sign-changeable interacting dark energy models
Let us consider a flat, homogeneous and isotropic universe in the framework of General Relativity and a cosmological scenario containing radiation (r), baryons (b), cold dark matter (c) and a dark energy (x). In the context of an interacting dark energy scenario the conservation equation can be separated into the following equations: where by convenience we use η = ln a as time variable and we define ( ) = d/dη (where a is the scale factor). Notice that Γ > 0 indicates an energy transfer from dark matter to dark energy and Γ < 0 indicates the opposite. In this work we study the following sign-changeable linear interactions, where α 1i and α 2j are constants, q is the deceleration parameter defined as q = −1 − H /H and H is the Hubble expansion rate. The subscripts T and d denote total energy density and dark sector energy density, (ρ c + ρ x ), respectively. It was noticed in Ref. [16] that interaction Γ 1T has analytical solutions for the dark sector energy densities, however the other scenarios do not, to the authors knowledge. We do not include Γ 2x in our analysis given that this model leads to a constant deceleration parameter in an interacting scenario. Notice that all the studied models reduces to the ΛCDM scenario when the interacting parameter α 1i or α 2j is set to zero.

Dynamical System Analysis
In this section we apply dynamical system methods [10,11] in order to identify the relevant cosmological eras in the studied models. In this sense, we rewrite the system of Eqs. (1)-(4) in terms of the density parameters Ω i = ρ i /3H 2 (dimensionless) and we use the Friedmann equation as a constraint among parameters, i.e., to reduce the system of equations (1)-(4) to: where interactions (5)-(6) are given in Table 1 in terms of the density parameters. The stability of the relevant cosmological epochs is found by calculating the eigenvalues of the linearized system at the critical points. For each cosmological interaction we follow the following scheme. First we consider the set of equations, where f i is a function of the density parameters Ω r , Ω c and Ω x and i = r, c, x. From this we find the critical points Ω * l of the system (11) by calculating Then we linearized the set of equations (11) around the critical points, where J l i = ∂f i ∂Ω l is the jacobian matrix. By analysing the jacobian matrix we can find regions in the parameter space where a specific behavior is presented. In particular, we have stable, saddle or unstable points when the real part of the eigenvalues are, all negatives, a mixture of positives and negatives or all positives, respectively.
In the following, we describe in detail each of the models in Table 1. First we present models Γ 1T , Γ 1x and Γ 2T which can be reduced to two-dimensional systems. Subsequently we describe models Γ 1d , Γ 1c , Γ 2d and Γ 2c , which are analysed in terms of three-dimensional systems. We describe each of the critical points with the focus on radiation, matter and dark energy dominance, corresponding to an unstable, saddle and stable points, respectively.

Case A: Γ 1T
For the model Γ 1T the set of Eqs.(8)-(10) becomes, We notice that the dynamical system analysis can be simplified given that Eqs. (12) and (14) become a closed system, then we can perform a 2-dimensional analysis. By considering Ω r = 0 and Ω x = 0 simultaneously we find the critical points of the system. In order to obtain the stability of the critical points we linearized the system around the critical points and we analyzed the corresponding eigenvalues.
For the model Γ 1T we find the following three critical points A i = {Ω r , Ω x }, where the matter contribution to the critical points Ω m = Ω b + Ω c is determined by the constraint In the physical analysis of these points we have considered as existence conditions Ω r ≥ 0 and Ω m ≥ 0 at the critical points. On the other hand, notice that the effective state parameter ω eff = p T ρ T (p T is the total pressure) is related to the deceleration parameter by q = 1 2 (1 + 3ω eff ), and in terms of the density parameters q is given by Table 2 shows the existence and stability conditions on the eigenvalues A 1 − A 3 , as well as the ω eff value. The critical point A 1 in Table 2 describes a radiation dominated epoch where radiation, matter and dark energy coexist for 0 < α < 4 9 . The matter contribution to this critical point is given by Ω m = 3α and in this case Ω x is always negative given the existence conditions. The signs of the eigenvalues of the linearized system indicate that A 1 is an unstable critical point for α < 4 9 . The critical point A 2 has ω eff > − 1 3 irrespective of the value of α. In the limit α → 0, this phase can be understood as a phase of matter domination, given that ω eff → 0 in this limit. In this phase matter and dark energy coexists. At this critical point the matter component is Ω m =  Table 2 Description of the critical points for scenario Γ 1T .
, which is positive irrespective of the α value. On the other hand, we have Ω x < 0 for α > 0 and Ω x > 0 for α < 0. A 2 is a saddle critical point for α < 4 9 or an unstable critical point for α > 4 9 . The critical point A 3 dynamically corresponds to a phase of dark energy domination (ω eff < − 1 3 irrespective of the α value), where matter and dark energy coexist. At this point Ω m = 1 , which is positive for α > 0. Under this existence condition Ω x is always positive. A 3 is a stable critical point for α > 0.  Figure 1 shows the 2-dimensional phase space Ω x − Ω r . The three critical points in Table 2 are shown as red dots, where we can clearly notice the unstable, saddle and stable points, these 3 points can be find in model Γ 1T for 0 ≤ α < 4 9 , nevertheless, in this range Ω x is negative at points A 1 and A 2 . We also show a possible trajectory going through the three critical points, starting at radiation domination and ending at dark energy domination. In the figure, the "Present" point is indicated considering {Ω m , Ω x } = {0.3111, 0.6889} [1]. Finally notice that the critical points A 1 − A 3 in table 2 and their existence and stability conditions reduce to the ones in the ΛCDM scenario in the limit α → 0. Table 3 Description of the critical points for scenario Γ 1x .

Case B: Γ 1x
For the model Γ 1x we can reduce the set of Eqs.(8)-(10) to a 2-dimensional set including Eq.(8) and the following equation, For this model we find three critical points B i = {Ω r , Ω x }, where the matter contribution to the critical points is determined by the constraint (15). The existence and stability conditions of the critical points B 1 − B 3 are shown in Table 3, as well as the ω eff parameter. The critical points B 1 and B 2 in Table 3 corresponds to a phase of radiation and matter domination, respectively. There are no conditions on the existence of these critical points and the stability analysis indicate that B 1 is unstable for α > − 4 3 and saddle for α < − 4 3 , meanwhile, B 2 is saddle for α > −2 and stable for α < −2.
At the critical point B 3 matter and dark energy coexist. At this point Ω m = 2α 2+3α which is positive for α < − 2 3 or α > 0. On the other hand, for α > 0, this point corresponds to dark energy domination with Ω x > 0 as a stable critical point. For α < − 2 3 we get ω eff < − 1 3 , Ω m > Ω x and Ω x > 0 for α < −2. The corresponding stability conditions in this case are given in Table 3. The figure 2 shows the 2-dimensional phase space Ω x − Ω r for the Γ 1x scenario, where the three critical points in Table 3 are shown as red dots and we can clearly notice the unstable, saddle and stable points. These points can be find in model Γ 1x for α ≥ 0, where Ω x is positive at the critical point B 3 . In the figure the black dot labeled as "Present" considers the values in Ref. [1]. As well as model Γ 1T , the critical points B 1 − B 3 reduce to the ΛCDM critical points in the limit α → 0.

Case C: Γ 2T
For the model Γ 2T we can reduce the set of Eqs.(8)-(10) to a 2-dimensional set including Eq.(8) and the following equation, For this model we find three critical points C i = {Ω r , Ω x }, where the matter contribution to the critical points is determined by the constraint (15).
The existence and stability conditions of the critical points are shown in Table 4, as well as the ω eff parameter.
At the critical point C 1 coexist radiation, matter and dark energy in a radiation dominated era (ω eff = 1 3 ). For this critical point we have Ω m = −12α which is positive for α < 0 but necessarily in this case Ω x < 0. At the critical point C 2 coexist matter and dark energy, ω eff ≥ − 1 3 for α > − 2 9 . At this point we have Ω m = 2+6α 2+9α which is positive for α < − 1 3 or α > − 2 9 . In the range α < − 2 9 or α > 0, the parameter Ω x is also positive. This point corresponds to matter domination in the limit α → 0. The point C 3 corresponds to a de-Sitter attractor for α > − 1 3 , where only the dark energy term contributes.
The figure 3 shows the 2-dimensional phase space Ω x − Ω r for the Γ 2T scenario. The three critical points in Table 4 are shown as red dots, we can clearly notice the unstable, saddle and stable points, which can be find in the model Γ 2T for − 1 9 < α ≤ 0. However, in this range, Ω x is negative at the critical points C 1 and C 2 . In the figure the "Present" dot considers the values in Ref. [1].
The critical points C 1 − C 3 corresponding to model Γ 2T in Table 4 reduce to the ΛCDM critical points in the limit α → 0.  Table 4 Description of the critical points for scenario Γ 2T .

Case D: Γ 1d
For the model Γ 1d the set of Eqs.(8)-(10) reduces to, , , For this model we find the following critical points where the Ω c contribution to the critical points is determined by the constraint (7).
In the physical analysis of these points we have considered as existence conditions Ω r ≥ 0, Ω b ≥ 0 and Ω c ≥ 0 at the critical points. The critical points D 1 and D 2 in Table 5 describe a radiation and a baryon dominated epoch, respectively, irrespective of the α value. The signs of the eigenvalues of the linearized system indicate that D 1 is an unstable critical point in the range − 1 4 ≤ α < 4 9 . On the other hand, the critical point D 2 corresponds to a saddle point in the above range but it represents a non-physical point because there is not a baryon dominated epoch in the actual evolution of our universe. We notice that scenarios Γ 1T and Γ 1d share the same critical points for the dark sector, D 3 and D 4 , with the same existence and stability conditions. The critical point D 3 can be understood as dark matter domination in the limit α → 0, this point (and its analysis) is equivalent to the point A 2 for the scenario Γ 1T . The critical point D 4 corresponds to dark energy domination for α > 0. This point and its analysis is equivalent to the points A 3 in scenario Γ 1T . The existence conditions in Table 5 indicate that these critical points (and the radiation one) can be find for 0 ≤ α < 4 9 , but in this range necessarily Ω x becomes negative at the critical point D 3 . In Fig.4 we show a projected 2-dimensional phase space Ω x − Ω c where we have fixed Ω r = 0. In this figure we show two critical points (physical) corresponding to the saddle point D 3 and the dark energy attractor D 4 .

Case E: Γ 1c
For the model Γ 1c the set of Eqs.(8)-(10) reduces to, For this model we find the following critical points where the contribution Ω c to the critical points is determined by the constraint (7) to be zero for each point but E 3 , where Ω c = 2(1−α) 2−3α . Critical Points ω eff Existence Stability Conditions stable for α > 1 saddle for α < 1 3 or 2 3 < α < 1 E 4 -1 α ∈ R stable for α < 1 saddle for α > 1 Table 6 Description of the critical points for scenario Γ 1c .
The critical points E 1 and E 2 in Table 6 describe a radiation and a baryon dominated epoch, respectively. The signs of the eigenvalues of the linearized system indicate that E 1 is an unstable critical point for α < 1 3 . On the other hand, the non-physical point E 2 (baryon domination) corresponds to a saddle point in the above range. The point E 3 represent a phase where the dark fluids coexist with Ω c > 0 for α < 2 3 or α > 1 and ω eff > − 1 3 for α < 2 3 , this point corresponds to dark matter domination in the limit α → 0. The point E 4 corresponds to a de-Sitter stage with the dominance of dark energy. This point is stable for α < 1.
From Table 6 we notice that the critical points corresponding to unstable, saddle and stable behavior can be find for α < 1 3 . Considering this range, Ω x remains positive for α < 0 at the critical point E 3 . In , For this model we find the following critical points F i = {Ω r , Ω b , Ω x }, where the contribution Ω c to the critical points is determined by the constraint (7) to be zero at each point but F 3 , where Ω c = 2+6α 2+9α . The critical points F 1 and F 2 in Table 7 describe a radiation and a baryon dominated era, respectively. The signs of the eigenvalues of the linearized system indicate that F 1 is an unstable critical point for α > − 1 9 . On the other hand, the non-physical point F 2 (baryon domination) corresponds to a saddle point in the above range. The point F 3 represents a phase where the dark fluids coexist with Ω c > 0 for α < − 1 3 or α > − 2 9 and ω eff > − 1 3 for α > − 2 9 , this point corresponds to dark matter domination in the limit α → 0. The point F 3 is saddle in the range α > − 1 9 . The point F 4 corresponds to a de-Sitter stage with the dominance of dark energy and this point is stable for α > − 1 3 . From Table 7 we notice that the critical points corresponding to unstable, saddle and stable can be find for α > − 1 9 . Notice that in this range, Ω x remains positive for α > 0 at the critical point F 3 . In Fig.6 we show the projected 2-dimensional phase space Ω x − Ω c where we have fixed Ω r = 0 and we have chosen α = 0.01. In this figure we show two critical points (physical) corresponding to the saddle point F 3 and the attractor point F 4 . Notice that the models Γ 2T and Γ 2d share the same critical points and existence/stability conditions in the points corresponding to the dark sector, C 2 − C 3 and saddle for 2 + 9α = 0 Table 7 Description of the critical points for scenario Γ 2d . , , For this model we find the following critical points G i = {Ω r , Ω b , Ω x }, where given Eq. (7), Ω c is zero at G 1 and G 2 , at G 4 . The critical points G 1 and G 2 in Table 8 describe a radiation and a baryon dominated epoch, respectively. The signs of the eigenvalues of the linearized system indicate that G 1 is an unstable critical point for α < − 1 3 or α > − 1 12 . The non-physical point G 2 (baryon domination) corresponds to a saddle point in the above range. At the critical point G 3 the dark sector coexist in the ranges α < − √ 3+2 3 , < α < 0 or α > 0. Inside these limits G 3 can be unstable for and this point reduces to dark matter domination in the limit 3 . At the point G 4 coexist the dark components in the range where it is a stable critical point. In the range where G 4 is a saddle point we have ω eff > − 1 3 , outside this range, in the allowed region, the effective fluid accelerate the universe's expansion.
The Ω x contribution of this point is positive 3 . Notice that the point G 4 diverges in the limit α → 0. The point G 5 corresponds to a de-Sitter stage with the dominance of dark energy. This point is stable for α < 1 3 and unstable outside this range. We have sketched the behavior of the critical points G 3 − G 5 in Fig. 8 < α < − 1 12 there is one unstable, one saddle and one stable point, G 3 , G 4 and G 5 , respectively. In the range − 1 12 < α < 0 we have two saddle points, G 3 and G 4 , and one stable point, G 5 . For 0 < α < 1 3 we have one saddle point and one stable point, G 3 and G 5 , respectively. For α > 1 3 we have two saddle points (G 3 and G 5 ) and no stable point. In Fig. 7 Fig. 7 Projected phase plots for model Γ 2c with α = 0.01 In summary, if we search for the models having an unstable, saddle and stable critical points during the evolution we get the following intervals, where existence/stability conditions are jointly met: α > 0 for Γ 1x , α < 0 for Γ 1c , α > 0 for Γ 2d and 0 < α < 1 3 for Γ 2c . For these models Ω x > 0 and the sign of the interacting parameter α becomes defined in each case.
On the other hand, the corresponding intervals where existence/stability conditions are jointly met for models where Ω x < 0 are: 0 < α < 4 9 for Γ 1T , − 1 9 < α < 0 for Γ 2T and 0 < α < 4 9 for Γ 1d . Models Γ 1T and Γ 2T have Ω x < 0 at radiation and dark matter domination, nevertheless Ω x changes sign and it becomes positive at dark energy domination. For model Γ 1d we have Ω x < 0 for dark matter domination, but it is positive for radiation and dark energy domination, meaning at least two change of sign of Ω x during evolution.
Finally, notice that the stability conditions of center type critical points have not been addressed in this work.

Observational analysis
In the observational analysis we use background data such as the local measurement of Hubble parameter (H 0 ) [18], cosmic chronometers [19], supernovae type Ia (SNe Ia) [20], baryon acoustic oscillations (BAO) [21], and the angular scale of the sound horizon at the last scattering [22]. In this section we briefly describe the used data, a detailed description can be find in Ref. [17].

Local determination of H 0
We use the H 0 value obtained by the SH0ES team using a local distance ladder method based on Cepheids, h = 0.7402 ± 0.0142 [18].

Cosmic chronometers
We use 24 cosmic chronometers obtained through the differential age method [23]. This procedure provides cosmological-independent direct measurements of the expansion history of the universe up to redshift 1.2 [24], see Table 3 in Ref. [17].

Supernovae Type Ia
We use the Pantheon Sample, a set of 1048 spectroscopically confirmed SNe Ia ranging from redshift 0.01 to 2.3 [20], along with the corresponding covariance matrix. The Pantheon catalog contains measurements of peak magnitudes in the B-band's rest frame, m B , related to the distance modulus by µ = m B + M B . M B is the absolute B-band magnitude of a fiducial SN Ia, a nuisance parameter. In our analysis the distance modulus at a given redshift is, where the luminosity distance d L is given by,

BAO data
We use isotropic measurements of the BAO signal from 6dFGS [25], MGS [26] and eBOSS [27]. These measurements are given in terms of the dimensionless ratio, where z d is the redshift at the drag epoch, with c the speed of light and R = 3Ω b 4Ω γ (1+z) [28]. Furthermore, we use the anisotropic BAO measurements from BOSS DR12 [29] and Lyα forest [30], which are defined in terms of D A and D H = c/H(z), as shown in table 2 of Ref. [31]. We use these data along with the corresponding covariance matrix in Ref. [31].

Cosmic microwave background data
We consider the angular scale of the sound horizon at the last scattering as the only contribution of CMB data, where the comoving size of the sound horizon is evaluated at z * = 1089.80 [32]. We compare the value obtained in our study with the one reported by the Planck collaboration in 2015, a = 301.63±0.15 [22].
To compute the maximum likelihood and posterior distributions we use the MultiNest algorithm 1 [33,34], with a global log-evidence tolerance of 0.01 as a convergence criterion and working with 800 live points to improve the accuracy.
We perform the observational analysis for the only model with an analytical solution, i.e. scenario Γ 1T . The posterior distribution for the parameters h, Ω c and α are shown in Fig. 9 and Table 9. In this analysis we fix the physical baryon density to Ω b h 2 = 0.022383 [32], the physical photon density to Ω γ h 2 = 2.469 × 10 −5 [36], we consider the the radiation density as Ω r = Ω γ 1 + 7 8 4 11 4/3 N eff , and the effective number of neutrinos N eff = 3.046 [35]. Notice that we have considered a flat prior for the α parameter consistent with the results of section 3.1, i.e., 0 < α < 4 9 .  Table 9 Best fit parameters and 1σ error for the observational analysis of the Γ 1T scenario. Fig. 10 Evolution of the density parameters for the model Γ 1T considering the best fit parameters in Table 9.
The results in Table 9 are consistent with a small but non-null interaction in the dark sector. Given that the interaction Γ 1T changes sign during evolution, our results indicate that today dark energy is transferred to dark matter but in the past the transference was the opposite. In Fig.10 we show the evolution of the density parameters consistent with the best fit parameters shown in Table 9. Here we notice that the dark energy density was negative in the past and around z = 7.2 it transitioned to a positive dark energy density, as shown in the subplot in Fig. 10. In the limit z → ∞, the dark energy density tends to Ω x → −0.00675 instead of 0, as in the ΛCDM model. Finally, notice that a negative dark energy density in the past is a intrinsic feature of model Γ 1T , which arise from imposing the existence conditions in the dynamical system analysis.

Statefinder parameters
The statefinder parameters r and s are introduced in order to compare the plethora of dark energy models available. In the definition of these parameters, derivatives of the scale factor exceed the second-order. We can explore the phase space by defining: In terms of the density parameters r and s are given by: where q is defined in (16). For each model in Table 1 we can obtain the parameters r and s and the r − q and r − s planes can be drawn. In order to do this we have to numerically integrate equations (8)-(10) and we have considered as border conditions the present values of the density parameters as {Ω m , Ω c , Ω x } = {0.3111, 0.26212, 0.6889} [1]. The trajectories in the planes r − q and r − s can exhibit quite different behaviors, even though the overall behavior of models may appear similar in the dynamical system analysis of the previous section. The r − q and r − s planes for the Γ 1T model are shown in Fig. 11, here the trajectories correspond to different values of α consistent with the results of section 3.1. We notice that in both planes the variation in α modify the the trajectories with respect to the ΛCDM model. In Fig. 12 we show the r−q and r−s planes for the models Γ 1T , Γ 1d and Γ 1x in the case α = 0.01. From the trajectories we see that all the curves pass through the points (r, q) = (1, 0) and (r, s) = (1, 0), and they present the same outcomes in the future. Notice that models Γ 1T and Γ 1d are very similar in both planes. Fig. 13 shows the r − q and r − s planes for models Γ 1c , Γ 2T , Γ 2d and Γ 2c , we consider α = −0.01 for models Γ 2T and Γ 2d and α = 0.01 for models Γ 2d and Γ 2c , consistent with the results in section 3. We see that all trajectories in Figure 13 pass through the points (r, q) = (1, 0) and they end at the points (r, q) = (1, −1) and (r, s) = (1, 0). Even though the interacting models trajectories end in the same point as ΛCDM in the plane r − s, the interacting models present a loop in their trajectories, passing through the point (r, s) = (1, 0) twice.  Fig. 13 The figure shows the planes r − q (left) and r − s (right) for the models Γ 2T and Γ 1c in the case α = −0.01 and models Γ 2c and Γ 2d in the case α = 0.01. Black dots indicate current values and gray arrows the evolution's direction. The solid black line represents ΛCDM, meanwhile, purple, red, orange and black dotted lines correspond to models Γ 2T , Γ 1c , Γ 2d and Γ 2c , respectively.
We notice that the ending point is different from ΛCDM for Figures 11 and 12, but is the same point for Fig. 13. The trajectories in Fig. 12 tend to an asymptotic point corresponding to a scaling solution and in Fig. 13 the trajectories evolve towards the de Sitter point. Notice that the trajectories in Fig.  13 have different shapes from those shown in Fig. 12, this separation seems not to be related with the classification of models in (5)- (6).
Finally, similar trajectories to those shown in Figures 12 and 13 are found in Ref. [14], where two different interacting models where studied, which do not have a change of sign in the interacting term.

Final Remarks
We have performed a dynamical system analysis and a statefinder analysis for a class of sign-changeable interactions in the dark sector [8]. We have shown that the studied models can explain the current acceleration of our universe as a late time attractor in a defined range for the α parameter (see section 3). For the studied models the evolution of the universe can go from a radiation dominated era, through a dark matter dominated era and a final stage of dark energy domination. These kind of evolution is found in each model for a defined range of the interaction parameter α (see section 3). We find that models, Γ 1T , Γ 2T , Γ 1d , need to have Ω x < 0 at some critical points in order to follow the universe's evolution, this necessarily implies at least one change of sign of Ω x during evolution, resulting in Ω x > 0 at the final attractor. On the other hand, the requirement of the existence and stability conditions for the critical points selects in each case a defined sign for the interacting parameter α.
On the other hand, the observational analysis of model Γ 1T with background data shows that the posterior distribution of the α parameter is consistent with the result from the dynamical system analysis. Besides, this model includes a change of sign for the dark energy density parameter, which must have been negative in the past. In order to distinguish interacting models from ΛCDM, we apply the statefinder analysis. The statefinder diagrams show that for some of the models, Γ 2T , Γ 1c , Γ 2d and Γ 2c , the trajectories will eventually approach to the ΛCDM fixed point in the r − s and r − q planes (see Fig.13). Nevertheless, models that asymptotically tend to a scaling state, Γ 1T , Γ 1d and Γ 1x , do not follow this lineament and have a different asymptotic behavior than ΛCDM in the r − q and r − s planes (see Fig.12). Furthermore, the statefinder diagrams indicate that the interaction has a significant effect on the evolution of the universe, and the interacting models can be distinguished from the ΛCDM model by this analysis.