e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>
11institutetext: Núcleo de Matemática, Física y Estadística, Universidad
Mayor, Temuco, Chile. 22institutetext: Departamento de Física, Universidad del
Bío-Bío, Casilla 5-C, Concepción, Chile. 33institutetext: Centro de Ciencias
Exactas, Universidad del Bío-Bío, Casilla 447, Chillán, Chile.

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

Fabiola Arévaloe1,addr1  Antonella Cide2,addr2,addr2b

(Received: date / Accepted: date)

###### Abstract

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.

##  1 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
Planck:2018vyg , but to the present understanding of physics, both components
have unknown nature Weinberg:1988cp ; Lin:2019uvt . Usually we consider the
standard cosmological scenario $\Lambda$CDM, where the dark components evolve
independently. A different approach is to consider that both components
interact coincidence ; interaction . Interacting scenarios allow to alleviate
the coincidence problem coincidence  and also are useful in addressing the
tension in the $H_{0}$ parameter H0 . 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 Wei:2010cs ; Arevalo:2019axj ; Pan:2019jqh .

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. Copeland:2006wr  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. Bahamonde:2017ize . Both works include interacting scenarios.

The authors of Ref. Sahni:2002fz  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 allows to
differentiate among dark energy models. For instance, in Ref. Zimdahl:2003wg
the authors use the statefinder parameters to characterize interacting
quintessence models of dark energy.

In Refs. Panotopoulos:2020kpo ; Huang:2021zgj  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 $\Lambda$CDM.

This work is organized as follows. In section 2 we describe interacting
scenarios and the class of sign-changeable 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.

##  2 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:

$\displaystyle\rho^{\prime}_{b}+3\rho_{b}$ $\displaystyle=$ $\displaystyle 0,$
(1) $\displaystyle\rho^{\prime}_{r}+4\rho_{r}$ $\displaystyle=$ $\displaystyle
0,$ (2) $\displaystyle\rho^{\prime}_{c}+3\rho_{c}$ $\displaystyle=$
$\displaystyle-3\Gamma,$ (3) $\displaystyle\rho^{\prime}_{x}$ $\displaystyle=$
$\displaystyle 3\Gamma,$ (4)

where by convenience we use $\eta=\ln a$ as time variable and we define
$(\,)^{\prime}=d/d\eta$ (where $a$ is the scale factor). Notice that
$\Gamma>0$ indicates an energy transfer from dark matter to dark energy and
$\Gamma<0$ indicates the opposite.

In this work we study the following sign-changeable linear interactions,

$\displaystyle\Gamma_{1i}$ $\displaystyle=$
$\displaystyle\alpha_{1i}\,q\,\rho_{i},\qquad i=T,d,c,x$ (5)
$\displaystyle\Gamma_{2j}$ $\displaystyle=$
$\displaystyle\alpha_{2j}\,q\,\rho^{\prime}_{j},\qquad j=T,d,c$ (6)

where $\alpha_{1i}$ and $\alpha_{2j}$ are constants, $q$ is the deceleration
parameter defined as $q=-1-H^{\prime}/H$ and $H$ is the Hubble expansion rate.
The subscripts $T$ and $d$ denote total energy density and dark sector energy
density, ($\rho_{c}+\rho_{x}$), respectively. It was noticed in
Ref.Arevalo:2016epc  that interaction $\Gamma_{1T}$ has analytical solutions
for the dark sector energy densities, however the other scenarios do not, to
the authors knowledge. We do not include $\Gamma_{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 $\Lambda$CDM
scenario when the interacting parameter $\alpha_{1i}$ or $\alpha_{2j}$ is set
to zero.

##  3 Dynamical System Analysis

In this section we apply dynamical system methods Copeland:2006wr ;
Bahamonde:2017ize  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 $\Omega_{i}=\rho_{i}/3H^{2}$ (dimensionless) and we
use the Friedmann equation as a constraint among parameters, i.e.,

$\Omega_{r}+\Omega_{b}+\Omega_{c}+\Omega_{x}=1,$ (7)

to reduce the system of equations (1)-(4) to:

$\displaystyle\Omega_{r}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{r}(-1+\Omega_{r}-3\Omega_{x}),$ (8)
$\displaystyle\Omega_{c}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{\Gamma}{H^{2}}+\Omega_{c}(\Omega_{r}-3\Omega_{x}),$ (9)
$\displaystyle\Omega_{x}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\Gamma}{H^{2}}+\Omega_{x}(3+\Omega_{r}-3\Omega_{x}),$ (10)

where interactions (5)–(6) are given in Table 1 in terms of the density
parameters.

$\Gamma_{1T}=\frac{3}{2}\alpha H^{2}\left(1+\Omega_{r}-3\Omega_{x}\right)$
---
$\Gamma_{1d}=\frac{3}{2}\alpha
H^{2}\left(1+\Omega_{r}-3\Omega_{x}\right)\left(\Omega_{c}+\Omega_{x}\right)$
$\Gamma_{1c}=\frac{3}{2}\alpha
H^{2}\left(1+\Omega_{r}-3\Omega_{x}\right)\Omega_{c}$
$\Gamma_{1x}=\frac{3}{2}\alpha
H^{2}\left(1+\Omega_{r}-3\Omega_{x}\right)\Omega_{x}$
$\Gamma_{2T}=-\frac{3}{2}\alpha
H^{2}\left(1+\Omega_{r}-3\Omega_{x}\right)\left(3+\Omega_{r}-3\Omega_{x}\right)$
$\Gamma_{2d}=-\frac{9}{2}\alpha H^{2}(1+\Omega_{r}-3\Omega_{x})\Omega_{c}$
$\Gamma_{2c}=-9\alpha
H^{2}\frac{\Omega_{c}(1+\Omega_{r}-3\Omega_{x})}{2+3\alpha(1+\Omega_{r}-3\Omega_{x})}$
Table 1:  Revised sign-changeable interactions

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,

$\displaystyle\Omega_{i}^{\prime}=f_{i}(\Omega_{l}),$ (11)

where $f_{i}$ is a function of the density parameters $\Omega_{r}$,
$\Omega_{c}$ and $\Omega_{x}$ and $i=r,c,x$. From this we find the critical
points $\Omega_{l}^{*}$ of the system (11) by calculating

$f_{i}(\Omega_{l}^{*})=0.$

Then we linearized the set of equations (11) around the critical points,

$\displaystyle\delta\Omega_{i}^{\prime}=J_{i}^{l}(\Omega_{j}^{*})\Omega_{l},$

where $J_{i}^{l}=\frac{\partial f_{i}}{\partial\Omega_{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 $\Gamma_{1T}$, $\Gamma_{1x}$ and $\Gamma_{2T}$ which can be
reduced to two-dimensional systems. Subsequently we describe models
$\Gamma_{1d}$, $\Gamma_{1c}$, $\Gamma_{2d}$ and $\Gamma_{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.

###  3.1 Case A: $\Gamma_{1T}$

For the model $\Gamma_{1T}$ the set of Eqs.(8)-(10) becomes,

$\displaystyle\Omega_{r}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{r}(-1+\Omega_{r}-3\Omega_{x}),$ (12)
$\displaystyle\Omega_{c}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{3\alpha}{2}\left(1+\Omega_{r}-3\Omega_{x}\right)+\Omega_{c}(\Omega_{r}-3\Omega_{x}),$
(13) $\displaystyle\Omega_{x}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{3\alpha}{2}\left(1+\Omega_{r}-3\Omega_{x}\right)+\Omega_{x}(3+\Omega_{r}-3\Omega_{x}).$
(14)

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 $\Omega_{r}^{\prime}=0$ and $\Omega_{x}^{\prime}=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 $\Gamma_{1T}$ we find the following three critical points
$A_{i}=\\{\Omega_{r},\Omega_{x}\\}$,

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\\{1-\frac{9}{4}\alpha,-\frac{3}{4}\alpha\\},$ $\displaystyle
A_{2}$ $\displaystyle=$
$\displaystyle\left\\{0,\frac{1}{4}\left(2-3\alpha-\sqrt{4-\alpha\left(4-9\alpha\right)}\right)\right\\},$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\left\\{0,\frac{1}{4}\left(2-3\alpha+\sqrt{4-\alpha\left(4-9\alpha\right)}\right)\right\\},$

where the matter contribution to the critical points
$\Omega_{m}=\Omega_{b}+\Omega_{c}$ is determined by the constraint

$\Omega_{r}+\Omega_{m}+\Omega_{x}=1.$ (15)

In the physical analysis of these points we have considered as existence
conditions $\Omega_{r}\geq 0$ and $\Omega_{m}\geq 0$ at the critical points.
On the other hand, notice that the effective state parameter $\omega_{\rm
eff}=\frac{p_{T}}{\rho_{T}}$ ($p_{T}$ is the total pressure) is related to the
deceleration parameter by $q=\frac{1}{2}(1+3\omega_{\rm eff})$, and in terms
of the density parameters $q$ is given by

$q=\frac{1}{2}(1+\Omega_{r}-3\Omega_{x}).$ (16)

Table 2 shows the existence and stability conditions on the eigenvalues
$A_{1}-A_{3}$, as well as the $\omega_{\rm eff}$ value. The critical point
$A_{1}$ in Table 2 describes a radiation dominated epoch where radiation,
matter and dark energy coexist for $0<\alpha<\frac{4}{9}$. The matter
contribution to this critical point is given by $\Omega_{m}=3\alpha$ and in
this case $\Omega_{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 $\alpha<\frac{4}{9}$.

Critical Points | $\omega_{\rm eff}$ | Existence | Stability Conditions
---|---|---|---
$A_{1}$ | $\frac{1}{3}$ | $0<\alpha<\frac{4}{9}$ |  unstable for $\alpha<\frac{4}{9}$
saddle for $\alpha>\frac{4}{9}$
$A_{2}$ | $-\frac{1}{2}+\frac{1}{2}\sqrt{1-\frac{\alpha\left(4-9\alpha\right)}{4}}+\frac{3\alpha}{4}$ | $\alpha\in\mathbb{R}$ | unstable for $\alpha>\frac{4}{9}$
saddle for $\alpha<\frac{4}{9}$
$A_{3}$ | $-\frac{1}{2}-\frac{1}{2}\sqrt{1-\frac{\alpha\left(4-9\alpha\right)}{4}}+\frac{3\alpha}{4}$ | $\alpha>0$ | stable
Table 2:  Description of the critical points for scenario $\Gamma_{1T}$.

The critical point $A_{2}$ has $\omega_{\rm eff}>-\frac{1}{3}$ irrespective of
the value of $\alpha$. In the limit $\alpha\rightarrow 0$, this phase can be
understood as a phase of matter domination, given that $\omega_{\rm
eff}\rightarrow 0$ in this limit. In this phase matter and dark energy
coexists. At this critical point the matter component is
$\Omega_{m}=\frac{1}{4}\left(2+3\alpha+\sqrt{4-\alpha(4-9\alpha)}\right)$,
which is positive irrespective of the $\alpha$ value. On the other hand, we
have $\Omega_{x}<0$ for $\alpha>0$ and $\Omega_{x}>0$ for $\alpha<0$. $A_{2}$
is a saddle critical point for $\alpha<\frac{4}{9}$ or an unstable critical
point for $\alpha>\frac{4}{9}$.

The critical point $A_{3}$ dynamically corresponds to a phase of dark energy
domination ($\omega_{\rm eff}<-\frac{1}{3}$ irrespective of the $\alpha$
value), where matter and dark energy coexist. At this point
$\Omega_{m}=\frac{1}{4}\left(2+3\alpha-\sqrt{4-\alpha(4-9\alpha)}\right)$,
which is positive for $\alpha>0$. Under this existence condition $\Omega_{x}$
is always positive. $A_{3}$ is a stable critical point for $\alpha>0$.

Figure 1: Phase plot for model $\Gamma_{1T}$ with $\alpha=0.01$.

Figure 1 shows the 2-dimensional phase space $\Omega_{x}-\Omega_{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 $\Gamma_{1T}$ for $0\leq\alpha<\frac{4}{9}$, nevertheless, in this range
$\Omega_{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
$\\{\Omega_{m},\Omega_{x}\\}=\\{0.3111,0.6889\\}$ Planck:2018vyg .

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 $\Lambda$CDM
scenario in the limit $\alpha\rightarrow 0$.

###  3.2 Case B: $\Gamma_{1x}$

For the model $\Gamma_{1x}$ we can reduce the set of Eqs.(8)-(10) to a
2-dimensional set including Eq.(8) and the following equation,

$\displaystyle\Omega_{x}^{\prime}=\Omega_{x}\left(\frac{3\alpha}{2}\left(1+\Omega_{r}-3\Omega_{x}\right)+3+\Omega_{r}-3\Omega_{x}\right).$

For this model we find three critical points
$B_{i}=\\{\Omega_{r},\Omega_{x}\\}$,

$\displaystyle B_{1}=\\{1,0\\},\quad B_{2}=\left\\{0,0\right\\},\quad
B_{3}=\left\\{0,\frac{2+\alpha}{2+3\alpha}\right\\},$

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 $\omega_{\rm eff}$ parameter.

Critical Points | $\omega_{\rm eff}$ | Existence | Stability Conditions
---|---|---|---
$B_{1}$ | $\frac{1}{3}$ | $\alpha\in\mathbb{R}$ |  unstable $\alpha>-\frac{4}{3}$
saddle $\alpha<-\frac{4}{3}$
$B_{2}$ | $0$ | $\alpha\in\mathbb{R}$ | stable for $\alpha<-2$
saddle for $\alpha>-2$
$B_{3}$ | $-\frac{2+\alpha}{2+3\alpha}$ | $\alpha>0$ or $\alpha<-\frac{2}{3}$ |  stable for $\alpha>-\frac{2}{3}$ or $-2<\alpha<-\frac{4}{3}$
saddle for $\alpha<-2$ or $-\frac{4}{3}<\alpha<-\frac{2}{3}$
Table 3:  Description of the critical points for scenario $\Gamma_{1x}$.

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 $\alpha>-\frac{4}{3}$ and saddle for
$\alpha<-\frac{4}{3}$, meanwhile, $B_{2}$ is saddle for $\alpha>-2$ and stable
for $\alpha<-2$.

At the critical point $B_{3}$ matter and dark energy coexist. At this point
$\Omega_{m}=\frac{2\alpha}{2+3\alpha}$ which is positive for
$\alpha<-\frac{2}{3}$ or $\alpha>0$. On the other hand, for $\alpha>0$, this
point corresponds to dark energy domination with $\Omega_{x}>0$ as a stable
critical point. For $\alpha<-\frac{2}{3}$ we get $\omega_{\rm
eff}<-\frac{1}{3}$, $\Omega_{m}>\Omega_{x}$ and $\Omega_{x}>0$ for
$\alpha<-2$. The corresponding stability conditions in this case are given in
Table 3.

Figure 2: Phase plot for model $\Gamma_{1x}$ with $\alpha=0.01$.

The figure 2 shows the 2-dimensional phase space $\Omega_{x}-\Omega_{r}$ for
the $\Gamma_{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 $\Gamma_{1x}$ for $\alpha\geq 0$,
where $\Omega_{x}$ is positive at the critical point $B_{3}$. In the figure
the black dot labeled as “Present” considers the values in Ref.Planck:2018vyg
.

As well as model $\Gamma_{1T}$, the critical points $B_{1}-B_{3}$ reduce to
the $\Lambda$CDM critical points in the limit $\alpha\rightarrow 0$.

###  3.3 Case C: $\Gamma_{2T}$

For the model $\Gamma_{2T}$ we can reduce the set of Eqs.(8)-(10) to a
2-dimensional set including Eq.(8) and the following equation,

$\displaystyle\Omega_{x}^{\prime}=\left(3+\Omega_{r}-3\Omega_{x}\right)\left(\Omega_{x}-\frac{3\alpha}{2}\left(1+\Omega_{r}-3\Omega_{x}\right)\right).$

For this model we find three critical points
$C_{i}=\\{\Omega_{r},\Omega_{x}\\}$,

$\displaystyle C_{1}=\\{1+9\alpha,3\alpha\\},\,\quad
C_{2}=\left\\{0,\frac{3\alpha}{2+9\alpha}\right\\},\,\quad C_{3}=\\{0,1\\},$

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 $\omega_{\rm eff}$ parameter.

At the critical point $C_{1}$ coexist radiation, matter and dark energy in a
radiation dominated era ($\omega_{\rm eff}=\frac{1}{3}$). For this critical
point we have $\Omega_{m}=-12\alpha$ which is positive for $\alpha<0$ but
necessarily in this case $\Omega_{x}<0$. At the critical point $C_{2}$ coexist
matter and dark energy, $\omega_{\rm eff}\geq-\frac{1}{3}$ for
$\alpha>-\frac{2}{9}$. At this point we have
$\Omega_{m}=\frac{2+6\alpha}{2+9\alpha}$ which is positive for
$\alpha<-\frac{1}{3}$ or $\alpha>-\frac{2}{9}$. In the range
$\alpha<-\frac{2}{9}$ or $\alpha>0$, the parameter $\Omega_{x}$ is also
positive. This point corresponds to matter domination in the limit
$\alpha\rightarrow 0$.

Figure 3: Phase plot for model $\Gamma_{2T}$ with $\alpha=-0.01$.

The point $C_{3}$ corresponds to a de-Sitter attractor for
$\alpha>-\frac{1}{3}$, where only the dark energy term contributes.

The figure 3 shows the 2-dimensional phase space $\Omega_{x}-\Omega_{r}$ for
the $\Gamma_{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 $\Gamma_{2T}$ for $-\frac{1}{9}<\alpha\leq 0$.
However, in this range, $\Omega_{x}$ is negative at the critical points
$C_{1}$ and $C_{2}$. In the figure the “Present” dot considers the values in
Ref.Planck:2018vyg .

Critical Points | $\omega_{\rm eff}$ | Existence | Stability Conditions
---|---|---|---
$C_{1}$ | $\frac{1}{3}$ | $-\frac{1}{9}<\alpha<0$ |  unstable for $\alpha>-\frac{1}{9}$
saddle for $\alpha<-\frac{1}{9}$
$C_{2}$ | $-\frac{3\alpha}{2+9\alpha}$ | $\alpha>-\frac{2}{9}$ or $\alpha<-\frac{1}{3}$ | unstable for $-\frac{2}{9}<\alpha<-\frac{1}{9}$
stable for $\alpha<-\frac{1}{3}$
saddle for $-\frac{1}{3}<\alpha<-\frac{2}{9}$ or $\alpha>-\frac{1}{9}$
$C_{3}$ | -1 | $\alpha\in\mathbb{R}$ |  stable for $\alpha>-\frac{1}{3}$
saddle for $\alpha<-\frac{1}{3}$
Table 4:  Description of the critical points for scenario $\Gamma_{2T}$.

The critical points $C_{1}-C_{3}$ corresponding to model $\Gamma_{2T}$ in
Table 4 reduce to the $\Lambda$CDM critical points in the limit
$\alpha\rightarrow 0$.

###  3.4 Case D: $\Gamma_{1d}$

For the model $\Gamma_{1d}$ the set of Eqs.(8)-(10) reduces to,

$\displaystyle\Omega_{r}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{r}(-1+\Omega_{r}-3\Omega_{x}),$
$\displaystyle\Omega_{c}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{c}(\Omega_{r}-3\Omega_{x})-\frac{3\alpha}{2}\left(1+\Omega_{r}-3\Omega_{x}\right)\left(\Omega_{c}+\Omega_{x}\right),$
$\displaystyle\Omega_{x}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{x}(3+\Omega_{r}-3\Omega_{x})+\frac{3\alpha}{2}\left(1+\Omega_{r}-3\Omega_{x}\right)\left(\Omega_{c}+\Omega_{x}\right).$

For this model we find the following critical points
$D_{i}=\\{\Omega_{r},\Omega_{b},\Omega_{x}\\}$,

$\displaystyle D_{1}=\\{1,0,0\\},\quad D_{2}=\\{0,1,0\\},$ $\displaystyle
D_{3}=\left\\{0,0,\frac{1}{4}\left(2-3\alpha-\sqrt{4-\alpha\left(4-9\alpha\right)}\right)\right\\},$
$\displaystyle
D_{4}=\left\\{0,0,\frac{1}{4}\left(2-3\alpha+\sqrt{4-\alpha\left(4-9\alpha\right)}\right)\right\\},$

where the $\Omega_{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 $\Omega_{r}\geq 0$, $\Omega_{b}\geq 0$ and $\Omega_{c}\geq 0$ at
the critical points.

Critical Points | $\omega_{\rm eff}$ | Existence | Stability Conditions
---|---|---|---
$D_{1}$ | $\frac{1}{3}$ | $\alpha\in\mathbb{R}$ |  unstable for $-\frac{1}{4}\leq\alpha<\frac{4}{9}$
saddle for $\alpha>\frac{4}{9}$
$D_{2}$ | $0$ | $\alpha\in\mathbb{R}$ | stable for $\alpha=\frac{8}{9}$
saddle for $-\frac{1}{2}\leq\alpha<\frac{8}{9}$ or $\alpha>\frac{8}{9}$
$D_{3}$ | $-\frac{1}{2}+\frac{1}{2}\sqrt{1-\frac{\alpha\left(4-9\alpha\right)}{4}}+\frac{3\alpha}{4}$ | $\alpha\in\mathbb{R}$ | unstable for $\alpha>\frac{4}{9}$
saddle for $\alpha<\frac{4}{9}$
$D_{4}$ | $-\frac{1}{2}-\frac{1}{2}\sqrt{1-\frac{\alpha\left(4-9\alpha\right)}{4}}+\frac{3\alpha}{4}$ | $\alpha>0$ | stable
Table 5:  Description of the critical points for scenario $\Gamma_{1d}$.

The critical points $D_{1}$ and $D_{2}$ in Table 5 describe a radiation and a
baryon dominated epoch, respectively, irrespective of the $\alpha$ value. The
signs of the eigenvalues of the linearized system indicate that $D_{1}$ is an
unstable critical point in the range $-\frac{1}{4}\leq\alpha<\frac{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 $\Gamma_{1T}$ and $\Gamma_{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 $\alpha\rightarrow 0$, this point (and its analysis) is equivalent to
the point $A_{2}$ for the scenario $\Gamma_{1T}$. The critical point $D_{4}$
corresponds to dark energy domination for $\alpha>0$. This point and its
analysis is equivalent to the points $A_{3}$ in scenario $\Gamma_{1T}$. The
existence conditions in Table 5 indicate that these critical points (and the
radiation one) can be find for $0\leq\alpha<\frac{4}{9}$, but in this range
necessarily $\Omega_{x}$ becomes negative at the critical point $D_{3}$. In
Fig.4 we show a projected 2-dimensional phase space $\Omega_{x}-\Omega_{c}$
where we have fixed $\Omega_{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}$.

Figure 4: Projected phase plot for model $\Gamma_{1d}$ with $\alpha=0.01$

###  3.5 Case E: $\Gamma_{1c}$

For the model $\Gamma_{1c}$ the set of Eqs.(8)-(10) reduces to,

$\displaystyle\Omega_{r}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{r}(-1+\Omega_{r}-3\Omega_{x}),$
$\displaystyle\Omega_{c}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{c}(\Omega_{r}-3\Omega_{x})-\frac{3}{2}\alpha\left(1+\Omega_{r}-3\Omega_{x}\right)\Omega_{c},$
$\displaystyle\Omega_{x}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{x}(3+\Omega_{r}-3\Omega_{x})+\frac{3}{2}\alpha\left(1+\Omega_{r}-3\Omega_{x}\right)\Omega_{c}.$

For this model we find the following critical points
$E_{i}=\\{\Omega_{r},\Omega_{b},\Omega_{x}\\}$,

$\displaystyle E_{1}=\\{1,0,0\\},\quad E_{2}=\\{0,1,0\\},\quad
E_{3}=\left\\{0,0,-\frac{\alpha}{2-3\alpha}\right\\},\quad
E_{4}=\left\\{0,0,1\right\\},$

where the contribution $\Omega_{c}$ to the critical points is determined by
the constraint (7) to be zero for each point but $E_{3}$, where
$\Omega_{c}=\frac{2(1-\alpha)}{2-3\alpha}$.

Critical Points | $\omega_{\rm eff}$ | Existence | Stability Conditions
---|---|---|---
$E_{1}$ | $\frac{1}{3}$ | $\alpha\in\mathbb{R}$ |  unstable for $\alpha<\frac{1}{3}$
saddle for $\alpha>\frac{1}{3}$
$E_{2}$ | $0$ | $\alpha\in\mathbb{R}$ | stable for $\alpha=\frac{2}{3}$
saddle for $\alpha\neq\frac{2}{3}$
$E_{3}$ | $\frac{\alpha}{2-3\alpha}$ | $\alpha<\frac{2}{3}$ or $\alpha>1$ | unstable for $\frac{1}{3}<\alpha<\frac{2}{3}$
stable for $\alpha>1$
saddle for $\alpha<\frac{1}{3}$ or $\frac{2}{3}<\alpha<1$
$E_{4}$ | -1 | $\alpha\in\mathbb{R}$ |  stable for $\alpha<1$
saddle for $\alpha>1$

Table 6: Description of the critical points for scenario $\Gamma_{1c}$. Figure
5: Projected phase plot for model $\Gamma_{1c}$ with $\alpha=-0.01$

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
$\alpha<\frac{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
$\Omega_{c}>0$ for $\alpha<\frac{2}{3}$ or $\alpha>1$ and $\omega_{\rm
eff}>-\frac{1}{3}$ for $\alpha<\frac{2}{3}$, this point corresponds to dark
matter domination in the limit $\alpha\rightarrow 0$. The point $E_{4}$
corresponds to a de-Sitter stage with the dominance of dark energy. This point
is stable for $\alpha<1$.

From Table 6 we notice that the critical points corresponding to unstable,
saddle and stable behavior can be find for $\alpha<\frac{1}{3}$. Considering
this range, $\Omega_{x}$ remains positive for $\alpha<0$ at the critical point
$E_{3}$. In Fig.5 we show the projected 2-dimensional phase space
$\Omega_{x}-\Omega_{c}$ where we have fixed $\Omega_{r}=0$ and we have chosen
$\alpha=-0.01$. In this figure we show two critical points (physical)
corresponding to the saddle point $E_{3}$ and the attractor point $E_{4}$.

###  3.6 Case F: $\Gamma_{2d}$

For the model $\Gamma_{2d}$ the set of Eqs.(8)-(10) reduces to,

$\displaystyle\Omega_{r}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{r}(-1+\Omega_{r}-3\Omega_{x}),$
$\displaystyle\Omega_{c}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{9}{2}\alpha(1+\Omega_{r}-3\Omega_{x})\Omega_{c}+\Omega_{c}(\Omega_{r}-3\Omega_{x}),$
$\displaystyle\Omega_{x}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{9}{2}\alpha(1+\Omega_{r}-3\Omega_{x})\Omega_{c}+\Omega_{x}(3+\Omega_{r}-3\Omega_{x}).$

For this model we find the following critical points
$F_{i}=\\{\Omega_{r},\Omega_{b},\Omega_{x}\\}$,

$\displaystyle F_{1}=\\{1,0,0\\},\quad F_{2}=\\{0,1,0\\},\quad
F_{3}=\left\\{0,0,\frac{3\alpha}{2+9\alpha}\right\\},\quad
F_{4}=\left\\{0,0,1\right\\},$

where the contribution $\Omega_{c}$ to the critical points is determined by
the constraint (7) to be zero at each point but $F_{3}$, where
$\Omega_{c}=\frac{2+6\alpha}{2+9\alpha}$.

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
$\alpha>-\frac{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
$\Omega_{c}>0$ for $\alpha<-\frac{1}{3}$ or $\alpha>-\frac{2}{9}$ and
$\omega_{\rm eff}>-\frac{1}{3}$ for $\alpha>-\frac{2}{9}$, this point
corresponds to dark matter domination in the limit $\alpha\rightarrow 0$. The
point $F_{3}$ is saddle in the range $\alpha>-\frac{1}{9}$. The point $F_{4}$
corresponds to a de-Sitter stage with the dominance of dark energy and this
point is stable for $\alpha>-\frac{1}{3}$.

From Table 7 we notice that the critical points corresponding to unstable,
saddle and stable can be find for $\alpha>-\frac{1}{9}$. Notice that in this
range, $\Omega_{x}$ remains positive for $\alpha>0$ at the critical point
$F_{3}$. In Fig.6 we show the projected 2-dimensional phase space
$\Omega_{x}-\Omega_{c}$ where we have fixed $\Omega_{r}=0$ and we have chosen
$\alpha=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 $\Gamma_{2T}$ and $\Gamma_{2d}$ share the same critical
points and existence/stability conditions in the points corresponding to the
dark sector, $C_{2}-C_{3}$ and $F_{3}-F_{4}$.

Critical Points | $\omega_{\rm eff}$ | Existence | Stability Conditions
---|---|---|---
$F_{1}$ | $\frac{1}{3}$ | $\alpha\in\mathbb{R}$ |  unstable for $\alpha>-\frac{1}{9}$
saddle for $\alpha<-\frac{1}{9}$
$F_{2}$ | $0$ | $\alpha\in\mathbb{R}$ | stable for $\alpha=-\frac{2}{9}$
saddle for $2+9\alpha\neq 0$
$F_{3}$ | $-\frac{3\alpha}{2+9\alpha}$ | $\alpha>-\frac{2}{9}$ or $\alpha<-\frac{1}{3}$ | unstable for $-\frac{2}{9}<\alpha<-\frac{1}{9}$
stable for $\alpha<-\frac{1}{3}$
saddle for $-\frac{1}{3}<\alpha<-\frac{2}{9}$ or $\alpha>-\frac{1}{9}$
$F_{4}$ | -1 | $\alpha\in\mathbb{R}$ |  stable for $\alpha>-\frac{1}{3}$
saddle for $\alpha<-\frac{1}{3}$
Table 7: Description of the critical points for scenario $\Gamma_{2d}$. Figure
6: Projected phase plot for model $\Gamma_{2d}$ with $\alpha=0.01$

###  3.7 Case G: $\Gamma_{2c}$

For the model $\Gamma_{2c}$ the set of Eqs.(8)-(10) reduces to,

$\displaystyle\Omega_{r}^{\prime}$ $\displaystyle=$
$\displaystyle\Omega_{r}(-1+\Omega_{r}-3\Omega_{x}),$
$\displaystyle\Omega_{c}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{9\alpha\Omega_{c}(1+\Omega_{r}-3\Omega_{x})}{2+3\alpha(1+\Omega_{r}-3\Omega_{x})}+\Omega_{c}(\Omega_{r}-3\Omega_{x}),$
$\displaystyle\Omega_{x}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{9\alpha\Omega_{c}(1+\Omega_{r}-3\Omega_{x})}{2+3\alpha(1+\Omega_{r}-3\Omega_{x})}+\Omega_{x}(3+\Omega_{r}-3\Omega_{x}).$

For this model we find the following critical points
$G_{i}=\\{\Omega_{r},\Omega_{b},\Omega_{x}\\}$,

$\displaystyle G_{1}=\\{1,0,0\\},\quad G_{2}=\\{0,1,0\\},\quad
G_{5}=\left\\{0,0,1\right\\},$ $\displaystyle
G_{3}=\left\\{0,0,\frac{1-\sqrt{1+3\alpha(4+3\alpha)}+6\alpha}{9\alpha}\right\\},$
$\displaystyle
G_{4}=\left\\{0,0,\frac{1+\sqrt{1+3\alpha(4+3\alpha)}+6\alpha}{9\alpha}\right\\},$

where given Eq. (7), $\Omega_{c}$ is zero at $G_{1}$ and $G_{2}$,
$\Omega_{c}=\frac{-1+\sqrt{1+3\alpha(4+3\alpha)}+3\alpha}{9\alpha}$ at
$G_{3}$, and
$\Omega_{c}=-\frac{1+\sqrt{1+3\alpha(4+3\alpha)}-3\alpha}{9\alpha}$ 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
$\alpha<-\frac{1}{3}$ or $\alpha>-\frac{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
$\alpha<-\frac{\sqrt{3}+2}{3}$, $\frac{\sqrt{3}-2}{3}<\alpha<0$ or $\alpha>0$.
Inside these limits $G_{3}$ can be unstable for
$\frac{\sqrt{3}-2}{3}<\alpha<-\frac{1}{12}$ or saddle for
$\alpha<-\frac{\sqrt{3}+2}{3}$ or $\alpha>-\frac{1}{12}$. We notice that
$\omega_{\rm eff}>-\frac{1}{3}$ for $\alpha>\frac{\sqrt{3}-2}{3}$ and this
point reduces to dark matter domination in the limit $\alpha\rightarrow 0$.
Besides, $\Omega_{x}>0$ for $\alpha>0$ or $\alpha\leq-\frac{\sqrt{3}+2}{3}$.
At the point $G_{4}$ coexist the dark components in the range
$\frac{\sqrt{3}-2}{3}<\alpha<0$, where this point is saddle and
$\alpha\leq-\frac{\sqrt{3}+2}{3}$ where it is a stable critical point. In the
range where $G_{4}$ is a saddle point we have $\omega_{\rm eff}>-\frac{1}{3}$,
outside this range, in the allowed region, the effective fluid accelerate the
universe’s expansion. The $\Omega_{x}$ contribution of this point is positive
for $\alpha<-\frac{\sqrt{3}+2}{3}$. Notice that the point $G_{4}$ diverges in
the limit $\alpha\rightarrow 0$. The point $G_{5}$ corresponds to a de-Sitter
stage with the dominance of dark energy. This point is stable for
$\alpha<\frac{1}{3}$ and unstable outside this range.

We have sketched the behavior of the critical points $G_{3}-G_{5}$ in Fig. 8,
where we notice that for $\alpha<-\frac{\sqrt{3}+2}{3}$ we have one saddle
point and two stable points corresponding to $G_{3}$ and $G_{4}$, $G_{5}$,
respectively. In the range
$-\frac{\sqrt{3}+2}{3}<\alpha<\frac{\sqrt{3}-2}{3}$, we have no saddle point
($G_{3}$ and $G_{4}$ do not exist) and one stable point, $G_{5}$. For
$\frac{\sqrt{3}-2}{3}<\alpha<-\frac{1}{12}$ there is one unstable, one saddle
and one stable point, $G_{3}$, $G_{4}$ and $G_{5}$, respectively. In the range
$-\frac{1}{12}<\alpha<0$ we have two saddle points, $G_{3}$ and $G_{4}$, and
one stable point, $G_{5}$. For $0<\alpha<\frac{1}{3}$ we have one saddle point
and one stable point, $G_{3}$ and $G_{5}$, respectively. For
$\alpha>\frac{1}{3}$ we have two saddle points ($G_{3}$ and $G_{5}$) and no
stable point.

In Fig. 7 we show the projected 2-dimensional phase space
$\Omega_{x}-\Omega_{c}$ where we have fixed $\Omega_{r}=0$ and we have chosen
$\alpha=0.01$. In this figure we show two critical points (physical)
corresponding to the saddle point $G_{3}$ and the attractor point $G_{5}$. The
critical point $G_{4}$ does not exist for $\alpha>0$.

Figure 7: Projected phase plots for model $\Gamma_{2c}$ with $\alpha=0.01$ Critical Points | $\omega_{\rm eff}$ | Existence | Stability Conditions
---|---|---|---
$G_{1}$ | $\frac{1}{3}$ | $\alpha\in\mathbb{R}$ |  unstable for $\alpha<-\frac{1}{3}$ or $\alpha>-\frac{1}{12}$
saddle for $-\frac{1}{3}<\alpha<-\frac{1}{12}$
$G_{2}$ | $0$ | $\alpha\in\mathbb{R}$ | stable $\alpha=-\frac{1}{6}$
saddle for $\alpha<-\frac{2}{3}$ or $-\frac{2}{3}<\alpha<-\frac{1}{6}$ or
$\alpha>-\frac{1}{6}$
$G_{3}$ | $-\frac{1+6\alpha-\sqrt{1+12\alpha+9\alpha^{2}}}{9\alpha}$ |  $\alpha>0$ or $\frac{\sqrt{3}-2}{3}\leq\alpha<0$ | unstable for $\frac{\sqrt{3}-2}{3}<\alpha<-\frac{1}{12}$
or $\alpha\leq-\frac{\sqrt{3}+2}{3}$ |  saddle for $\alpha<-\frac{\sqrt{3}+2}{3}$ or $\alpha>-\frac{1}{12}$
$G_{4}$ | $-\frac{1+6\alpha+\sqrt{1+12\alpha+9\alpha^{2}}}{9\alpha}$ | $\frac{\sqrt{3}-2}{3}\leq\alpha<0$ |  saddle for $\frac{\sqrt{3}-2}{3}<\alpha<0$
$\alpha\leq-\frac{\sqrt{3}+2}{3}$ |  stable for $\alpha<-\frac{\sqrt{3}+2}{3}$ or $\alpha>0$
$G_{5}$ | -1 | $\alpha\in\mathbb{R}$ |  stable for $\alpha<\frac{1}{3}$
saddle for $\alpha>\frac{1}{3}$
Table 8: Description of the critical points for scenario $\Gamma_{2c}$.

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: $\alpha>0$ for $\Gamma_{1x}$,
$\alpha<0$ for $\Gamma_{1c}$, $\alpha>0$ for $\Gamma_{2d}$ and
$0<\alpha<\frac{1}{3}$ for $\Gamma_{2c}$. For these models $\Omega_{x}>0$ and
the sign of the interacting parameter $\alpha$ becomes defined in each case.

On the other hand, the corresponding intervals where existence/stability
conditions are jointly met for models where $\Omega_{x}<0$ are:
$0<\alpha<\frac{4}{9}$ for $\Gamma_{1T}$, $-\frac{1}{9}<\alpha<0$ for
$\Gamma_{2T}$ and $0<\alpha<\frac{4}{9}$ for $\Gamma_{1d}$. Models
$\Gamma_{1T}$ and $\Gamma_{2T}$ have $\Omega_{x}<0$ at radiation and dark
matter domination, nevertheless $\Omega_{x}$ changes sign and it becomes
positive at dark energy domination. For model $\Gamma_{1d}$ we have
$\Omega_{x}<0$ for dark matter domination, but it is positive for radiation
and dark energy domination, meaning at least two change of sign of
$\Omega_{x}$ during evolution.

Finally, notice that the stability conditions of center type critical points
have not been addressed in this work.

Figure 8: Scheme of the behavior of critical points $G_{3}-G_{5}$ in terms of
the interacting parameter $\alpha$, for model $\Gamma_{2c}$. The pale blue and
gray regions represent a saddle and stable critical point, respectively. The
boxed pale yellow regions represent an interval where the existence conditions
are not met.

##  4 Observational analysis

In the observational analysis we use background data such as the local
measurement of Hubble parameter ($H_{0}$) Riess:2019cxk , cosmic chronometers
Moresco:2016mzx , supernovae type Ia (SNe Ia) Scolnic:2017caz , baryon
acoustic oscillations (BAO) Eisenstein:2005su , and the angular scale of the
sound horizon at the last scattering Ade:2015rim . In this section we briefly
describe the used data, a detailed description can be find in Ref. Cid:2018ugy
.

###  4.1 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\pm 0.0142$ Riess:2019cxk .

###  4.2 Cosmic chronometers

We use 24 cosmic chronometers obtained through the differential age method
Jimenez:2001gg . This procedure provides cosmological-independent direct
measurements of the expansion history of the universe up to redshift 1.2
Verde:2014qea , see Table 3 in Ref. Cid:2018ugy .

###  4.3 Supernovae Type Ia

We use the Pantheon Sample, a set of 1048 spectroscopically confirmed SNe Ia
ranging from redshift $0.01$ to $2.3$ Scolnic:2017caz , 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 $\mu=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,

$\mu(z)=5\log d_{L}(z)+25,$ (17)

where the luminosity distance $d_{L}$ is given by,

$d_{L}=(1+z)\int^{z}_{0}\frac{H_{0}\,dz^{\prime}}{H(z^{\prime})}\,{\rm[Mpc]}.$
(18)

###  4.4 BAO data

We use isotropic measurements of the BAO signal from 6dFGS Beutler:2011hx ,
MGS Ross:2014qpa  and eBOSS Ata:2017dya . These measurements are given in
terms of the dimensionless ratio,

$d_{z}(z)=D_{V}(z)/r_{s}(z_{d}),$ (19)

where $z_{d}$ is the redshift at the drag epoch,

$\displaystyle D_{V}(z)$ $\displaystyle=$
$\displaystyle\left((1+z)^{2}D_{A}(z)^{2}\frac{c\,z}{H(z)}\right)^{1/3},$ (20)
$\displaystyle D_{A}(z)$ $\displaystyle=$
$\displaystyle\frac{c}{(1+z)}\int^{z}_{0}\frac{dz}{H(z)},$ (21) $\displaystyle
r_{s}(z)$ $\displaystyle=$
$\displaystyle\int^{\infty}_{z}\frac{c}{\sqrt{3(1+R)}}\frac{dz}{H(z)},$ (22)

with $c$ the speed of light and $R=\frac{3\Omega_{b}}{4\Omega_{\gamma}(1+z)}$
Eisenstein:1997ik .

Furthermore, we use the anisotropic BAO measurements from BOSS DR12
Alam:2016hwk  and Ly$\alpha$ forest duMasdesBourboux:2017mrl , which are
defined in terms of $D_{A}$ and $D_{H}=c/H(z)$, as shown in table 2 of Ref.
Evslin:2017qdn . We use these data along with the corresponding covariance
matrix in Ref. Evslin:2017qdn .

###  4.5 Cosmic microwave background data

We consider the angular scale of the sound horizon at the last scattering as
the only contribution of CMB data,

$\ell_{a}=\frac{\pi(1+z_{*})D_{A}(z_{*})}{r_{s}(z_{*})},$ (23)

where the comoving size of the sound horizon is evaluated at $z_{*}=1089.80$
Aghanim:2018eyx . We compare the value obtained in our study with the one
reported by the Planck collaboration in 2015, $\ell_{a}=301.63\pm 0.15$
Ade:2015rim .

To compute the maximum likelihood and posterior distributions we use the
MultiNest algorithm111https://github.com/JohannesBuchner/MultiNest
Feroz:2007kg ; Feroz:2008xx , 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 $\Gamma_{1T}$. The posterior distribution for the
parameters $h$, $\Omega_{c}$ and $\alpha$ are shown in Fig. 9 and Table 9.

Figure 9: Contour plots for the model $\Gamma_{1T}$, $1\sigma$ and $2\sigma$ confidence levels are shown.  $h$ | $\Omega_{c}$ | $\alpha$ | $\chi^{2}_{\rm min}$
---|---|---|---
$0.7204^{+0.0098}_{-0.011}$ | $0.2533\pm 0.0070$ | $0.0090\pm 0.0027$ | $1044.278$
Table 9: Best fit parameters and $1\sigma$ error for the observational
analysis of the $\Gamma_{1T}$ scenario.

In this analysis we fix the physical baryon density to
$\Omega_{b}h^{2}=0.022383$ Aghanim:2018eyx , the physical photon density to
$\Omega_{\gamma}h^{2}=2.469\times 10^{-5}$ Komatsu:2010fb , we consider the
the radiation density as
$\Omega_{r}=\Omega_{\gamma}\left(1+\frac{7}{8}\left(\frac{4}{11}\right)^{4/3}N_{\rm
eff}\right)$, and the effective number of neutrinos $N_{\rm eff}=3.046$
Lesgourgues:2018ncw . Notice that we have considered a flat prior for the
$\alpha$ parameter consistent with the results of section 3.1, i.e.,
$0<\alpha<\frac{4}{9}$.

Figure 10: Evolution of the density parameters for the model $\Gamma_{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 $\Gamma_{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\rightarrow\infty$, the dark energy density tends to
$\Omega_{x}\rightarrow-0.00675$ instead of 0, as in the $\Lambda$CDM model.
Finally, notice that a negative dark energy density in the past is a intrinsic
feature of model $\Gamma_{1T}$, which arise from imposing the existence
conditions in the dynamical system analysis.

##  5 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:

$\displaystyle r\equiv\frac{\dddot{a}}{aH^{3}},\quad
s\equiv\frac{r-1}{3\left(q-\frac{1}{2}\right)}.$ (24)

In terms of the density parameters $r$ and $s$ are given by:

$\displaystyle r$ $\displaystyle=$
$\displaystyle\left(2+\Omega_{r}-3\Omega_{x}\right)q-q^{\prime},$ (25)
$\displaystyle s$ $\displaystyle=$
$\displaystyle\frac{2\left(-1+q(2+\Omega_{r}-3\Omega_{x})-q^{\prime}\right)}{3(\Omega_{r}-3\Omega_{x})},$
(26)

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 $\\{\Omega_{m},\Omega_{c},\Omega_{x}\\}=\\{0.3111,0.26212,0.6889\\}$
Planck:2018vyg . 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.

Figure 11: The figures show the space parameter $r-q$ (left) and $r-s$ (right)
for the $\Gamma_{1T}$ model and several values of the $\alpha$ parameter.
Black dots indicate current values and gray arrows the evolution’s direction.
The solid black line represents $\Lambda$CDM, meanwhile, blue solid, dashed
and dotted lines indicate $\alpha=0.01,0.03,0.08$, respectively.

Figure 12: The figure shows the planes $r-q$ (left) and $r-s$ (right) for the
$\Gamma_{1T}$, $\Gamma_{1d}$ and $\Gamma_{1x}$ models in the case
$\alpha=0.01$. Black dots indicate current values and gray arrows the
evolution’s direction. The solid black line represents $\Lambda$CDM,
meanwhile, blue solid, blue dot-dashed and brown lines indicate $\Gamma_{1T}$,
$\Gamma_{1d}$ and $\Gamma_{1x}$, respectively.

Figure 13: The figure shows the planes $r-q$ (left) and $r-s$ (right) for the
models $\Gamma_{2T}$ and $\Gamma_{1c}$ in the case $\alpha=-0.01$ and models
$\Gamma_{2c}$ and $\Gamma_{2d}$ in the case $\alpha=0.01$. Black dots indicate
current values and gray arrows the evolution’s direction. The solid black line
represents $\Lambda$CDM, meanwhile, purple, red, orange and black dotted lines
correspond to models $\Gamma_{2T}$, $\Gamma_{1c}$, $\Gamma_{2d}$ and
$\Gamma_{2c}$, respectively.

The $r-q$ and $r-s$ planes for the $\Gamma_{1T}$ model are shown in Fig. 11,
here the trajectories correspond to different values of $\alpha$ consistent
with the results of section 3.1. We notice that in both planes the variation
in $\alpha$ modify the the trajectories with respect to the $\Lambda$CDM
model.

In Fig. 12 we show the $r-q$ and $r-s$ planes for the models $\Gamma_{1T}$,
$\Gamma_{1d}$ and $\Gamma_{1x}$ in the case $\alpha=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 $\Gamma_{1T}$ and $\Gamma_{1d}$ are very similar in both planes.

Fig. 13 shows the $r-q$ and $r-s$ planes for models $\Gamma_{1c}$,
$\Gamma_{2T}$, $\Gamma_{2d}$ and $\Gamma_{2c}$, we consider $\alpha=-0.01$ for
models $\Gamma_{2T}$ and $\Gamma_{2d}$ and $\alpha=0.01$ for models
$\Gamma_{2d}$ and $\Gamma_{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 $\Lambda$CDM in the
plane $r-s$, the interacting models present a loop in their trajectories,
passing through the point $(r,s)=(1,0)$ twice.

We notice that the ending point is different from $\Lambda$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. Panotopoulos:2020kpo , where two different interacting models where
studied, which do not have a change of sign in the interacting term.

##  6 Final Remarks

We have performed a dynamical system analysis and a statefinder analysis for a
class of sign-changeable interactions in the dark sector Arevalo:2019axj . 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 $\alpha$
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 $\alpha$ (see section 3). We find that models, $\Gamma_{1T}$,
$\Gamma_{2T}$, $\Gamma_{1d}$, need to have $\Omega_{x}<0$ at some critical
points in order to follow the universe’s evolution, this necessarily implies
at least one change of sign of $\Omega_{x}$ during evolution, resulting in
$\Omega_{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 $\alpha$.

On the other hand, the observational analysis of model $\Gamma_{1T}$ with
background data shows that the posterior distribution of the $\alpha$
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 $\Lambda$CDM, we apply the
statefinder analysis. The statefinder diagrams show that for some of the
models, $\Gamma_{2T}$, $\Gamma_{1c}$, $\Gamma_{2d}$ and $\Gamma_{2c}$, the
trajectories will eventually approach to the $\Lambda$CDM fixed point in the
$r-s$ and $r-q$ planes (see Fig.13). Nevertheless, models that asymptotically
tend to a scaling state, $\Gamma_{1T}$, $\Gamma_{1d}$ and $\Gamma_{1x}$, do
not follow this lineament and have a different asymptotic behavior than
$\Lambda$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 $\Lambda$CDM model by this analysis.

Acknowledgements. AC acknowledges the support of Vicerrectoría de
Investigación y Postgrado of Universidad del Bío-Bío through grant no. 2120247
IF/R.

## References

  * (1) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020)
  * (2) S. Weinberg, Rev. Mod. Phys. 61, 1-23 (1989)
  * (3) T. Lin, “Dark matter models and direct detection,” PoS 333, 009 (2019)
  * (4) L.P. Chimento, A.S. Jakubi, D. Pavon, W.Zimdahl, Phys. Rev. D 67 083513, (2003); D. Pavon and W. Zimdahl, Phys. Lett. B 628, 206 (2005); Z.-K. Guo and Y.-Z. Zhang, Phys. Rev. D 71, 023501 (2005)
  * (5) W. Zimdahl, Int. J. Mod. Phys. D 14, 2319 (2005); J. D. Barrow and T. Clifton, Phys. Rev. D 73, 103520 (2006); W. Zimdahl, Int. J. Mod. Phys. D 17, 651 (2008); F. Arevalo, A. P. R. Bacalhau and W. Zimdahl, Class. Quant. Grav. 29, 235001 (2012); R. F. vom Marttens, W. S. Hipolito-Ricaldi and W. Zimdahl, JCAP 1408, 004 (2014); F. Arevalo, P. Cifuentes, S. Lepe and F. Pena, Astrophys. Space Sci. 352, 899 (2014); B. Wang, E. Abdalla, F. Atrio-Barandela and D. Pavon, Rept. Prog. Phys. 79, no.9, 096901 (2016) R. F. vom Marttens, L. Casarini, W. Zimdahl, W. S. Hipolito-Ricaldi and D. F. Mota, Phys. Dark Univ. 15, 114 (2017); R. F. vom Marttens, L. Casarini, W. S. Hipolito-Ricaldi and W. Zimdahl, JCAP 1701, no. 01, 050 (2017). A. Cid, B. Santos, C. Pigozzo, T. Ferreira and J. Alcaniz, JCAP 03, 030 (2019); A. Cid, C. Rodriguez-Benites, M. Cataldo and G. Casanova, Eur. Phys. J. C 81, no.1, 31 (2021)
  * (6) S. Kumar and R. C. Nunes, Phys. Rev. D 96, no.10, 103511 (2017); E. Di Valentino, A. Melchiorri and O. Mena, Phys. Rev. D 96, no.4, 043503 (2017); W. Yang, S. Pan, E. Di Valentino, R. C. Nunes, S. Vagnozzi and D. F. Mota, JCAP 09, 019 (2018); W. Yang, A. Mukherjee, E. Di Valentino and S. Pan, Phys. Rev. D 98, no.12, 123527 (2018); S. Kumar, R. C. Nunes and S. K. Yadav, Eur. Phys. J. C 79, no.7, 576 (2019)
  * (7) H. Wei, Commun. Theor. Phys. 56, 972-980 (2011)
  * (8) F. Arevalo, A. Cid, L. P. Chimento and P. Mella, Eur. Phys. J. C 79, no.4, 355 (2019)
  * (9) S. Pan, W. Yang, C. Singha and E. N. Saridakis, Phys. Rev. D 100, no.8, 083539 (2019)
  * (10) E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753-1936 (2006)
  * (11) S. Bahamonde, C. G. Böhmer, S. Carloni, E. J. Copeland, W. Fang and N. Tamanini, Phys. Rept. 775-777, 1-122 (2018)
  * (12) V. Sahni, T. D. Saini, A. A. Starobinsky and U. Alam, JETP Lett. 77, 201-206 (2003)
  * (13) W. Zimdahl and D. Pavon, Gen. Rel. Grav. 36, 1483-1491 (2004)
  * (14) G. Panotopoulos, Á. Rincón, G. Otalora and N. Videla, Eur. Phys. J. C 80, no.3, 286 (2020)
  * (15) Q. Huang, H. Huang, B. Xu, F. Tu and J. Chen, Eur. Phys. J. C 81 no.8, 686 (2021)
  * (16) F. Arevalo, A. Cid and J. Moya, Eur. Phys. J. C 77, no.8, 565 (2017)
  * (17) A. Cid, B. Santos, C. Pigozzo, T. Ferreira and J. Alcaniz, JCAP 03, 030 (2019)
  * (18) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri and D. Scolnic, Astrophys. J. 876, no.1, 85 (2019)
  * (19) M. Moresco et al., JCAP 1605, no. 05, 014 (2016)
  * (20) D. M. Scolnic et al., Astrophys. J. 859, no. 2, 101 (2018)
  * (21) D. J. Eisenstein et al. [SDSS Collaboration], Astrophys. J. 633, 560 (2005)
  * (22) P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. 594, A14 (2016)
  * (23) R. Jimenez and A. Loeb, Astrophys. J. 573, 37 (2002)
  * (24) L. Verde, P. Protopapas and R. Jimenez, Phys. Dark Univ. 5-6, 307 (2014)
  * (25) F. Beutler et al., Mon. Not. Roy. Astron. Soc. 416, 3017 (2011)
  * (26) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden and M. Manera, Mon. Not. Roy. Astron. Soc. 449, no. 1, 835 (2015)
  * (27) M. Ata et al., Mon. Not. Roy. Astron. Soc. 473, no. 4, 4773 (2018)
  * (28) D. J. Eisenstein and W. Hu, Astrophys. J. 496, 605 (1998)
  * (29) S. Alam et al. [BOSS Collaboration], Mon. Not. Roy. Astron. Soc. 470, no. 3, 2617 (2017)
  * (30) H. du Mas des Bourboux, J. M. Le Goff, M. Blomqvist, N. G. Busca, J. Guy, J. Rich, C. Yèche, J. E. Bautista, E. Burtin and K. S. Dawson, et al. Astron. Astrophys. 608, A130 (2017)
  * (31) J. Evslin, A. A. Sen and Ruchika, Phys. Rev. D 97, no. 10, 103511 (2018)
  * (32) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020)
  * (33) F. Feroz and M. Hobson, Mon. Not. Roy. Astron. Soc. 384, 449 (2008)
  * (34) F. Feroz, M. Hobson and M. Bridges, Mon. Not. Roy. Astron. Soc. 398, 1601-1614 (2009)
  * (35) J. Lesgourgues, G. Mangano, G. Miele and S. Pastor, “Neutrino Cosmology,” Cambridge University Press (2013)
  * (36) E. Komatsu et al. [WMAP Collaboration], Astrophys. J. Suppl. 192, 18 (2011)