# Phase space analysis of sign-shifting interacting dark energy models

Sudip Halder<EMAIL_ADDRESS>Department of Mathematics, Presidency
University, 86/1 College Street, Kolkata 700073, India  Jaume de Haro
<EMAIL_ADDRESS>Departament de Matemàtiques, Universitat Politècnica de
Catalunya, Diagonal 647, 08028 Barcelona, Spain  Tapan Saha
<EMAIL_ADDRESS>Department of Mathematics, Presidency University,
86/1 College Street, Kolkata 700073, India  Supriya Pan
<EMAIL_ADDRESS>Department of Mathematics, Presidency
University, 86/1 College Street, Kolkata 700073, India  Institute of Systems
Science, Durban University of Technology, PO Box 1334, Durban 4000, Republic
of South Africa

###### Abstract

The theory of non-gravitational interaction between a pressure-less dark
matter (DM) and dark energy (DE) is a phenomenologically rich cosmological
domain which has received magnificent attention in the community. In the
present article we have considered some interacting scenarios with some novel
features: the interaction functions do not depend on the external parameters
of the universe, rather, they depend on the intrinsic nature of the dark
components; the assumption of unidirectional flow of energy between DM and DE
has been extended by allowing the possibility of bidirectional energy flow
characterized by some sign shifting interaction functions; and the DE equation
of state has been considered to be either constant or dynamical in nature.
These altogether add new ingredients in this context, and, we performed the
phase space analysis of each interacting scenario in order to understand their
global behaviour. According to the existing records in the literature, this
combined picture has not been reported elsewhere. From the analyses, we
observed that the DE equation of state as well as the coupling parameter(s) of
the interaction models can significantly affect the nature of the critical
points. It has been found that within these proposed sign shifting interacting
scenarios, it is possible to obtain stable late time attractors which may act
as global attractors corresponding to an accelerating expansion of the
universe. The overall outcomes of this study clearly highlight that the sign
shifting interaction functions are quite appealing in the context of
cosmological dynamics and they deserve further attention.

###### pacs:

98.80.-k, 95.36.+x, 95.35.+d, 98.80.Es

##  I Introduction

Over the last two decades, dynamics of the universe has been surprisingly
thrilling due to the availability of a large amount of observational data. At
the end of nineties, observations from Type Ia supernovae first reported that
our universe is passing through a phase of accelerated expansion [1, 2] and
this accelerated expansion is supposed to be driven by the presence of some
exotic matter sector in our universe sector having large negative pressure.
This exotic matter can be described in various ways. Two well known approaches
are the modification of the matter sector of the universe in the context of
Einstein’s General Relativity, dubbed as dark energy (DE) [3, 4, 5], or, the
modification of the gravitational sector of the universe in various ways [6,
7, 8, 9, 10, 11, 12, 13, 14, 15], known as geometrical DE. Apart from DE or
geometrical DE, our universe sector also contains a non-luminous dark matter
(DM) fluid responsible for the observed structure formation of the universe.
And according to the high precision data from several astronomical missions,
nearly, 68% of the total energy density of the universe is occupied by either
DE or geometrical DE and more or less 28% of the total energy density of the
universe is occupied by DM, that means, nearly 96% of the total energy budget
of the universe is comprised by DE and DM. Thus, the dynamics of the universe
is heavily dependent on the dark sector (DE$+$DM) of the universe. However,
despite many astronomical missions, the nature, origin and the evolution of
the dark sector have remained mysterious so far and probing the physics of the
dark sector has been one of the challenges for modern cosmology at the present
moment. In order to describe the present universe, several cosmological models
have been proposed and investigated by several investigators. Among these
models, the $\Lambda$-Cold Dark Matter ($\Lambda$CDM) cosmological model
constructed in the framework of GR in which the cosmological constant
$\Lambda$ plugged into the Einstein’s gravitational equations acts as the
source of DE and DM is cold (pressure-less), has been found to be extremely
successful in the light of a number of observational datasets. Nevertheless,
$\Lambda$CDM faces several theoretical and observational challenges, and
therefore, a revision of the standard $\Lambda$CDM cosmology has been
suggested in recent times.

A popular revision of the $\Lambda$CDM cosmology is the theory of non-
gravitational interaction between DE and DM where an energy exchange
phenomenon between these dark sectors, widely known as the interacting DE-DM,
also known as interacting DE (IDE) or Coupled DE (CDE) cosmology. This
particular theory received massive attention in the scientific community for
many interesting consequences [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46,
47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 76, 79, 80, 81, 82, 83,
84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96] (also see [97, 98, 99]),
such as, alleviating the cosmic coincidence problem [16, 20, 100, 101, 102,
103]; crossing the phantom divide line without invoking any scalar field with
negative sign in the kinetic part [104, 105, 106]; weakening/solving the
Hubble constant tension [107, 108, 63, 70, 109, 75] between Planck (within
$\Lambda$CDM paradigm) [110] and SH0ES (Supernovae and $H_{0}$ for the
Equation of State of dark energy) [111, 112]; and the clustering tension [113,
114, 109] between Planck (within $\Lambda$CDM model) and other astronomical
probes at low redshifts, e.g. weak gravitational lensing and galaxy clustering
[115, 116, 117, 118, 119, 120, 121, 122, 123, 124]. In IDE, the dynamics of
the dark sector is mainly governed by the choice of a coupling/interacting
function, $Q$, that controls the transfer of energy between DE and DM and this
coupling function is taken from the phenomenological ground111Some attempts
have been made to derive the coupling functions from an action integral [125,
126, 127, 128, 129, 79], however, the final destination is yet to be
discovered. Thus, at this moment, there is no reason to exclude any possible
approach to study the theory of DM-DE interaction, even the approach adopts a
phenomenological route.. Now, in the choice of the coupling functions, as they
represent the transfer of energy and/or momentum between the DM and DE
sectors, $Q$ is usually assumed to be the functions of the energy densities of
DE and DM. In general, two varieties of the interaction functions are
considered in the literature: (i) the interaction functions where the Hubble
rate, $H$, of the Friedmann-Lemaître-Robertson-Walker (FLRW) universe
explicitly appears, see for instance, [20, 21, 22]; (ii) the interaction
functions where $H$ does not appear explicitly, e.g. [29, 43, 130]. Concerning
the above two approaches, even though the interaction between these dark
sectors is viewed as a local phenomenon [29], and the presence of the global
expansion factor may be avoided, however, this debate is still unending and it
is very hard to prefer the first approach over the other (see [80]). On the
other hand, one can put another question mark on the direction of energy flow
between the dark sectors which is characterized by the choice of the
interaction function. In a large class of interaction models, the flow of
energy between the dark sectors is assumed to be unidirectional, that means,
throughout the period of energy exchange mechanism between the dark sectors,
the energy transfer can happen either from “DE to DM” or from “DM to DE”.
According to the theoretical and observational grounds, there are evidences of
the energy transfer from DE to DM [46, 70, 96], while the direction of energy
flow can be reversed in future [131], but this conclusion depends on the
underlying interaction model, properties of DE, and the observational data
[132, 92]222We note that in Refs. [132, 92], the authors have considered
various interaction models and constrained them using different datasets and
they reported that both the possibilities, that means the transfer of energy
from DE to DM and from DM to DE are allowed according to the observational
datasets. We further mention that the properties of DE (i.e. whether it is
quintessence or phantom) and the direction of energy transfer between the dark
components are connected with the stability of the interaction model at the
level of perturbations, see for instance, Ref. [83]. , hence, this is one of
the interesting questions in the context of interacting DE scenarios. Although
these unidirectional interacting scenarios are simple by construction and they
have been widely used in the community, however, it is very natural to examine
whether the direction of energy transfer may alter during the course of energy
exchange between the dark sectors. These kind of interaction models are known
as sign changeable or sign shifting interaction functions and such models are
appealing since they allow us to investigate whether the cosmologies with sign
changeable interaction models are physically viable. However, because of some
unknown reasons, sign changeable interaction models did not get much attention
in the community [133, 134, 135, 136, 137, 138, 139, 140], but such models are
worth investigating in the light of current cosmological tensions [139].
Interestingly, model independent inference on the interaction between the dark
sectors as performed in [141] hints for a sign shifting nature of the
interaction function. This gives enough motivation to allow a sign changeable
nature in the interaction functions and investigate the consequences.

In this article we therefore consider some sign shifting interaction models
where the interaction functions depend only on the intrinsic nature of the
dark fluids and perform their phase space analysis. As the choice of the
interaction functions are not unique, thus, we have considered a variety of
interaction functions that have been constructed using the known interaction
functions in the literature. On the other hand, as the nature of DE is another
unknown character to be discovered (hopefully) with the help of the upcoming
astronomical surveys, therefore, in this work, in order to be inclusive we
have considered that the equation of state of DE could be either constant or
dynamical. Now, focusing on the dynamical equation of state of DE, one may
have a cluster of possibilities since there is no unique route to determine
its expression. Keeping this issue in mind, we have adopted a very well known
dynamical equation of state of DE which depends only on its energy density and
having only one free parameter which characterizes the nature of the DE (i.e.
where DE is quintessential or phantom) through its sign. This equation of
state has been extensively investigated in the cosmological dynamics and it
recovers the usual barotropic equation of state of DE as a special case. So
far we are aware of the literature, the phase space analysis of the proposed
sign shifting interacting functions considering both the constant and
dynamical equation of state parameters of DE has never been performed in the
literature. This is the first time we are reporting the results in the
literature.

The article is structured as follows. In section II we introduce the basic
equations of an interacting DM-DE model and then propose the models of
interaction that we wish to study in this work. In section III we construct
the autonomous system corresponding to each interaction function and discuss
the nature of the critical points obtained from the interaction functions and
also their qualitative behaviour in terms of the cosmological parameters.
Finally, in section IV we close the article summarizing the key findings.

##  II Interacting Dark Energy

We consider the homogeneous and isotropic universe where its gravitational
sector is described by the Einstein’s General Relativity (GR) and its matter
distribution is minimally coupled to gravity. The matter sector consists of
two heavy fluids of the universe, namely, a pressure-less (or cold) DM and a
DE fluid which are interacting with each other in a non-gravitational way. In
order to proceed with the mathematical structure of such scenario, we consider
the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric

$\displaystyle ds^{2}=-{dt}^{2}+{a^{2}(t)}{d\bf{x}}^{2},$ (1)

where $t$ is the co-moving time; $a(t)$ is the expansion scale factor of the
universe; ${d\bf x}^{2}$ represents the 3-dimensional flat space line element.
The Friedmann equations for the above line element can be written as

$\displaystyle 3H^{2}={\kappa}^{2}(\rho_{c}+\rho_{d}),$ (2) $\displaystyle
2\dot{H}+3H^{2}=-{\kappa}^{2}(p_{c}+p_{d}),$ (3)

where an overhead dot denotes the derivative with respect to the cosmic time;
$\kappa^{2}=8\pi G$, is the Einstein’s gravitational constant;
$H\equiv\dot{a}(t)/a(t)$, is the Hubble rate of the FLRW universe;
${\rho_{d},~{}p_{d}}$ are respectively the energy density and pressure of the
DE fluid obeying the barotropic equation-of-state $w_{d}=p_{d}/\rho_{d}<-1/3$;
$\rho_{c}$, $p_{c}$ are respectively the energy density and pressure of DM in
the form of dust, i.e., $p_{c}=0$, henceforth, we call this DM as cold DM,
abbreviated as CDM. As CDM and DE are interacting with each other, therefore,
the conservation equations of these dark fluids can be represented as

$\displaystyle\dot{\rho}_{c}+3H\rho_{c}=-Q(\rho_{c},\rho_{d}),$ (4)
$\displaystyle\dot{\rho}_{d}+3H(1+w_{d})\rho_{d}=+Q(\rho_{c},\rho_{d}),$ (5)

where $Q(\rho_{c},\rho_{d})$ denotes the real valued interaction function
(also known as the interaction rate) that corresponds to the transfer of
energy and (or) momentum between these dark fluids. For
$Q(\rho_{c},\rho_{d})>0$, energy flow occurs from DM to DE, and for
$Q(\rho_{c},\rho_{d})<0$, energy flow occurs in the reverse direction, that
means from DE to DM. The interaction function $Q(\rho_{c},\rho_{d})$ is the
key ingredient of this scenario because it controls the dynamics of the
universe at the background and perturbation levels. We notice that the
conservation equations (4) and (5) can be put in a different format leading to

$\displaystyle\dot{\rho}_{c}+3H\left(1+w_{\rm c,eff}\right)\rho_{c}=0,$ (6)
$\displaystyle\dot{\rho}_{d}+3H\left(1+w_{\rm d,eff}\right)\rho_{d}=0,$ (7)

which represent a non-interacting scenario of DM and DE with the effective
equation-of-state parameters

$\displaystyle w_{\rm c,eff}=\frac{Q(\rho_{c},\rho_{d})}{3H\rho_{c}},\quad
w_{\rm d,eff}=w_{d}-\frac{Q(\rho_{c},\rho_{d})}{3H\rho_{d}},$ (8)

from which one can notice that the effective nature of the DM equation of
state could be non-cold in the sense that the effective equation of state of
DM could be non-zero (i.e. $w_{\rm c,eff}\neq 0$) for
$Q(\rho_{c},\rho_{d})\neq 0$, see for instance [142], and additionally, the
effective nature of the DE equation of state could be either quintessential
($w_{\rm d,eff}>-1$) or phantom ($w_{\rm d,eff}<-1$) depending on the sign of
$Q(\rho_{c},\rho_{d})$.

Now since the interaction function affects the evolution of both CDM and DE,
henceforth, the expansion rate of the universe $H$ will be equally affected
and as a result the cosmological parameters will be influenced as well. We
introduce the equation of state of the total fluid $w_{\rm tot}=\frac{{\rm
total~{}pressure}}{{\rm
total~{}energy~{}density}}=\frac{p_{c}+p_{d}}{\rho_{c}+\rho_{d}}=\frac{p_{d}}{\rho_{c}+\rho_{d}}$
(since $p_{c}=0$) and the deceleration parameter of the universe,
$q=-(1+\dot{H}/H^{2})$, which take the following forms

$\displaystyle w_{\rm tot}=w_{d}\Omega_{d},\quad\quad
q=\frac{1}{2}(1+3w_{d}\Omega_{d}),$ (9)

where $\Omega_{d}=\kappa^{2}\rho_{d}/3H^{2}$ is the density parameter for DE
and from the Friedmann equation (2), one can derive the density parameter for
CDM, $\Omega_{c}~{}(=\kappa^{2}\rho_{c}/3H^{2})$ as $\Omega_{c}=1-\Omega_{d}$.

###  II.1 Models

In this work we propose several interaction functions having the sign changing
property during the evolution of the universe. One of the important features
of all the interaction functions that we are going to propose in this section
is that, all of them do not depend on the external parameters of the universe,
rather they all depend on the intrinsic nature of the dark sector. The first
interaction function in this series has the following form

$\displaystyle\mbox{Model I:}\quad Q_{\rm I}=\Gamma(\rho_{c}-\rho_{d}),$ (10)

where $\Gamma$ is the coupling parameter of the interaction function measuring
the strength of the interaction and it has the dimension of the Hubble rate
$H$. As argued, $Q$ does not depend on the external parameters of the
universe, e.g. the scale factor of the universe or its expansion rate, rather
it depends on the intrinsic properties of DM and DE, namely, their energy
densities, $\rho_{c}$ and $\rho_{d}$. Hence, one may expect that this
interaction model could offer some inherent nature of the dark components. We
further note that $Q_{\rm I}$ may change its sign depending on the dominating
role played by one of the fluids, that means, if the dominating role played by
DM over DE ($\rho_{c}>\rho_{d}$) is suddenly altered, i.e., DE starts
dominating over DM as in the late time accelerating phase
($\rho_{d}>\rho_{c}$), then $Q_{\rm I}$ shifts its sign.

We generalize the sign-shifting interaction function of eqn. (10) as follows

$\displaystyle\mbox{Model II:}\quad Q_{\rm
II}=\Gamma_{c}\rho_{c}-\Gamma_{d}\rho_{d},$ (11)

where $\Gamma_{c}$ and $\Gamma_{d}$ are the coupling parameters of the
interaction function having the dimension equal to the dimension of the Hubble
parameter. For the interaction function $Q_{\rm II}$ to be sign changeable,
both the coupling parameters $\Gamma_{c}$, $\Gamma_{d}$ should have the same
sign, that means either $\Gamma_{c}>0$, $\Gamma_{d}>0$ or $\Gamma_{c}<0$,
$\Gamma_{d}<0$, but never $\Gamma_{c}\Gamma_{d}<0$. Similar to the earlier
interaction function, this model also does not include any external parameters
of the universe.

The next model in this series that we introduce has the following form

$\displaystyle\mbox{Model III:}\quad Q_{\rm
III}=\Gamma\left(\rho_{c}-\rho_{d}-\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}\right),$
(12)

where as already noted, $\Gamma$ is the coupling parameter of the interaction
function having the dimension equal to the dimension of the Hubble parameter.
Notice that the interaction function (12) is obtained by including a new
function $Q_{\rm new}=-\Gamma\rho_{c}\rho_{d}(\rho_{c}+\rho_{d})^{-1}$ with
the model of eqn. (10), that means $Q_{\rm III}=Q_{\rm I}+Q_{\rm new}$.
Similar to the earlier two interaction functions, here too, we notice that
this interaction function depends only on the intrinsic nature of the dark
fluids.

Lastly, we introduce two new interactions of the forms

$\displaystyle\mbox{Model IV:}\quad Q_{\rm
IV}=\Gamma_{c}\rho_{c}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}},$
(13)

and

$\displaystyle\mbox{Model V:}\quad Q_{\rm
V}=\Gamma_{d}\rho_{d}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}},$
(14)

where $\Gamma_{c}$, $\Gamma_{d}$ and $\Gamma_{cd}$ are the coupling parameters
and they all have the dimension of the Hubble parameter. One can easily notice
that the last term of both Model IV (eqn. (13)) and Model V (eqn. (14)) are
same but the models differ in their first terms containing $\rho_{c}$ and
$\rho_{d}$, respectively. We further note that $\Gamma_{c}$, $\Gamma_{d}$ and
$\Gamma_{cd}$ are all constants in such a way so that the models allow sign
shifting property. That means, for Model IV (eqn. (13)), $\Gamma_{c}$ and
$\Gamma_{cd}$ will enjoy the same sign and for Model V (eqn. (14)),
$\Gamma_{d}$ and $\Gamma_{cd}$ will enjoy the same sign but
$\Gamma_{c}\neq\Gamma_{cd}$ for Model IV and $\Gamma_{d}\neq\Gamma_{cd}$ for
Model V, otherwise $Q_{\rm IV}$ and $Q_{\rm V}$ will represent the energy flow
only in one direction. As the choice of the interaction function is not
unique, therefore, one can construct a variety of such models, however, it
should be kept in mind that the interaction functions may lead to negative
energy densities of the dark sectors as argued in [130], hence, working with
an arbitrary interaction function needs precaution.

##  III Dynamical Analysis of Interacting models

The dynamical analysis of the interaction models is the heart of this work.
The dynamical analysis plays a crucial role in understanding the local and
global dynamics of the underlying cosmological scenarios. In order to perform
the dynamical analysis of the underlying interacting scenarios, one needs to
define a new set of dimensionless variables in terms of which one can study
the behaviour of the system of differential equations. We refer [143, 144,
145, 146, 27, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
160, 161, 162, 163, 164, 165] (also see the review article [166] and the
references therein) where several interacting scenarios have been studied
through the dynamical system analysis. In this section we shall describe the
dynamical systems for the proposed sign-shifting interaction models.
Additionally, we shall also show that the choice of the dimensionless
variables is extremely important because for a wrong choice of these
variables, one may not be able explore the whole space of critical points.

###  III.1 Model I

We begin our analysis with the first interaction function $Q_{\rm I}$ of eqn.
(10) and we define the following dimensionless variables

$\displaystyle x=\frac{\kappa^{2}\rho_{c}}{3H^{2}},\quad\quad
y=\frac{H_{0}}{H},$ (15)

where $H_{0}$ ($>0$) is a constant and here it denotes the present value of
the Hubble parameter.333Note that instead of taking $H_{0}$, one may consider
any $\widetilde{H}$ which is also constant so that $y=\widetilde{H}/H$ becomes
dimensionless.  The density parameters for CDM, DE, the equation state of the
total fluid, $w_{\rm tot}$, and the deceleration parameter $q$ can be
expressed in terms of the dimensionless parameters as follows

$\displaystyle\Omega_{c}=x,\quad\quad\quad\quad\quad\quad\Omega_{d}=1-x,$
$\displaystyle w_{\rm tot}=w_{d}(1-x),\quad\quad
q=\frac{1}{2}[1+3w_{d}(1-x)].$

Now, inserting the dimensionless variables of eqn. (15) into the equations of
motion (2), (3), (4), and using the interaction function $Q_{\rm I}$ of eqn.
(10), we have the following autonomous system:

$\displaystyle x^{\prime}=-\gamma y(2x-1)+3w_{d}x(1-x),$ (16) $\displaystyle
y^{\prime}=\frac{3}{2}y\left(1+w_{d}(1-x)\right),$ (17)

where a prime over a variable denotes its differentiation with respect to
$N=\ln(a/a_{0})$ (here $a_{0}$ is the present value of the scale factor) and
$\gamma$ is defined as $\gamma=\Gamma/H_{0}$, hence, $\gamma$ becomes a
dimensionless parameter. The critical points of the autonomous system (16) $-$
(17) are obtained by solving the equations $x^{\prime}=0$ and $y^{\prime}=0$.
In this case we have three critical points $A_{1}=(0,0)$, $A_{2}=(1,0)$ and
$A_{3}=\left(\frac{1+w_{d}}{w_{d}},-\frac{3(1+w_{d})}{\gamma(2+w_{d})}\right)$.

In order to better understand the dynamical system when the Hubble rate is
small, it will be useful to use the variables $(x,H)$. Now, using the cosmic
time derivative, the dynamical system for the variables $(x,H)$ becomes

$\displaystyle\left\\{\begin{array}[]{ccc}\dot{x}&=&-\Gamma(2x-1)+3w_{d}Hx(1-x),\\\
\dot{H}&=&-\frac{3}{2}H^{2}(1+w_{d}(1-x)).\end{array}\right.$ (20)

On can see that the autonomous system (20) has only two critical points,
namely, $A_{0}=(\frac{1}{2},0)$ and
$\left(\frac{1+w_{d}}{w_{d}},-\frac{\Gamma(2+w_{d})}{3(1+w_{d})}\right)$ which
is actually $A_{3}$ if we consider this critical point in terms of $(x,y)$
coordinates. Here, it is important to understand that the point $A_{0}$ does
not appear in the coordinates $(x,y)$, because it corresponds to $y=\infty$,
and for the same reason the points $A_{1}$ and $A_{2}$ do not appear in the
coordinates $(x,H)$, because they correspond to $H=\infty$. This clearly
indicates that the above two autonomous systems are not giving the complete
information.

For this reason it is important to get a new coordinate system where all the
critical points appear. This could be done by introducing a new variable

$\displaystyle z=\frac{y}{1+y}=\frac{H_{0}}{H+H_{0}},$ (21)

where $H=0$ corresponds to $z=1$ and $H=\infty$ corresponds to $z=0$. So, in
the new coordinates $(x,z)$ the dynamical system (16) $-$ (17) becomes

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma\left(\frac{z}{1-z}\right)(2x-1)+3w_{d}x(1-x),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z(1+w_{d}(1-x)).\end{array}\right.$ (24)

But, note that, the dynamical system (24) is singular at $z=1$ (i.e., $H=0$).
Thus, in order to regularize it, and taking into account that the phase
portrait does not change topologically when one multiplies the vector field by
a positive function, we regularize it by multiplying the factor $(1-z)$, which
leads to the regular dynamical system

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma{z}(2x-1)+3(1-z)w_{d}x(1-x),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z(1+w_{d}(1-x)).\end{array}\right.$ (27)

Thus, finally one can investigate the dynamical system (27) through the
stability analysis of the critical points. The physical domain, namely $R$, is
the square $R=[0,1]^{2}$, which in order to be positively invariant, i.e., to
ensure that a solution of the dynamical system with initial conditions in $R$
never leaves it, one has to impose that $\gamma>0$ (i.e., $\Gamma>0$).
Effectively, it follows from the equation (24) that close to $x=0$, one has
$x^{\prime}\cong\gamma z>0$ and close to $x=1$, one has
$x^{\prime}\cong-\gamma z<0$. For $z=0$ and $z=1$, one has $z^{\prime}=0$,
meaning that, for $\gamma>0$, the dynamical system never crosses the lines
$x=0$, $x=1$, $z=0$ and $z=1$, that is, the domain $R$ is positively
invariant. As the dynamical system contains the DE equation of state which
could be either constant or dynamical, therefore, we aim to investigate both
the cases separately. In the following we present our analyses for constant
$w_{d}$ and dynamical $w_{d}$.

####  III.1.1 Constant $w_{d}$

Considering that the DE equation of state, $w_{d}$, is constant, in Table 1 we
present the critical points of the dynamical system (27), their existence,
stability, and the values of the cosmological parameters evaluated at those
critical points. Now depending on the nature of $w_{d}$, it is comprehensible
that the nature of the critical points will certainly be affected. Thus, in
order to be precise, we divide the entire parameter space of $w_{d}$ into
three disjoint regions, namely, quintessence or non-phantom ($w_{d}>-1$),
cosmological constant ($w_{d}=-1$) and phantom ($w_{d}<-1$). In the following
we present how $w_{d}$ affects the nature of the critical points.

  1. 1.

For non-phantom dark energy, i.e., when $w_{d}>-1$, the point $A_{3}$ does not
belong to the physical domain $R$. In addition, since $w_{\rm
tot}=w_{d}(1-x)$, this means that $w_{\rm tot}>-1$ and taking into account
that $z^{\prime}=\frac{3}{2}(1-z)^{2}z(1+w_{d}(1-x))$, we have $z^{\prime}>0$.
So, $A_{1}$ and $A_{2}$ are unstable. On the line $z=1$, we have
$z^{\prime}=0$, $x^{\prime}<0$ for $x>1/2$ and $x^{\prime}>0$ for $x<1/2$.
Thus, the critical point $A_{0}$ is stable, in fact, it is a global attractor.
Therefore, at late times, the universe accelerates with $w_{\rm
tot}=w_{d}/2>-1$, when $-1<w_{d}<-2/3$, and it decelerates when
$-2/3<w_{d}<-1/3$. In addition, since $w_{\rm tot}>-1$ the Hubble rate
decreases to zero, and $\Omega_{d}=\Omega_{c}=1/2$. The phase plot is shown in
Fig. 1. This finishes the study for a non-phantom dark energy.

  2. 2.

When $w_{d}=-1$, one has $A_{1}=A_{3}$. In this case, the linearization does
not decide the nature of the critical point $A_{1}$. In fact, $A_{1}$ and
$A_{2}$ are unstable because $z^{\prime}=\frac{3}{2}(1-z)^{2}zx>0$, for
$0<z<1$ and $x>0$. Moreover, only the unphysical orbit $z=0$ ($H=\infty$)
converges to $A_{1}$. Now, on $z=1$ line, we obtain $z^{\prime}=0$,
$x^{\prime}<0$ for $x>1/2$ and $x^{\prime}>0$ for $x<1/2$. Consequently, in
this case $A_{0}$ is a global attractor. Again, the phase plot is given in
Fig. 1.

  3. 3.

For a phantom dark fluid, two cases arise:

    1. (a)

When $-2<w_{d}<-1$: The critical point $A_{3}$ belongs to the physical region
$R$. Clearly, we obtain $0<(1+w_{d})/w_{d}<1/2$. Whenever $x<(1+w_{d})/w_{d}$,
we get $1+w_{\rm tot}<0$ and so, $z^{\prime}$ is negative. Similarly, we have
$1+w_{\rm tot}>0$ for $x>(1+w_{d})/w_{d}$ and as a result $z^{\prime}$ is
positive. On $z=1$ line, we obtain $z^{\prime}=0$, $x^{\prime}<0$ for $x>1/2$
and $x^{\prime}>0$ for $x<1/2$. Again, on $z=0$ line, $z^{\prime}=0$ and
$x^{\prime}=3w_{d}x(1-x)$ which is negative. Therefore, our domain $R$ is
divided into four regions. Thus, an orbit in the regions I and II of Fig. 2,
at late times, converges to $A_{1}$. For an orbit in the regions III and IV,
at late time, it converges to $A_{0}$. Note that $A_{0}$ means $H=0$ with
$\Omega_{c}=\Omega_{d}=1/2$ and $w_{\rm tot}=w_{d}/2>-1$. On the contrary,
$A_{1}$ means $H=\infty$ with $\Omega_{d}=1$ and $w_{\rm tot}=w_{d}<-1$. Fig.
3 shows the evolution of the density parameters, namely, $\Omega_{c}$,
$\Omega_{d}$, and the total equation of state parameter, $w_{\rm tot}$.

    2. (b)

When $w_{d}\leq-2$: For the particular case with $w_{d}=-2$, which represents
a very high phantom regime, $A_{3}=A_{0}$. Here, $(1+w_{d})/w_{d}=1/2$, so we
have $z^{\prime}<0$ for $x<1/2$ and $z^{\prime}>0$ for $x>1/2$. On the line
$z=0$, we obtain $z^{\prime}=0$ and $x^{\prime}<0$. Therefore, $A_{0}$,
$A_{2}$ are unstable and $A_{1}$ becomes a global attractor. When $-2>w_{d}$,
the critical point $A_{3}$ does not belong to the physical region $R$. Again,
$(1+w_{d})/w_{d}>1/2$ which implies that $z^{\prime}<0$ if $x<(1+w_{d})/w_{d}$
and $z^{\prime}>0$ if $x>(1+w_{d})/w_{d}$. Thus, $A_{0}$, $A_{2}$ are unstable
and $A_{1}$ is a global attractor. Fig. 4 exhibits the behavior.

Point |  $x$ |  $z$ |  Existence |  Stability |  Acceleration |  $\Omega_{c}$ |  $\Omega_{d}$ |  $w_{\rm tot}$
---|---|---|---|---|---|---|---|---
$A_{0}$ | $\frac{1}{2}$ | 1 | for all $\gamma>0$ and $w_{d}<-1/3$ | stable for $-2<w_{d}<-1/3$ | $w_{d}<-2/3$ | $1/2$ | $1/2$ | $w_{d}/2$
$A_{1}$ | $0$ | $0$ | for all $\gamma>0$ and $w_{d}<-\frac{1}{3}$ | stable for $w_{d}<-1$ | $w_{d}<-\frac{1}{3}$ | $0$ | $1$ | $w_{d}$
$A_{2}$ | $1$ | $0$ | for all $\gamma>0$ and $w_{d}<-\frac{1}{3}$ | unstable | no | $1$ | $0$ | $0$
$A_{3}$ | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{3(1+w_{d})}{2\gamma-3+w_{d}(\gamma-3)}$ | for all $\gamma>0$ and $-2\leq w_{d}\leq-1$ | unstable | yes | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{1}{w_{d}}$ | $-1$
Table 1: The critical points, their existence, stability, and the values of
the cosmological parameters evaluated at those points for the interacting
scenario driven by the interaction function $Q_{\rm
I}=\Gamma(\rho_{c}-\rho_{d})$ of eqn. (10) are summarized.  Figure 1: The
phase portrait plot describing Model I (eqn. (10)) with $w_{d}\geq-1$ and
$\gamma>0$. In this case we have taken $w_{d}=-0.95$ and $\gamma=0.4$. We note
that one can take any value of $w_{d}\geq-1$ and any positive value of
$\gamma$ in order to get similar graphics. Here, the yellow shaded region
represents the accelerated region (i.e. $q<0$) and the pink shaded region
corresponds to the decelerated region (i.e. $q>0$).  Figure 2: The phase
portrait plot describing Model I (eqn. (10)) for $-2<w_{d}<-1$ and $\gamma>0$.
In this case we have taken $w_{d}=-1.3$ and $\gamma=0.8$. We note that one can
take any specific value of $\gamma~{}(>0)$ to draw the plot, however, as long
as $\gamma$ decreases, the regions I and IV become very small and they look
indistinguishable from one another.  Figure 3: We show the evolution of the
CDM density parameter ($\Omega_{c}$), dark energy density parameter
($\Omega_{d}$) and the total equation of state (EoS) parameter ($w_{\rm tot}$)
for Model I (eqn. (10)) for $-2<w_{d}<-1$. We have taken $w_{d}=-1.3$,
$\gamma=0.8$ with the initial conditions $x~{}(N=0)=0.25$, $z~{}(N=0)=0.05$
taken from region II of Fig. 2. For the initial condition on $x(N)$ and $z(N)$
from the region I of Fig. 2, again we shall obtain $\Omega_{c}=0$ and
$\Omega_{d}=1$ at late time. If we choose the initial conditions on $x(N)$ and
$y(N)$ from region III and IV of Fig. 2, we shall reach
$\Omega_{c}=\Omega_{d}=1/2$ in an asymptotic way.  Figure 4: The phase
portrait plot describing Model I of (10) with $w_{d}\leq-2$ and $\gamma>0$. In
this case we have taken $w_{d}=-2.5$ and $\gamma=0.8$. We note that one can
take any value of $w_{d}\leq-2$ and any positive value of $\gamma$ in order to
get similar graphics. Here, the yellow shaded region represents the
accelerated region (i.e. $q<0$) and the pink shaded region corresponds to the
decelerated region (i.e. $q>0$)

####  III.1.2 Dynamical $w_{d}$

The case with dynamical equation of state of DE, $w_{d}$, is interesting for
two reasons. First of all, the autonomous system (27) with constant $w_{d}$ is
a special case of the dynamical $w_{d}$ case. On the other hand and most
importantly, the dynamical $w_{d}$ scenario may offer a bigger space of
critical points and hence one may expect new results in this context. In this
article, we shall consider a parametric form for $w_{d}$ to investigate the
autonomous system (27). The choice of the dynamical $w_{d}$ is not unique, and
one can consider a variety of choices. The question then arises, what should
be a possible choice for $w_{d}$ to proceed with the analysis? A possible
choice for dynamical $w_{d}$ capturing a wide variety of models in this
direction may take the following form [167, 168]

$\displaystyle p_{d}=-\rho_{d}-f(\rho_{d}),$ (28)

where $f$ is any analytic function of $\rho_{d}$. Notice from (28) that,
$w_{d}=p_{d}/\rho_{d}=-1-f(\rho_{d})/\rho_{d}$, describes a deviation from the
cosmological constant $w_{d}=-1$ through the dynamical term
$f(\rho_{d})/\rho_{d}$. A general choice of this equation of state could be
$p_{d}=-\rho_{d}-A\rho_{d}^{n}$, where $n$ and $A$ are constants in which $n$
is a dimensionless constant but $A$ has dimension. We restrict ourselves to
$n=2$ in this article for which we have

$\displaystyle p_{d}=-\rho_{d}-A\rho_{d}^{2}\quad\Longleftrightarrow\quad
w_{d}=-1-A\rho_{d}.$ (29)

Now, for the equation of state (29), the autonomous system (24) becomes,

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma\left(\frac{z}{1-z}\right)(2x-1)-3x(1-x)\left(1+\nu\frac{(1-z)^{2}(1-x)}{z^{2}}\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z\left(x-\nu\frac{(1-z)^{2}(1-x)^{2}}{z^{2}}\right),\end{array}\right.$
(32)

where $\nu=\frac{3AH_{0}^{2}}{\kappa^{2}}$. We regularize the autonomous
system (32) by multiplying the factor $z^{2}(1-z)$ on the right hand sides of
(32) and finally obtain

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma{z^{3}}(2x-1)-3(1-z)x(1-x)\left(z^{2}+\nu(1-z)^{2}(1-x)\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z\left(xz^{2}-\nu(1-z)^{2}(1-x)^{2}\right).\end{array}\right.$
(35)

Note that the qualitative behaviour for both the autonomous systems (32) and
(35) remain topologically equivalent. From the autonomous system (35), one can
now find the critical points by solving the equations $x^{\prime}=0$ and
$z^{\prime}=0$ and the critical points of the system (35) are

  * •

$\bar{A}_{0}=(\frac{1}{2},1)$, $\bar{A}_{1}=(0,0)$, $\bar{A}_{2}=(1,0)$,
$S=\left\\{\left(x_{c},\frac{3x_{c}}{3x_{c}-\gamma(2x_{c}-1)}\right)\right\\}$,

where $S$ is the set of critical points in which $x_{c}$ denotes a real root
of $f(x)\equiv 9x^{3}-\nu{\gamma}^{2}(2x-1)^{2}(1-x)^{2}=0$.444We note that
$f(x)$ can be obtained from the following two nullclines:  $\displaystyle
xz^{2}-\nu(1-z)^{2}(1-x)^{2}=0,$ (36) $\displaystyle-\gamma
z^{3}(2x-1)-3x(1-x)(1-z)\left[z^{2}+\nu(1-z)^{2}(1-x)\right]=0.$ (37) Since,
$f(x)$ is a fourth degree equation in $x$, the set $S$ may contain maximum 4
critical points.

Now, as the physical domain in our case is $R=[0,1]^{2}$, therefore, we are
interested to investigate the number of roots of $f(x)$ in $[0,1]$. As
$f(0)=-{\gamma}^{2}\nu<0$ (for $\nu>0$) and $f(1)=9>0$, hence, by the
Bolzano’s theorem555Bolzano’s theorem: Let $f$ be a real valued and continuous
function in a compact interval $[k,l]$ in $\mathbb{R}$ and suppose that
$f(k)$, $f(l)$ have opposite signs, that means $f(k)f(l)<0$. Then there is at
least one point $m$ in $(k,l)$ such that $f(m)=0$. [169], $f(x)$ will have at
least one real root in $(0,1)$. Here we argue that $f$ will have only one real
root in $(0,1)$. It also follows that since $0<x_{c}<1$, therefore, the
condition $z=\frac{3x_{c}}{3x_{c}-\gamma(2x_{c}-1)}\leq 1$ leads to $x_{c}\leq
1/2$. Hence, our domain is slightly reduced and we need to check the number of
roots of $f(x)$ in $[0,1/2]$. Since, we have $f(1/2)=9/8>0$, it follows from
the Bolzano’s theorem [169], that there is at least one real root of $f(x)$ in
$(0,1/2)$. Now, looking at the derivative of $f(x)$ with respect to $x$, given
by
$f^{\prime}(x)=-16{\gamma}^{2}\nu(x-1)\left(x-\frac{3}{4}\right)\left(x-\frac{1}{2}\right)+27x^{2}$,
one can check that $f^{\prime}(x)>0$ for all $x\in(0,1/2)$. This shows that
the function $f(x)$ is strictly increasing in $(0,1/2)$ and consequently
$f(x)$ has only one root in $(0,1/2)$. This concludes that the set $S$ has
only one critical point in the physical domain $R$ and we label this critical
point as $\bar{A}_{3}$. Now, in this case, we observe that the point
$\bar{A}_{3}$ behaves qualitatively same as the point $A_{3}$ described
earlier and correspondingly, the phase portrait is same as Fig. 2.

On the other hand, for $\nu<0$, the algebraic curve represented by the eqn.
(36) has no branches in the positive quadrant. Hence, for $\nu<0$, $S$ is an
empty set, and as a result, the autonomous system for $\nu<0$ admits only
three critical points, namely, $\bar{A}_{0}$, $\bar{A}_{1}$, $\bar{A}_{2}$.
The corresponding phase plot will be similar to Fig. 1.

###  III.2 Model II

We now consider the second interaction model in this series, i.e. $Q_{\rm II}$
of (11). Notice that $Q_{\rm II}$ has two coupling parameters $\Gamma_{c}$ and
$\Gamma_{d}$. As demonstrated in section III.1, here we shall consider the
following dimensionless variables

$\displaystyle x=\frac{\kappa^{2}\rho_{c}}{3H^{2}},\quad\quad
z=\frac{H_{0}}{H+H_{0}},$ (38)

in order to understand the dynamics of the interacting scenario. With these
choice of the dynamical variables, the autonomous system for this interaction
function can be expressed as

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\frac{z}{1-z}[\alpha
x-\beta(1-x)]+3w_{d}x(1-x),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z(1+w_{d}(1-x)),\end{array}\right.$ (41)

where $\alpha,\beta$ are the dimensionless parameters defined as
$\alpha=\Gamma_{c}/H_{0}$, $\beta=\Gamma_{d}/H_{0}$. Now, regularizing (41) we
get

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-z[\alpha
x-\beta(1-x)]+3(1-z)w_{d}x(1-x),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z(1+w_{d}(1-x)).\end{array}\right.$ (44)

The physical domain, namely $R$, is the square $R=[0,1]^{2}$, and it follows
from the autonomous system (44) that along the lines $x=0$ and $x=1$, one has
$x^{\prime}=\beta z$ and $x^{\prime}=-\alpha z$, and also the lines $z=0$ and
$z=1$ remain invariant. Thus, to ensure that the physical domain $R$ is
positively invariant we have to restrict our attention on $\alpha>0$ and
$\beta>0$ (i.e., $\Gamma_{c}>0$ and $\Gamma_{d}>0$), and therefore, $Q_{\rm
II}$ may allow a sign change during the evolution of the universe without
exhibiting any unphysical properties in the energy densities of the dark
sector.

Now, in a similar fashion we focus on two cases, namely, the constant $w_{d}$
and dynamical $w_{d}$. In the following we consider both the possibilities.

####  III.2.1 Constant $w_{d}$

Considering $w_{d}$ as a constant, in Table 2 we summarize the critical points
of the autonomous system (44), their existence, stability and as well as the
cosmological parameters evaluated at those critical points. Now we consider
three different regions of $w_{d}$ as follows: if $w_{d}$ has quintessential
nature (i.e. $w_{d}>-1$); if $w_{d}$ mimics a cosmological constant (i.e.
$w_{d}=-1$); if $w_{d}$ has a phantom character ($w_{d}<-1$). In what follows
we investigate each case.

  1. 1.

When $w_{d}>-1$, the point $B_{3}$ does not belong to the physical domain $R$.
Now, we can see that on $z=1$ line, we have $z^{\prime}=0$, $x^{\prime}>0$ for
$x<\frac{\beta}{\alpha+\beta}$ and $x^{\prime}<0$ for
$x>\frac{\beta}{\alpha+\beta}$. Also, we have $w_{\rm tot}=w_{d}(1-x)>-1$
which implies $z^{\prime}=\frac{3}{2}z(1-z)^{2}(1+w_{\rm tot})>0$. So,
$B_{1}$, $B_{2}$ are unstable critical points and $B_{0}$ is a global
attractor. Note that, $B_{0}$ corresponds to $H=0$,
$\Omega_{c}=\frac{\beta}{\alpha+\beta}$ and
$\Omega_{d}=\frac{\alpha}{\alpha+\beta}$. The qualitative behaviour is
displayed in Fig. 5.

  2. 2.

When $w_{d}=-1$, we have $B_{3}=B_{1}$. Again, in this case, on the line
$z=1$, one has $z^{\prime}=0$, $x^{\prime}>0$ for
$x<\frac{\beta}{\alpha+\beta}$ and $x^{\prime}<0$ for
$x>\frac{\beta}{\alpha+\beta}$. Noting that
$z^{\prime}=\frac{3}{2}z(1-z)^{2}x$ is positive. Thus, $B_{1}$, $B_{2}$ are
unstable, and $B_{0}$ is a global attractor. Again, Fig. 5 shows the phase
plot.

  3. 3.

Now dealing with a phantom dark energy fluid, we have the following results:

    1. (a)

When $-1-\frac{\beta}{\alpha}<w_{d}<-1$, the point $B_{3}$ belongs to the
physical domain and we get
$0<\frac{1+w_{d}}{w_{d}}<\frac{\beta}{\alpha+\beta}$. Also, we have
$z^{\prime}<0$ for $x<\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}>0$ for
$x>\frac{1+w_{d}}{w_{d}}$. On $z=0$ line, $x^{\prime}$ is negative. Again, on
the line $z=1$, $x^{\prime}$ is positive for $x<\frac{\beta}{\alpha+\beta}$
and $x^{\prime}$ is negative for $x>\frac{\beta}{\alpha+\beta}$. Thus, we have
two “invariant stable orbits” (see Fig. 6) which divide the physical domain in
two parts. The orbits below these “invariant manifolds” converge to $B_{1}$
and the orbits above them converge to $B_{0}$. Fig. 7 displays the evolution
of the density parameters, namely, $\Omega_{c}$, $\Omega_{d}$, and the total
equation of state parameter, $w_{\rm tot}$.

    2. (b)

When $w_{d}=-1-\frac{\beta}{\alpha}$, we have $B_{0}=B_{3}$ and
$\frac{1+w_{d}}{w_{d}}=\frac{\beta}{\alpha+\beta}$. For
$x<\frac{1+w_{d}}{w_{d}}$, one has $z^{\prime}<0$ and for
$x>\frac{1+w_{d}}{w_{d}}$, $z^{\prime}$ is positive. Since, $x^{\prime}<0$ on
$z=0$, $B_{1}$ is a global attractor. The qualitative behavior is given in
Fig. 8.

    3. (c)

When $w_{d}<-1-\frac{\beta}{\alpha}$, the point $B_{3}$ does not belong to the
physical domain $R$. Here, we obtain
$\frac{1+w_{d}}{w_{d}}>\frac{\beta}{\alpha+\beta}$, $z^{\prime}<0$ if
$x<\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}>0$ if $x>\frac{1+w_{d}}{w_{d}}$. On
$z=0$ line, we have $x^{\prime}=3w_{d}x(1-x)$ which is negative. Therefore,
$B_{0}$, $B_{1}$ are unstable critical points and once again $B_{1}$ is a
global attractor. Fig. 8 exhibits the nature of phase portrait.

Point | $x$ | $z$ | Existence | Stability | Acceleration | $\Omega_{c}$ | $\Omega_{d}$ | $w_{\rm tot}$
---|---|---|---|---|---|---|---|---
$B_{0}$ | $\frac{\beta}{\alpha+\beta}$ | 1 |  $\alpha~{}(>0),\beta~{}(>0)$ and $w_{d}<-\frac{1}{3}$ | $-1-\frac{\beta}{\alpha}<w_{d}<-\frac{1}{3}$ | $w_{d}<-\frac{1}{3}(1+\frac{\beta}{\alpha})$ | $\frac{\beta}{\alpha+\beta}$ | $\frac{\alpha}{\alpha+\beta}$ | $\frac{\alpha w_{d}}{\alpha+\beta}$
$B_{1}$ | $0$ | $0$ |  $\alpha~{}(>0),\beta~{}(>0)$ and $w_{d}<-\frac{1}{3}$ | $w_{d}<-1$ | $w_{d}<-\frac{1}{3}$ | $0$ | $1$ | $w_{d}$
$B_{2}$ | $1$ | $0$ |  $\alpha~{}(>0),\beta~{}(>0)$ and $w_{d}<-\frac{1}{3}$ | unstable | no | $1$ | $0$ | $0$
$B_{3}$ | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{3(1+w_{d})}{(\alpha-3)w_{d}+(\alpha+\beta-3)}$ |  $\alpha~{}(>0),\beta~{}(>0)$ and $-1-\frac{\beta}{\alpha}\leq w_{d}\leq-1$ | unstable | yes | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{1}{w_{d}}$ | $-1$
Table 2: The critical points, their existence, stability, and the values of
the cosmological parameters evaluated at those points for the interacting
scenario driven by the interaction function $Q_{\rm
II}=\Gamma_{c}\rho_{c}-\Gamma_{d}\rho_{d}$ of eqn. (11) are summarized.
Figure 5: Phase portrait plot depicting Model II (eqn. (11)) with
$\alpha>0,\beta>0$ and $w_{d}\geq-1$. In this case we have taken $w_{d}=-0.9$,
$\alpha=0.5$ and $\beta=0.6$. We note that one can take any value of
$w_{d}\geq-1$ and any positive value of $\alpha$ and $\beta$ in order to get
similar graphics. Here, the yellow shaded region represents the accelerated
region (i.e. $q<0$) and the pink shaded region corresponds to the decelerated
region (i.e. $q>0$).  Figure 6: Phase plot for Model II (eqn. (11)) with
$\alpha>0,\beta>0$ and $-1-\frac{\beta}{\alpha}<w_{d}<-1$. Here, we have taken
$w_{d}=-1.4$, $\alpha=0.8$ and $\beta=0.9$. Note that one can take any
specific value of $\alpha~{}(>0)$ and $\beta~{}(>0)$ to draw the plot,
however, as long as $\alpha$ and $\beta$ decrease, the regions I and IV become
very small and they look indistinguishable from one another.  Figure 7: We
show the evolution of the CDM density parameter ($\Omega_{c}$), dark energy
density parameter ($\Omega_{d}$) and the total equation of state (EoS)
parameter ($w_{\rm tot}$) for the Model II (eqn. (11)). We have taken the
following values of the parameters: $w_{d}=-1.4$, $\alpha=0.8$, $\beta=0.9$
and the following initial conditions $x~{}(N=0)=0.27$, $z~{}(N=0)=0.08$ from
the region II of Fig. 6. For the values of initial conditions on $x(N)$ and
$y(N)$ from the region I, in a similar fashion we shall obtain $\Omega_{c}=0$
and $\Omega_{d}=1$ at late time. Again if we take initial conditions on $x(N)$
and $y(N)$ from the regions III and IV of Fig. 6, we shall reach
$\Omega_{c}=\frac{\beta}{\alpha+\beta}$ and
$\Omega_{d}=\frac{\alpha}{\alpha+\beta}$ in an asymptotic fashion.  Figure 8:
Phase portrait plot for Model II (eqn. (11)) with $\alpha>0,\beta>0$ and
$w_{d}\leq-1-\frac{\beta}{\alpha}$. In this case we have taken $w_{d}=-2.4$,
$\alpha=0.6$ and $\beta=0.8$. We note that one can take any value of
$w_{d}\leq-1-\frac{\beta}{\alpha}$ and any positive value of $\alpha$ and
$\beta$ in order to get similar graphics. Here, the yellow shaded region
represents the accelerated region (i.e. $q<0$) and the pink shaded region
corresponds to the decelerated region (i.e. $q>0$).

####  III.2.2 Dynamical $w_{d}$

Now, consider the dynamical $w_{d}$ as in eqn. (29), for which the autonomous
system (41) takes the form:

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\left(\frac{z}{1-z}\right)\left(\alpha
x-\beta(1-x)\right)-3x(1-x)\left(1+\nu\frac{(1-z)^{2}(1-x)}{z^{2}}\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z\left(x-\nu\frac{(1-z)^{2}(1-x)^{2}}{z^{2}}\right),\end{array}\right.$
(47)

where $\nu=\frac{3AH_{0}^{2}}{\kappa^{2}}$. We regularize the autonomous
system (47) by multiplying the factor $z^{2}(1-z)$ on the right hand sides of
(47) and finally obtain the following autonomous system which is topologically
equivalent to (47).

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-{z^{3}}\left(\alpha
x-\beta(1-x)\right)-3(1-z)x(1-x)\left(z^{2}+\nu(1-z)^{2}(1-x)\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z\left(xz^{2}-\nu(1-z)^{2}(1-x)^{2}\right).\end{array}\right.$
(50)

Now, the critical points of the system (50) are

  * •

$\bar{B}_{0}=\left(\frac{\beta}{\alpha+\beta},1\right)$, $\bar{B}_{1}=(0,0)$,
$\bar{B}_{2}=(1,0)$,
$S=\left\\{\left(x_{c},\frac{3x_{c}}{3x_{c}-(\alpha+\beta)x+\beta}\right)\right\\},$

where $S$ represents the set of critical points in which $x_{c}$ denotes a
real root of $h(x)\equiv
9x^{3}-\nu((\alpha+\beta)x-\beta)^{2}(1-x)^{2}=0$.666Note that $h(x)$ can be
obtained from the following two nullclines:  $\displaystyle
xz^{2}-\nu(1-z)^{2}(1-x)^{2}=0,$ (51) $\displaystyle-z^{3}[\alpha
x-\beta(1-x)]-3x(1-x)(1-z)\left[z^{2}+\nu(1-z)^{2}(1-x)\right]=0.$ (52) As
$h(x)$ is a fourth degree equation in $x$, therefore, the set $S$ may contain
maximum four critical points. Now, since our physical domain is $R$,
therefore, we are interested to investigate the number of roots of $h(x)$
within the interval $[0,1]$. As $h(0)=-{\beta}^{2}\nu<0$ (for $\nu>0$) and
$h(1)=9>0$, hence, by the Bolzano’s theorem [169], $h(x)$ will have at least
one real root in $(0,1)$. Note that the $z$ component of the critical point
should satisfy $0\leq z\leq 1$. Now, for any $x_{c}$ in $(0,1)$, the condition
on the $z$ component, $\frac{3x_{c}}{3x_{c}-(\alpha+\beta)x+\beta}\leq 1$
leads to $x_{c}\leq\beta/(\alpha+\beta)$. Consequently, we need to check the
number of roots of $h(x)$ in $[0,\frac{\beta}{\alpha+\beta}]$. Again, we
notice that $h(0)<0$ and
$h\left(\beta/(\alpha+\beta)\right)=9\left((\beta/(\alpha+\beta)\right)^{3}>0$
(since $\alpha>0$, $\beta>0$). Thus, from the Bolzano’s theorem [169], we
claim that there is at least one root of $h(x)$ in
$\left(0,\beta/(\alpha+\beta)\right)$. Now, looking at the derivative of
$h(x)$ given by

$\displaystyle h^{\prime}(x)$
$\displaystyle=-4(\alpha+\beta)^{2}\nu(x-1)\left(x-\frac{\beta}{\alpha+\beta}\right)\left(x-\frac{\alpha+2\beta}{2\alpha+2\beta}\right)+27x^{2},$

we can see that $h^{\prime}(x)$ is positive in
$\left(0,\beta/(\alpha+\beta)\right)$ i.e., $h(x)$ is strictly increasing in
$\left(0,\beta/(\alpha+\beta)\right)$. Hence, $h(x)$ has only one root in
$\left(0,\beta/(\alpha+\beta)\right)$ and correspondingly, the set $S$
contains only one critical point and we label this critical point as
$\bar{B}_{3}$. In this case, we also observe that the qualitative nature of
the critical point $\bar{B}_{3}$ is same as $B_{3}$ which has been described
earlier and therefore, the phase portrait will be same as Fig. 6.

On the other hand, for $\nu<0$, the algebraic curve represented by the eqn.
(51) has no branches in the positive quadrant. Hence, for $\nu<0$, $S$ does
not have any critical point, that means, the autonomous system in this case
has only three critical points, namely, $\bar{B}_{0}$, $\bar{B}_{1}$,
$\bar{B}_{2}$ and the phase plot will be similar to Fig. 5.

###  III.3 Model III

In this section we describe the dynamical analysis for the interacting
scenario driven by the interaction function $Q_{\rm III}$ of (12). Using the
same dynamical variables ($x$, $z$) defined as

$\displaystyle x=\frac{\kappa^{2}\rho_{c}}{3H^{2}},\quad\quad
z=\frac{H_{0}}{H+H_{0}},$ (53)

the autonomous system for this interacting scenario takes the form

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma\left(\frac{z}{1-z}\right)\left(x^{2}+x-1\right)+3w_{d}x(1-x),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z\left(1+w_{d}(1-x)\right),\end{array}\right.$
(56)

where $\gamma=\Gamma/H_{0}$ is the dimensionless parameter. Now, regularizing
the vector fields, as we have described in III.1, the autonomous system (56)
can be reduced to the form

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma
z\left(x^{2}+x-1\right)+3w_{d}(1-z)x(1-x),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z\left(1+w_{d}(1-x)\right).\end{array}\right.$
(59)

Here, $R=[0,1]^{2}$ is the physical domain and proceeding as in earlier we see
that $R$ will be positively invariant if we restrict the parameter $\gamma$ by
$\gamma>0$ (i.e. $\Gamma>0$). We, now, investigate the autonomous system (59)
in terms of the nature of the critical points and their implications for both
constant and dynamical $w_{d}$.

####  III.3.1 Constant $w_{d}$

For constant $w_{d}$, the critical points of the autonomous system (59), their
existence, stability and as well as the cosmological parameters evaluated at
those critical points are summarized in Table 3. In the following we
investigate the nature of the critical points for three different regions of
$w_{d}$, namely, quintessence (i.e. $w_{d}>-1$); cosmological constant (i.e.
$w_{d}=-1$); phantom (i.e. $w_{d}<-1$).

  1. 1.

When $w_{d}>-1$, the point $C_{3}$ does not belong to the physical domain and
$w_{\rm tot}=w_{d}(1-x)>-1$, which means
$z^{\prime}=\frac{3}{2}z(1-z)^{2}(1+w_{\rm tot})>0$. Now, on $z=1$,
$x^{\prime}$ is positive for $x<\frac{\sqrt{5}-1}{2}$ and $x^{\prime}$ is
negative for $x>\frac{\sqrt{5}-1}{2}$. So, $C_{1}$, $C_{2}$ are unstable and
$C_{0}$ is a global attractor. Note that $C_{0}$ corresponds to $H=0$,
$\Omega_{c}=\frac{\sqrt{5}-1}{2}$ and $\Omega_{d}=\frac{3-\sqrt{5}}{2}$. The
phase plot is displayed in Fig. 9.

  2. 2.

When $w_{d}=-1$, we have $C_{1}=C_{3}$, and
$z^{\prime}=\frac{3}{2}z(1-z)^{2}x$ which is positive. On $z=1$ line, one has
$x^{\prime}>0$ for $x<\frac{\sqrt{5}-1}{2}$ and $x^{\prime}<0$ for
$x>\frac{\sqrt{5}-1}{2}$. Once again, $C_{1}$, $C_{2}$ are unstable critical
points and the point $C_{0}$ continues being a global attractor. Again, the
phase plot is shown in Fig. 9.

  3. 3.

For $w_{d}<-1$, the parameter space can be categorized in two ways:

     * •

When $-\frac{\sqrt{5}+3}{2}<w_{d}<-1$, the point $C_{3}$ enters in the
physical domain and we obtain $0<\frac{1+w_{d}}{w_{d}}<\frac{\sqrt{5}-1}{2}$.
Now, if $x<\frac{1+w_{d}}{w_{d}}$, implies $1+w_{\rm tot}<0$ which gives
$z^{\prime}=\frac{3}{2}z(1-z)^{2}(1+w_{\rm tot})<0$. Similarly, $z^{\prime}>0$
for $x>\frac{1+w_{d}}{w_{d}}$. Again, on $z=1$, we obtain $x^{\prime}>0$ for
$x<\frac{\sqrt{5}-1}{2}$ and $x^{\prime}<0$ for $x>\frac{\sqrt{5}-1}{2}$.
Also, $x^{\prime}$ is negative on $z=0$. Thus, the physical region $R$ is
divided into four regions. Trajectories from regions I and II converge to
$C_{1}$ and trajectories from regions III and IV converge to $C_{0}$. The Fig.
10 shows the qualitative nature and Fig. 11 displays the evolution of
$\Omega_{c}$, $\Omega_{d}$ and $w_{\rm tot}$.

     * •

When $w_{d}\leq-\frac{\sqrt{5}+3}{2}$: For the special case with
$w_{d}=-\frac{\sqrt{5}+3}{2}$, one has $C_{0}=C_{3}$ and
$\frac{1+w_{d}}{w_{d}}=\frac{\sqrt{5}-1}{2}$. Now, $z^{\prime}$ is positive
for $x>\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}$ is negative for
$x<\frac{1+w_{d}}{w_{d}}$. On $z=0$, $x^{\prime}$ is negative which shows
$C_{1}$ is a global attractor. Again for $w_{d}<-\frac{\sqrt{5}+3}{2}$, the
point $C_{3}$ leaves the physical domain and one obtains
$\frac{\sqrt{5}-1}{2}<\frac{1+w_{d}}{w_{d}}<1$. Also, one has $z^{\prime}>0$
for $x>\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}<0$ for
$x<\frac{1+w_{d}}{w_{d}}$. Thus, $C_{1}$ is a global attractor. In this case,
the qualitative behavior is given in the Fig. 12.

Point | $x$ | $z$ | Existence | Stability | Acceleration | $\Omega_{c}$ | $\Omega_{d}$ | $w_{\rm tot}$
---|---|---|---|---|---|---|---|---
$C_{0}$ | $\frac{\sqrt{5}-1}{2}$ | 1 |  $\gamma~{}(>0)$ and $w_{d}<-\frac{1}{3}$ | $-\frac{3+\sqrt{5}}{2}<w_{d}<-\frac{1}{3}$ | $w_{d}<-\frac{3+\sqrt{5}}{6}$ | $\frac{\sqrt{5}-1}{2}$ | $\frac{3-\sqrt{5}}{2}$ | $\frac{w_{d}(3-\sqrt{5})}{2}$
$C_{1}$ | $0$ | $0$ |  $\gamma~{}(>0)$ and $w_{d}<-\frac{1}{3}$ | $w_{d}<-1$ | $w_{d}<-\frac{1}{3}$ | $0$ | $1$ | $w_{d}$
$C_{2}$ | $1$ | $0$ |  $\gamma~{}(>0)$ and $w_{d}<-\frac{1}{3}$ | unstable | no | $1$ | $0$ | $0$
$C_{3}$ | $\frac{1+w_{d}}{w_{d}}$ | $\frac{3w_{d}(1+w_{d})}{(3-\gamma){w_{d}}^{2}+3w_{d}(1-\gamma)-\gamma}$ |  $\gamma~{}(>0)$ and $-\frac{3+\sqrt{5}}{2}\leq w_{d}\leq-1$ | unstable | yes | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{1}{w_{d}}$ | $-1$
Table 3: The critical points, their existence, stability, and the values of
the cosmological parameters evaluated at those points for the interacting
scenario driven by the interaction function $Q_{\rm
III}=\Gamma\left(\rho_{c}-\rho_{d}-\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}\right)$
of eqn. (12) are summarized.  Figure 9: Phase plot for Model III (eqn. (12))
with $\gamma>0$ and $w_{d}\geq-1$. In this case we have taken $w_{d}=-0.9$ and
$\gamma=0.5$. We note that one can take any value of $w_{d}\geq-1$ and any
positive value of $\gamma$ in order to get similar graphics. Here, the yellow
shaded region represents the accelerated region (i.e. $q<0$) and the pink
shaded region corresponds to the decelerated region (i.e. $q>0$).  Figure 10:
Phase plot for Model III (eqn. (12)) with $\gamma>0$ and
$-\frac{3+\sqrt{5}}{2}<w_{d}<-1$. In this case we have taken $w_{d}=-1.35$ and
$\gamma=0.8$. We note that one can take any specific value of $\gamma~{}(>0)$
to draw the plot, however, as long as $\gamma$ decreases, the regions I and IV
become very small and they look indistinguishable from one another.  Figure
11: We display the evolution of the CDM density parameter $(\Omega_{c})$, dark
energy density parameter $(\Omega_{d})$ and the total equation of state
parameter $(w_{\rm tot})$ for the Model III (eqn. (12)). We have chosen the
following values of the parameters: $w_{d}=-1.35$, $\gamma=0.8$ and the
following initial conditions: $x~{}(N=0)=0.25$, $z~{}(N=0)=0.07$ from the
region II of Fig. 10. In addition, if we start any trajectory from region I of
Fig. 10, it converges to the critical point $C_{1}$. So, for any initial
conditions from region I, we shall get $\Omega_{c}=0$ and $\Omega_{d}=1$ at
late time. Any trajectories starting from regions III and IV of Fig. 10 will
converge to the critical point $C_{0}$. Therefore, if we take initial
conditions on $x(N)$ and $z(N)$ from regions III and IV, we shall reach
$\Omega_{c}=\frac{\sqrt{5}-1}{2}$ and $\Omega_{d}=\frac{3-\sqrt{5}}{2}$ in an
asymptotic fashion.  Figure 12: Phase plot for Model III (eqn. (12)) with
$\gamma>0$ and $w_{d}\leq-\frac{3+\sqrt{5}}{2}$. In this case we have taken
$w_{d}=-2.9$ and $\gamma=0.7$. We note that one can take any value of $w_{d}$
satisfying $w_{d}\leq-\frac{3+\sqrt{5}}{2}$ and any positive value of $\gamma$
in order to get similar graphics. Here, the yellow shaded region represents
the accelerated region (i.e. $q<0$) and the pink shaded region corresponds to
the decelerated region (i.e. $q>0$).

####  III.3.2 Dynamical $w_{d}$

We consider the dynamical $w_{d}$ of the form (29), for which the autonomuous
system (56) becomes:

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma\left(\frac{z}{1-z}\right)\left(x^{2}+x-1\right)-3x(1-x)\left(1+\nu\frac{(1-z)^{2}(1-x)}{z^{2}}\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z\left(x-\nu\frac{(1-z)^{2}(1-x)^{2}}{z^{2}}\right),\end{array}\right.$
(62)

where $\nu=\frac{3AH_{0}^{2}}{\kappa^{2}}$. The topologically equivalent
autonomous system after regularization is given by the following:

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\gamma{z^{3}}\left(x^{2}+x-1\right)-3(1-z)x(1-x)\left(z^{2}+\nu(1-z)^{2}(1-x)\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z\left(xz^{2}-\nu(1-z)^{2}(1-x)^{2}\right).\end{array}\right.$
(65)

The critical points of the system (65) are

  * •

$\bar{C}_{0}=\left(\frac{\sqrt{5}-1}{2},1\right)$, $\bar{C}_{1}=(0,0)$,
$\bar{C}_{2}=(1,0)$,
$S=\left\\{\left(x_{c},\frac{3x_{c}}{3x_{c}-\gamma\left({x_{c}}^{2}+x_{c}-1\right)}\right)\right\\},$

where $S$ is the set of critical points in which $x_{c}$ is a real root of
$g(x)\equiv
9x^{3}-{\gamma}^{2}{\nu}\left(x^{2}+x-1\right)^{2}(1-x)^{2}=0$.777Here $g(x)$
can be obtained from the following two nullclines:  $\displaystyle
xz^{2}-\nu(1-z)^{2}(1-x)^{2}=0,$ (66) $\displaystyle-\gamma
z^{3}\left(x^{2}+x-1\right)-3x(1-x)(1-z)\left[z^{2}+\nu(1-z)^{2}(1-x)\right]=0.$
(67) Note that as $g(x)$ is a sixth degree equation, therefore, the set $S$
may have maximum six critical points. Thus, it is now important to calculate
the number of critical points in order to understand the phase space of the
interacting scenario for this dynamical $w_{d}$. As $g(0)=-{\gamma}^{2}\nu<0$
(for $\nu>0$) and $g(1)=9>0$, it follows by the Bolzano’s theorem [169] that
$g(x)$ has at least one real root in $(0,1)$. Now, for any $x_{c}$ lying in
$(0,1)$, the condition
$z_{c}=\frac{3x_{c}}{3x_{c}-\gamma\left({x_{c}}^{2}+x_{c}-1\right)}\leq 1$
requires $x_{c}\leq\frac{\sqrt{5}-1}{2}$. We also observe that
$g\left(\frac{\sqrt{5}-1}{2}\right)=9(-2+\sqrt{5})>0$, and hence, from the
Bolzano’s theorem [169], there is at least one root of $g(x)$ in
$(0,\frac{\sqrt{5}-1}{2})$. We now claim that $g(x)$ has only one root in the
interval $[0,\frac{\sqrt{5}-1}{2}]$, and in that case, it is straightforward
to conclude that there is only one critical point in the set $S$. Now, from
the derivative of $g(x)$ given below

$\displaystyle g^{\prime}(x)$
$\displaystyle=-6{\gamma}^{2}\nu\left(x-1\right)\left(x-\sqrt{\frac{2}{3}}~{}\right)\left(x+\sqrt{\frac{2}{3}}~{}\right)\left(x-\frac{\sqrt{5}-1}{2}\right)\left(x+\frac{\sqrt{5}+1}{2}\right)+27x^{2},$

we notice that $g^{\prime}(x)>0$ for $x\in(0,\frac{\sqrt{5}-1}{2})$ which
means that $g(x)$ is strictly increasing in $(0,\frac{\sqrt{5}-1}{2})$.
Therefore, the set $S$ contains only one critical point and we label this
critical point as $\bar{C}_{3}$. In this case, the topological nature of the
critical point $\bar{C}_{3}$ is same as $C_{3}$ described earlier and hence,
the phase plot is same as the Fig. 10.

On the other hand, looking at eqn. (66), and following the similar arguments
as in earlier for $\nu<0$, we claim that the set $S$ does not have any
critical point for $\nu<0$. Thus, in this case, we have only three critical
points, namely, $\bar{C}_{0}$, $\bar{C}_{1}$, and $\bar{C}_{2}$, and the phase
portrait looks like same as the Fig. 9.

###  III.4 Model IV

In this section we describe the dynamical analysis for the interaction
function $Q_{\rm IV}$ of eqn. (13). Considering the $(x,z)$ variables defined
as

$\displaystyle x=\frac{\kappa^{2}\rho_{c}}{3H^{2}},\quad\quad
z=\frac{H_{0}}{H+H_{0}},$ (68)

the autonomous system for this interaction model takes the form

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\left(\frac{z}{1-z}\right)\left(\alpha
x-\beta x(1-x)\right)+3w_{d}x(1-x),\\\
z^{\prime}&=&\frac{3}{2}z(1-z)\left(1+w_{d}(1-x)\right),\end{array}\right.$
(71)

where $\alpha$, $\beta$ are defined as $\alpha=\Gamma_{c}/H_{0}$ and
$\beta=\Gamma_{cd}/H_{0}$ respectively and they are dimensionless. Note again
that $\alpha\neq\beta$ for the interaction model to be sign shifting. Here, we
shall regularize the autonomous system (71) similar to the procedure we have
applied in (III.1). The regularized autonomous system will be of the form

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-z\left(\alpha x-\beta
x(1-x)\right)+3w_{d}(1-z)x(1-x),\\\
z^{\prime}&=&\frac{3}{2}z(1-z)^{2}\left(1+w_{d}(1-x)\right).\end{array}\right.$
(74)

Also we note that $R=[0,1]^{2}$ is the physical domain and the lines $x=0$,
$z=0$, $z=1$ are invariant under the flow generated by the autonomous system
(74). We also see that along the line $x=1$, one obtains $x^{\prime}=-\alpha
z$. Hence, the physical domain $R$ will be positively invariant if we restrict
the parameter $\alpha$ by $\alpha>0$ (i.e., $\Gamma_{c}>0$). It has been
mentioned in (II.1) that the interaction function $Q_{\rm IV}$ will exhibit
sign shifting property provided $\Gamma_{c}$ and $\Gamma_{cd}$ are of the same
sign. Consequently, to be consistent with the positively invariance of the
physical domain $R$ and the sign shifting nature of the interaction function
$Q_{\rm IV}$, we will restrict our attention here on $\alpha>0,\beta>0$.

In the following we investigate the nature of the critical points and their
implications for the constant and non-constant nature of $w_{d}$.

####  III.4.1 Constant $w_{d}$

In this section we explore the nature of the critical points for different
regions of $w_{d}$. In Table 4, we summarize the critical points of the
autonomous system (74), their existence, stability and as well as the
cosmological parameters evaluated at those critical points. In the following
we investigate the nature of the critical points for three different regions
of $w_{d}$, namely, quintessence (i.e. $w_{d}>-1$); cosmological constant
(i.e. $w_{d}=-1$); phantom (i.e. $w_{d}<-1$).

  * •

When $\beta>\alpha>0$, the following situations arise:

    1. 1.

$w_{d}>-1$: When $w_{d}>-1$, the point $D_{3}$ does not belong to the physical
region $R$ and also $w_{\rm tot}=w_{d}(1-x)>-1$ which implies
$z^{\prime}=\frac{3}{2}z(1-z)^{2}(1+w_{\rm tot})>0$. In addition as
$\beta>\alpha>0$, the point $D_{0}$ enters in the region $R$. On $z=1$ line,
the value of $x^{\prime}$ is positive if $x<\frac{\beta-\alpha}{\beta}$ and
$x^{\prime}$ is negative if $x>\frac{\beta-\alpha}{\beta}$. So, $D_{00}$,
$D_{1}$ and $D_{2}$ are unstable critical points and $D_{0}$ is a global
attractor. At the point $D_{0}$, we have $H=0$,
$\Omega_{c}=\frac{\beta-\alpha}{\beta}$ and $\Omega_{d}=\frac{\alpha}{\beta}$.
Left plot of Fig. 13 shows the behavior.

    2. 2.

$w_{d}=-1$: When $w_{d}=-1$, all points in the $OZ$ axis are critical points
and also point $D_{3}$ lies on $OZ$ axis. Again, one gets
$z^{\prime}=\frac{3}{2}z(1-z)^{2}x$ which is positive. At $z=1$, one has
$x^{\prime}>0$ for $x<\frac{\beta-\alpha}{\beta}$ and $x^{\prime}<0$ for
$x>\frac{\beta-\alpha}{\beta}$. In this case, $D_{0}$ is an attractor but it
is not a global attractor. The qualitative nature is displayed in right plot
of Fig. 13.

    3. 3.

When $-\frac{\beta}{\alpha}<w_{d}<-1$: The critical point $D_{3}$ enters in
the physical domain. We have
$0<\frac{1+w_{d}}{w_{d}}<\frac{\beta-\alpha}{\beta}$. At $z=0$,
$x^{\prime}=3w_{d}x(1-x)$ which is negative and at $z=1$, we obtain
$x^{\prime}>0$ for $x<\frac{\beta-\alpha}{\beta}$ and $x^{\prime}<0$ for
$x>\frac{\beta-\alpha}{\beta}$. If $x<\frac{1+w_{d}}{w_{d}}$, implies
$1+w_{\rm tot}<0$. So, $z^{\prime}$ is negative. Similarly, $z^{\prime}$ is
positive for $x>\frac{1+w_{d}}{w_{d}}$. Thus, all orbits in the regions I and
II of left plot of Fig. 14, at late time, converge to $D_{1}$. For an orbit in
the regions III and IV of left plot Fig. 14, at late time, it converges to
$D_{0}$. Note that $D_{0}$ means $H=0$ with
$\Omega_{c}=\frac{\beta-\alpha}{\beta},\Omega_{d}=\frac{\alpha}{\beta}$ and
$w_{\rm tot}=\frac{\alpha w_{d}}{\beta}>-1$. On the contrary, $D_{1}$ means
$H=\infty$ with $\Omega_{d}=1$ and $w_{\rm tot}=w_{d}<-1$. The Fig. 15
presents the evolution of $\Omega_{c}$, $\Omega_{d}$, and $w_{\rm tot}$.

    4. 4.

$w_{d}\leq-\frac{\beta}{\alpha}$: When $w_{d}=-\frac{\beta}{\alpha}$, one
obtains $D_{0}=D_{3}$. Now, one has $z^{\prime}<0$ for
$x<\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}>0$ for $x>\frac{1+w_{d}}{w_{d}}$.
Also, $x^{\prime}$ is negative on $z=0$. Thus, $D_{1}$ is a global attractor.
For $w_{d}<-\frac{\beta}{\alpha}$, the critical point $D_{3}$ leaves the
physical region and one gets
$\frac{\beta-\alpha}{\beta}<\frac{1+w_{d}}{w_{d}}<1$ which gives
$z^{\prime}<0$ for $x<\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}>0$ for
$x>\frac{1+w_{d}}{w_{d}}$. Here also, we have $x^{\prime}<0$ on $z=0$. Once
again, $D_{1}$ is a global attractor. Right plot of Fig. 14 exhibits the phase
plot.

  * •

$\alpha>\beta>0$: In this case the critical points $D_{0}$ and $D_{3}$ never
belong to the physical domain.

    * –

When $w_{d}>-1$, $z^{\prime}$ is positive as $1+w_{\rm tot}>0$. At $z=1$ line,
$x^{\prime}$ is negative. Thus, $D_{00}$ is a global attractor. The point
$D_{00}$ corresponds to $H=0$, $\Omega_{c}=0$ and $\Omega_{d}=1$.

    * –

When $w_{d}=-1$, all points in the $OZ$ axis are critical points. Now, one has
$z^{\prime}=\frac{3}{2}z(1-z)^{2}x$ which is positive. Also, $x^{\prime}$ is
negative on both $z=0$ and $z=1$ lines. Here, we have no attractor.

    * –

For $w_{d}<-1$, we have $0<\frac{1+w_{d}}{w_{d}}<1$. On $z=0$ line,
$x^{\prime}=3w_{d}x(1-x)$ which is negative. Here, we have $z^{\prime}<0$ for
$x<\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}>0$ for $x>\frac{1+w_{d}}{w_{d}}$.
Thus, $D_{1}$ is the only global attractor. In all the cases, qualitative
behaviors are shown in Fig. 16.

Point | $x$ | $z$ | Existence | Stability | Acceleration | $\Omega_{c}$ | $\Omega_{d}$ | $w_{\rm tot}$
---|---|---|---|---|---|---|---|---
$D_{0}$ | $\frac{\beta-\alpha}{\beta}$ | 1 |  $\beta>\alpha>0$ and $w_{d}<-\frac{1}{3}$ | $-\frac{\beta}{\alpha}<w_{d}<-\frac{1}{3}$ | $w_{d}<-\frac{\beta}{3\alpha}$ | $\frac{\beta-\alpha}{\beta}$ | $\frac{\alpha}{\beta}$ | $\frac{\alpha w_{d}}{\beta}$
$D_{00}$ | $0$ | $1$ |  $\alpha>0$, $\beta>0$ and $w_{d}<-\frac{1}{3}$ |  $\alpha~{}(>0)>\beta$ with $-1<w_{d}<-\frac{1}{3}$ | $w_{d}<-\frac{1}{3}$ | $0$ | $1$ | $w_{d}$
$D_{1}$ | $0$ | $0$ |  $\alpha>0$, $\beta>0$ and $w_{d}<-\frac{1}{3}$ | $w_{d}<-1$ | $w_{d}<-\frac{1}{3}$ | $0$ | $1$ | $w_{d}$
$D_{2}$ | $1$ | $0$ |  $\alpha>0$, $\beta>0$ and $w_{d}<-\frac{1}{3}$ | unstable | no | $1$ | $0$ | $0$
$D_{3}$ | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{3w_{d}}{(\alpha-3)w_{d}+\beta}$ |  $\beta>\alpha>0$ with $-\frac{\beta}{\alpha}\leq w_{d}\leq-1$ | unstable for $w_{d}\neq-1$ | yes | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{1}{w_{d}}$ | $-1$
Table 4: The critical points, their existence, stability, and the values of
the cosmological parameters evaluated at those points for the interacting
scenario driven by the interaction function $Q_{\rm
IV}=\Gamma_{c}\rho_{c}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}$
of eqn. (13) are summarized.

Figure 13: Left plot: Phase plot for Model IV (eqn. (13)) with $w_{d}>-1$ and
$\beta>\alpha>0$. In this case we have used $w_{d}=-0.92$, $\alpha=0.3$ and
$\beta=0.8$. Right plot: Phase plot for Model IV (eqn. (13)) with $w_{d}=-1$
and $\beta>\alpha>0$. Particularly we have taken $w_{d}=-1$, $\alpha=0.31$ and
$\beta=0.9$. We note that the yellow shaded region represents the accelerated
region (i.e. $q<0$) and the pink shaded region corresponds to the decelerated
region (i.e. $q>0$).

Figure 14: Left Plot: Phase plot for Model IV (eqn. (13)) with
$-\frac{\beta}{\alpha}<w_{d}<-1$ and $\beta>\alpha>0$. For, numerical
simulation we have used $w_{d}=-2.2$, $\alpha=1$ and $\beta=6$. Right Plot:
Phase plot for Model IV (eqn. (13)) with $w_{d}\leq-\frac{\beta}{\alpha}$ and
$\beta>\alpha>0$. For numerical simulation, we have used $w_{d}=-4.1$,
$\alpha=0.2$ and $\beta=0.8$ from the parameter space. In this plot, the
yellow shaded region corresponds to the accelerated region (i.e. $q<0$) and
the pink shaded region corresponds to the decelerated region (i.e. $q>0$).
Figure 15: We display the evolution of the CDM density parameter
$(\Omega_{c})$, dark energy density parameter $(\Omega_{d})$ and the total
equation of state parameter $(w_{\rm tot})$ for the Model IV (eqn. (13)). We
have chosen the following values of the parameters: $w_{d}=-2.2$, $\alpha=1$,
$\beta=6$ and the following initial conditions: $x~{}(N=0)=0.24$,
$z~{}(N=0)=0.01$ from the region II of the left plot of Fig. 14. In addition,
if we start any trajectory from region I of the left plot of Fig. 14, it
converges to the critical point $D_{1}$. So, for any initial conditions from
region I, we shall get $\Omega_{c}=0$ and $\Omega_{d}=1$ at late time. Any
trajectory starting from regions III and IV of the left plot of Fig. 14 will
converge to the critical point $D_{0}$. Therefore, if we take initial
conditions on $x(N)$ and $z(N)$ from regions III and IV, we shall reach
$\Omega_{c}=\frac{\beta-\alpha}{\beta}$ and $\Omega_{d}=\frac{\alpha}{\beta}$
in an asymptotic fashion.

Figure 16: Left plot: Phase plot for Model IV (eqn. (13)) with $w_{d}>-1$ and
$\alpha>\beta>0$. In this case we have used $w_{d}=-0.7$, $\alpha=0.2$ and
$\beta=0.1$. Middle plot: Phase plot for Model IV (eqn. (13)) with $w_{d}=-1$
and $\alpha>\beta>0$. Particularly, we have chosen $w_{d}=-1$, $\alpha=0.08$
and $\beta=0.05$. Right plot: Phase plot for Model IV (eqn. (13)) with
$w_{d}<-1$ and $\alpha>\beta>0$. For numerical simulation we have used
$w_{d}=-3.3$, $\alpha=0.2$ and $\beta=0.1$. Here, the yellow shaded region
represents the accelerated region (i.e. $q<0$) and the pink shaded region
corresponds to the decelerated region (i.e. $q>0$).

####  III.4.2 Dynamical $w_{d}$

For the equation of state (29), the autonomuous system (71) becomes:

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\left(\frac{z}{1-z}\right)x\left(\alpha-\beta(1-x)\right)-3x(1-x)\left(1+\nu\frac{(1-z)^{2}(1-x)}{z^{2}}\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z\left(x-\nu\frac{(1-z)^{2}(1-x)^{2}}{z^{2}}\right),\end{array}\right.$
(77)

where $\nu=\frac{3AH_{0}^{2}}{\kappa^{2}}$. Proceeding as in the earlier
cosmological models for non-constant $w_{d}$, we regularize the autonomous
system (77) and obtain the following:

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-{z^{3}}x\left(\alpha-\beta(1-x)\right)-3(1-z)x(1-x)\left(z^{2}+\nu(1-z)^{2}(1-x)\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z\left(xz^{2}-\nu(1-z)^{2}(1-x)^{2}\right),\end{array}\right.$
(80)

The critical points of the system (80) are

  * •

$\bar{D}_{0}=\left(\frac{\beta-\alpha}{\beta},1\right)$, $\bar{D}_{00}=(0,1)$,
$\bar{D}_{1}=(0,0)$, $\bar{D}_{2}=(1,0)$,
$S=\left\\{\left(x_{c},\frac{3}{3-\alpha+\beta(1-x_{c})}\right)\right\\},$

where $x_{c}$ is a root of $\Theta(x)\equiv
9x-\nu(-\alpha+\beta(1-x))^{2}(1-x)^{2}=0$.888Note again that $\Theta(x)$ is
obtained from the following two nullclines:  $\displaystyle
xz^{2}-\nu(1-z)^{2}(1-x)^{2}=0,$ (81)
$\displaystyle-z^{3}x[\alpha-\beta(1-x)]-3x(1-x)(1-z)\left[z^{2}+\nu(1-z)^{2}(1-x)\right]=0.$
(82) Since $\Theta(x)$ represents a fourth degree equation in $x$, therefore,
$S$ may contain maximum four critical points. In a similar fashion, we proceed
to investigate the number of roots of $\Theta(x)$ in $[0,1]$. Now, since
$\Theta(0)=\nu(\beta-\alpha)^{2}<0$ (for $\nu>0$) and $\Theta(1)=9>0$, from
the Bolzano’s theorem [169], we can say that $\Theta(x)$ has at least one root
in the interval $(0,1)$. However, as the physical domain in this case is $R$,
therefore, for valid critical point in this domain, we have to check whether
the $z$ component of the critical point lies in $[0,1]$. Now, for any $x_{c}$
in $[0,1]$, the condition $\frac{3}{3-\alpha+\beta(1-x_{c})}\leq 1$ demands
that $\beta x_{c}\leq\beta-\alpha$. Now, as $\alpha>0$, $\beta>0$, we consider
the following cases:

  * •

$\beta>\alpha>0$: In this case, the condition $\beta x_{c}\leq\beta-\alpha$
reduces to $x_{c}\leq\frac{\beta-\alpha}{\beta}$. We see that $\Theta(0)<0$
and
$\Theta\left(\frac{\beta-\alpha}{\beta}\right)=9\frac{\beta-\alpha}{\beta}>0$,
hence, from the Bolzano’s theorem [169], there is at least one root of
$\Theta(x)$ in $\left(0,\frac{\beta-\alpha}{\beta}\right)$. We also observe
that

$\displaystyle\Theta^{\prime}(x)$
$\displaystyle=-4{\beta}^{2}\nu\left(x-1\right)\left(x-\frac{\beta-\alpha}{\beta}\right)\left(x-\frac{2\beta-\alpha}{2\beta}\right)+9>0,$

for $x\in\left(0,\frac{\beta-\alpha}{\beta}\right)$. Hence, $\Theta(x)$ is
strictly increasing in $x\in\left(0,\frac{\beta-\alpha}{\beta}\right)$ and as
a result $\Theta(x)$ has only one root in
$x\in\left(0,\frac{\beta-\alpha}{\beta}\right)$. Thus, the set $S$ contains
only one critical point which we call as $\bar{D}_{3}$. The critical point
$\bar{D}_{3}$ behaves qualitatively same as the point $D_{3}$ and
consequently, the phase portrait is topologically same as the left plot of
Fig. 14.

Again, following the same arguments as in the earlier models for $\nu<0$, we
can show that $\Theta(x)$ has no root in $(0,1)$ for $\nu<0$. Hence, in this
case the autonomous system has only four critical points, namely,
$\bar{D}_{0}$, $\bar{D}_{00}$, $\bar{D}_{1}$, $\bar{D}_{2}$. The phase space
stability analysis has been shown before. The phase portrait is same as the
left plot of Fig. 13.

  * •

$\alpha>\beta>0$: In this case, we have $\beta-\alpha<0$, and hence,
$\bar{D}_{0}$ does not belong to the physical domain $R$. Also, the condition
$\beta x_{c}\leq\beta-\alpha$, leads to the claim that $x_{c}$ does not belong
to the domain $[0,1]$. Therefore, we have only three critical points,
$\bar{D}_{00}$, $\bar{D}_{1}$, $\bar{D}_{2}$ and the phase portrait is same as
that of the phase portrait of the right plot of Fig. 16. On the other hand,
for $\nu<0$, we do not have any critical point in $S$. Here also we have only
three critical points, namely, $\bar{D}_{00}$, $\bar{D}_{1}$, $\bar{D}_{2}$
and we see that the corresponding phase portrait looks similar to left plot of
Fig. 16.

Figure 17: Left plot: Phase plot for Model V (eqn. (14)) with $w_{d}\geq-1$
and $\beta<\alpha<0$. In this case we have chosen $w_{d}=-0.9$, $\alpha=-0.3$
and $\beta=-0.6$. Middle plot: Phase plot for Model V (eqn. (14)) with
$w_{d}\leq\frac{\beta}{\alpha-\beta}$ and $\beta<\alpha<0$. In this case we
have chosen $w_{d}=-2.05$, $\alpha=-0.3$ and $\beta=-0.6$. Right plot: Phase
plot for Model V (eqn. (14)) with $w_{d}\geq-1$ and $\alpha<\beta<0$. In this
case we have chosen $w_{d}=-0.9$, $\alpha=-0.6$ and $\beta=-0.3$. Here, the
yellow shaded region represents the accelerated region (i.e. $q<0$) and the
pink shaded region corresponds to the decelerated region (i.e. $q>0$).

Figure 18: Left plot: Phase plot for Model V (eqn. (14)) with $w_{d}<-1$ and
$\alpha<\beta<0$. In this case we have chosen $w_{d}=-2$, $\alpha=-2$ and
$\beta=-0.1$. Right plot: Phase plot for Model V (eqn. (14)) with
$\frac{\beta}{\alpha-\beta}<w_{d}<-1$ and $\beta<\alpha<0$. In this case we
have taken $w_{d}=-1.38$, $\alpha=-0.6$ and $\beta=-0.8$. We note that one can
take any specific value of $\alpha~{}(<0)$ and $\beta~{}(<0)$ with
$\beta<\alpha<0$ to draw the plot, however, as long as $\alpha$ and $\beta$
increase, the regions I and IV become very small and they look
indistinguishable from one another.  Figure 19: We display the evolution of
the CDM density parameter $(\Omega_{c})$, dark energy density parameter
$(\Omega_{d})$ and the total equation of state parameter $(w_{\rm tot})$ for
Model V (eqn. (14)). We have chosen the following values of the parameters:
$w_{d}=-1.38$, $\alpha=-0.6$, $\beta=-0.8$ and the following initial
conditions: $x~{}(N=0)=0.3$, $z~{}(N=0)=0.1$ from the region II of the right
plot of Fig. 18. In addition, if we start any trajectory from region I of the
right plot of Fig. 18, it converges to the critical point $E_{1}$. So, for any
initial conditions from region I, we shall get $\Omega_{c}=0$ and
$\Omega_{d}=1$ at late time. Any trajectory starting from regions III and IV
of the right plot of Fig. 18 will converge to the critical point $E_{0}$.
Therefore, if we take initial conditions on $x(N)$ and $z(N)$ from regions III
and IV, we shall reach $\Omega_{c}=\frac{\alpha}{\beta}$ and
$\Omega_{d}=\frac{\beta-\alpha}{\beta}$ in an asymptotic fashion.

###  III.5 Model V

In this section we discuss the dynamical analysis for the interacting scenario
driven by the interaction function $Q_{\rm V}$ of eqn. (14). Using the
dimensionless variables $(x,z)$ defined as

$\displaystyle x=\frac{\kappa^{2}\rho_{c}}{3H^{2}},\quad\quad
z=\frac{H_{0}}{H+H_{0}},$ (83)

we obtain the following autonomous system for the prescribed interacting
scenario:

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\left(\frac{z}{1-z}\right)(1-x)\left(\alpha-\beta
x\right)+3w_{d}x(1-x),\\\
z^{\prime}&=&\frac{3}{2}z(1-z)\left(1+w_{d}(1-x)\right),\end{array}\right.$
(86)

where $\alpha$, $\beta$ are defined as $\alpha=\Gamma_{d}/H_{0}$ and
$\beta=\Gamma_{cd}/H_{0}$ respectively and $\alpha\neq\beta$. After
regularizing, in a similar way we have performed in (III.1), the autonomous
system (86) can be written of the form

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-z(1-x)\left(\alpha-\beta
x\right)+3w_{d}(1-z)x(1-x),\\\
z^{\prime}&=&\frac{3}{2}z(1-z)^{2}\left(1+w_{d}(1-x)\right).\end{array}\right.$
(89)

The physical region is $R$ which is the square $R=[0,1]^{2}$ and following the
similar arguments as in the case of earlier models, we observe that $R$ is
positively invariant if $\alpha<0$ (i.e. $\Gamma_{d}<0$). The interaction
function $Q_{\rm V}$ as mentioned in (14) is of sign shifting nature provided
that the coupling parameters $\Gamma_{d}$ and $\Gamma_{cd}$ are of the same
sign. Hence, for the dynamical analysis driven by the sign shifting
interaction function $Q_{\rm V}$, and to keep the physical domain $R$
positively invariant, we assume the parametric condition $\alpha<0,\beta<0$.

####  III.5.1 Constant $w_{d}$

For constant $w_{d}$, the critical points of the autonomous system (89), their
existence, stability and as well as the cosmological parameters evaluated at
those critical points are summarized in Table 5. Similar to the earlier cases,
here we consider various cases of $w_{d}$ depending on its parameter space. In
what follows we consider various cases.

  * •

We consider the first case where the dimensionless coupling parameter satisfy
the relation $\beta<\alpha<0$:

    1. 1.

When $w_{d}>-1$, $E_{3}$ leaves the domain $R$ and one has $w_{\rm
tot}=w_{d}(1-x)>-1$ which implies $z^{\prime}>0$. On $z=1$ line, one gets
$x^{\prime}>0$ for $x<\frac{\alpha}{\beta}$ and $x^{\prime}<0$ for
$x>\frac{\alpha}{\beta}$. Thus, $E_{0}$ is a global attractor. Left plot of
Fig. 17 shows the behavior.

    2. 2.

For $w_{d}=-1$, one obtains $E_{3}=E_{1}$ and
$z^{\prime}=\frac{3}{2}z(1-z)^{2}x>0$. On $z=1$ line, one has $x^{\prime}>0$
for $x<\frac{\alpha}{\beta}$ and $x^{\prime}<0$ for $x>\frac{\alpha}{\beta}$.
So, $E_{0}$ is a global attractor. Again, the phase plot is given in left plot
of Fig. 17.

    3. 3.

When $\frac{\beta}{\alpha-\beta}<w_{d}<-1$, then $E_{3}$ enters in the
physical region and one can get $\frac{1+w_{d}}{w_{d}}<\frac{\alpha}{\beta}$.
On $z=0$, one obtains $x^{\prime}<0$. Also, on $z=1$, one has $x^{\prime}>0$
for $x<\frac{\alpha}{\beta}$ and $x^{\prime}<0$ for $x>\frac{\alpha}{\beta}$.
Again, $z^{\prime}$ is negative whenever $x<\frac{1+w_{d}}{w_{d}}$ and
$z^{\prime}$ is positive whenever $x>\frac{1+w_{d}}{w_{d}}$. Thus, all orbits
in the regions I and II of the right plot of Fig. 18, at late time, converge
to $E_{1}$. For an orbit in the regions III and IV of the right plot of Fig.
18, at late time, it converges to $E_{0}$. In Fig. 19 we show the evolution of
$\Omega_{c}$, $\Omega_{d}$ and $w_{\rm tot}$.

    4. 4.

When $w_{d}=\frac{\beta}{\alpha-\beta}$, one can obtain $E_{3}=E_{0}$ and
$x^{\prime}$ is negative on $z=0$. Also $z^{\prime}$ is negative for
$x<\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}$ is positive for
$x>\frac{1+w_{d}}{w_{d}}$. Thus, $E_{1}$ is a global attractor. The
qualitative behavior is displayed in middle plot of Fig. 17.

    5. 5.

When $w_{d}<\frac{\beta}{\alpha-\beta}$, $E_{3}$ leaves the physical domain
and one has $\frac{1+w_{d}}{w_{d}}>\frac{\alpha}{\beta}$. As before,
$x^{\prime}$ is negative on $z=0$. Again, $z^{\prime}$ is negative in left
side of $x=\frac{1+w_{d}}{w_{d}}$ and $z^{\prime}$ is positive in right side
of $x=\frac{1+w_{d}}{w_{d}}$. Therefore, $E_{1}$ is a global attractor. The
middle plot of Fig. 17 exhibits the qualitative nature.

  * •

We consider the case when the dimensionless coupling parameters satisfy
$\alpha<\beta<0$. In this case, $E_{0}$ leaves the domain as
$\frac{\alpha}{\beta}>1$. In the following we discuss the nature of the
critical points for different values of $w_{d}$.

    1. 1.

For $w_{d}>-1$, $E_{3}$ leaves the domain $R$. Now since $w_{\rm
tot}=w_{d}(1-x)>-1$, hence, it implies that $z^{\prime}$ is positive. On $z=1$
line, $x^{\prime}$ is positive. As a result, $E_{1}$, $E_{2}$ are unstable and
$E_{00}$ is a global attractor. Right plot of Fig. 17 shows the qualitative
behavior.

    2. 2.

When $w_{d}=-1$, we obtain $E_{1}=E_{3}$ and
$z^{\prime}=\frac{3}{2}z(1-z)^{2}x$ which is positive. At $z=1$ line,
$x^{\prime}$ is positive. Thus, $E_{00}$ is again a global attractor. Phase
plot is displayed in right plot of Fig. 17.

    3. 3.

When $w_{d}<-1$, $E_{3}$ enters in the physical region $R$ and for this
$w_{d}$, we also have $0<\frac{1+w_{d}}{w_{d}}<1$. On $z=0$ line,
$x^{\prime}=3w_{d}x(1-x)$ which is negative and on $z=1$ line, $x^{\prime}$ is
positive. Now left side of $x=\frac{1+w_{d}}{w_{d}}$, $z^{\prime}$ is negative
and right side of $x=\frac{1+w_{d}}{w_{d}}$, $z^{\prime}$ is positive. Hence,
the physical region $R$ is divided into four regions. Thus, all orbits in the
regions I and II of the left plot of Fig. 18, at late time, converge to
$E_{1}$. For an orbit in the regions III and IV of the left plot of Fig. 18,
at late time, it converges to $E_{00}$.

Point | $x$ | $z$ | Existence | Stability | Acceleration | $\Omega_{c}$ | $\Omega_{d}$ | $w_{\rm tot}$
---|---|---|---|---|---|---|---|---
$E_{0}$ | $\frac{\alpha}{\beta}$ | 1 |  $\beta<\alpha<0$ and $w_{d}<-\frac{1}{3}$ | $\frac{\beta}{\alpha-\beta}<w_{d}<-\frac{1}{3}$ | $w_{d}<-\frac{1}{3}\frac{\beta}{\beta-\alpha}$ | $\frac{\alpha}{\beta}$ | $\frac{\beta-\alpha}{\beta}$ | $\frac{(\beta-\alpha)w_{d}}{\beta}$
$E_{00}$ | $1$ | $1$ |  $\alpha<0,\beta<0$ and $w_{d}<-\frac{1}{3}$ |  $\alpha<\beta$ with $w_{d}<-\frac{1}{3}$ | no | $1$ | $0$ | $0$
$E_{1}$ | $0$ | $0$ |  $\alpha<0,\beta<0$ and $w_{d}<-\frac{1}{3}$ | $w_{d}<-1$ | $w_{d}<-\frac{1}{3}$ | $0$ | $1$ | $w_{d}$
$E_{2}$ | $1$ | $0$ |  $\alpha<0,\beta<0$ and $w_{d}<-\frac{1}{3}$ | unstable | no | $1$ | $0$ | $0$
$E_{3}$ | $\frac{1+w_{d}}{w_{d}}$ | $\frac{3w_{d}(1+w_{d})}{3{w_{d}}^{2}+(\alpha-\beta+3)w_{d}-\beta}$ |  $\beta<\alpha<0$ with $\frac{\beta}{\alpha-\beta}\leq w_{d}\leq-1$ |  | yes | $\frac{1+w_{d}}{w_{d}}$ | $-\frac{1}{w_{d}}$ | $-1$
|  |  | and | unstable |  |  |  |
|  |  |  $\alpha<\beta<0$ with $w_{d}\leq-1$ |  |  |  |  |
Table 5: The critical points, their existence, stability, and the values of
the cosmological parameters evaluated at those points for the interacting
scenario driven by the interaction function $Q_{\rm
V}=\Gamma_{d}\rho_{d}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}$
of eqn. (14) are summarized.

Constant $w_{d}$
---
|  | Stability of Critical Points | Deceleration |
Interaction Models | Critical Points | Attractor | Global Attractor | Repeller | Parameter ($q$) | Figures
|  |  |  |  | $q(A_{0})=\frac{1}{2}\left(1+\frac{3w_{d}}{2}\right)<0$ |
|  |  |  $A_{0}$ is a global attractor |  $A_{0}:w_{d}\leq-2$, $\gamma>0$ | if $w_{d}<-2/3$; |
$Q_{\rm I}=\Gamma(\rho_{c}-\rho_{d})$ |  $A_{0}$, $A_{1}$ |  $A_{0}$, $A_{1}$ are attractors | for $w_{d}\geq-1$ and $\gamma>0$; |  $A_{1}:w_{d}\geq-1$, $\gamma>0$ | $q(A_{1})=\frac{1}{2}(1+3w_{d})<0$ | Figs. 1, 2, 3, 4.
|  $A_{2}$, $A_{3}$ | if $-2<w_{d}<-1$ and $\gamma>0$ |  $A_{1}$ is a global attractor for |  $A_{2}:$ always ($\gamma>0$) | if $w_{d}<-1/3$; |
|  |  |  $w_{d}\leq-2$ and $\gamma>0$ |  $A_{3}:-2\leq w_{d}\leq-1$, $\gamma>0$ |  $q(A_{2})=1/2$; $q(A_{3})=-1$ |
|  |  |  | For $\alpha>0$ and $\beta>0$ | $q(B_{0})=\frac{1}{2}\left(1+\frac{3\alpha w_{d}}{\alpha+\beta}\right)<0$ |
$Q_{\rm II}=\Gamma_{c}\rho_{c}-\Gamma_{d}\rho_{d}$ |  $B_{0}$, $B_{1}$ |  $B_{0}$, $B_{1}$ are attractors |  $B_{0}$ is a global attractor | $B_{0}:w_{d}\leq-\left(1+\frac{\beta}{\alpha}\right)$ | if $w_{d}<-\frac{1}{3}\left(1+\frac{\beta}{\alpha}\right)$; |
|  $B_{2}$, $B_{3}$ | for $-\left(1+\frac{\beta}{\alpha}\right)<w_{d}<-1$, | for $w_{d}\geq-1$, $\alpha>0$ and $\beta>0$; | $B_{1}:w_{d}\geq-1$ | $q(B_{1})=\frac{1}{2}(1+3w_{d})<0$ | Figs. 5, 6, 7, 8.
|  |  $\alpha>0$ and $\beta>0$ |  $B_{1}$ is a global attractor for |  $B_{2}:$ always | if $w_{d}<-1/3$; |
|  |  |  $w_{d}\leq-\left(1+\frac{\beta}{\alpha}\right)$, $\alpha>0$ and $\beta>0$ | $B_{3}:-\left(1+\frac{\beta}{\alpha}\right)\leq w_{d}\leq-1$ |  $q(B_{2})=1/2$; $q(B_{3})=-1$ |
|  |  |  | For $\gamma>0$ | $q(C_{0})=\frac{1}{2}\left(1+\frac{3w_{d}(3-\sqrt{5})}{2}\right)<0$ |
$Q_{\rm III}=\Gamma\left(\rho_{c}-\rho_{d}-\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}\right)$ |  $C_{0}$, $C_{1}$ |  $C_{0}$, $C_{1}$ are attractors |  $C_{0}$ is a global attractor | $C_{0}:w_{d}\leq-(3+\sqrt{5})/2$ | if $w_{d}<-(3+\sqrt{5})/6$; |
|  $C_{2}$, $C_{3}$ | for $-(3+\sqrt{5})/2<w_{d}<-1$ | for $w_{d}\geq-1$ and $\gamma>0$; | $C_{1}:w_{d}\geq-1$ | $q(C_{1})=\frac{1}{2}(1+3w_{d})<0$ | Figs. 9, 10, 11, 12.
|  | and $\gamma>0$ |  $C_{1}$ is a global attractor |  $C_{2}:$ always | if $w_{d}<-1/3$; |
|  |  | for $w_{d}\leq-(3+\sqrt{5})/2$ | $C_{3}:-(3+\sqrt{5})/2\leq w_{d}\leq-1$ |  $q(C_{2})=1/2$; $q(C_{3})=-1$ |
|  |  |  | For $\beta>\alpha>0$ |  |
$Q_{\rm IV}=\Gamma_{c}\rho_{c}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}$ |  $D_{0}$, $D_{00}$, |  $D_{0}$ is an attractor |  $D_{0}$ is a global attractor | $D_{0}:w_{d}\leq-\beta/\alpha$ | $q(D_{0})=\frac{1}{2}\left(1+\frac{3\alpha w_{d}}{\beta}\right)<0$ |
|  $D_{1}$, $D_{2}$, | if $-{\beta}/{\alpha}<w_{d}\leq-1$ | if $w_{d}>-1$ and $\beta>\alpha>0$; | $D_{00}:w_{d}\neq-1$ | if $w_{d}<-\frac{\beta}{3\alpha}$; | Figs. 13, 14, 15, 16.
| $D_{3}$ | and $\beta>\alpha>0$; |  $D_{1}$ is a global attractor | $D_{1}:w_{d}>-1$ | $q(D_{00})=\frac{1}{2}(1+3w_{d})<0$ |
|  |  $D_{1}$ is an attractor | if $w_{d}\leq-\beta/\alpha$ with $\beta>\alpha>0$ |  $D_{2}:$ always | if $w_{d}<-1/3$; |
|  | if $-\beta/\alpha<w_{d}<-1$ | and if $w_{d}<-1$ with $\alpha>\beta>0$; | $D_{3}:-\beta/\alpha\leq w_{d}<-1$ | $q(D_{1})=\frac{1}{2}(1+3w_{d})<0$ |
|  | and $\beta>\alpha>0$ |  $D_{00}$ is a global attractor | For $\alpha>\beta>0$ | if $w_{d}<-1/3$; |
|  |  | if $w_{d}>-1$ and $\alpha>\beta>0$ | $D_{00}:w_{d}<-1$ |  $q(D_{2})=1/2$; |
|  |  |  | $D_{1}:w_{d}>-1$ | $q(D_{3})=-1$ |
|  |  |  |  $D_{2}:$ always |  |
|  |  |  | For $\beta<\alpha<0$ |  |
$Q_{\rm V}=\Gamma_{d}\rho_{d}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}$ |  $E_{0}$, $E_{00}$ |  $E_{0}$, $E_{1}$ are attractors |  $E_{0}$ is a global attractor | $E_{0}:w_{d}\leq\beta/(\alpha-\beta)$ | $q(E_{0})=\frac{1}{2}\left(1+3w_{d}(1-\frac{\alpha}{\beta})\right)<0$ |
|  $E_{1}$, $E_{2}$, | if $\beta/(\alpha-\beta)<w_{d}<-1$ | if $w_{d}\geq-1$ with $\beta<\alpha<0$; |  $E_{00}:$ always | if $w_{d}<-\frac{1}{3}\frac{\beta}{\beta-\alpha}$; | Figs. 17, 18, 19.
| $E_{3}$ | with $\beta<\alpha<0$; |  $E_{1}$ is a global attractor | $E_{1}:w_{d}\geq-1$ |  $q(E_{00})=1/2$; |
|  |  $E_{00}$, $E_{1}$ are attractors | if $w_{d}\leq\beta/(\alpha-\beta)$ and $\beta<\alpha<0$; |  $E_{2}:$ always | $q(E_{1})=\frac{1}{2}(1+3w_{d})<0$ |
|  | if $w_{d}<-1$ with $\alpha<\beta<0$ |  $E_{00}$ is a global attractor | $E_{3}:\beta/(\alpha-\beta)\leq w_{d}\leq-1$ | if $w_{d}<-1/3$; |
|  |  | if $w_{d}\geq-1$ with $\alpha<\beta<0$ | For $\alpha<\beta<0$ |  $q(E_{2})=1/2$; |
|  |  |  | $E_{1}:w_{d}\geq-1$ | $q(E_{3})=-1$ |
|  |  |  |  $E_{2}:$ always |  |
|  |  |  | $E_{3}:w_{d}\leq-1$ |  |
Dynamical $w_{d}$
|  | Stability of Critical Points | Deceleration |
Interaction Models | Critical Points | Attractor | Global Attractor | Repeller | Parameter ($q$) | Figures
|  |  |  |  $\bar{A}_{1}:\nu<0$ and $\gamma>0$ |  $q(\bar{A}_{0})=-1/4$; |
$Q_{\rm I}=\Gamma(\rho_{c}-\rho_{d})$ |  $\bar{A}_{0}$, $\bar{A}_{1}$ |  $\bar{A}_{0}$, $\bar{A}_{1}$ are attractors |  $\bar{A}_{0}$ is a global attractor |  $\bar{A}_{2}:$ always ($\gamma>0$) |  $q(\bar{A}_{1})\longrightarrow-\infty$ if $\nu>0$ | Figs. 1, 2.
|  $\bar{A}_{2}$, $\bar{A}_{3}$ | if $\nu>0$ and $\gamma>0$ | if $\nu<0$ and $\gamma>0$ |  $\bar{A}_{3}:\nu>0$ and $\gamma>0$ | and $q(\bar{A}_{1})\longrightarrow+\infty$ if $\nu<0$; |
|  |  |  |  |  $q(\bar{A}_{2})=1/2$; $q(\bar{A}_{3})=-1$ |
|  |  |  |  | $q(\bar{B}_{0})=\frac{1}{2}\left(\frac{\beta-2\alpha}{\alpha+\beta}\right)<0$ |
$Q_{\rm II}=\Gamma_{c}\rho_{c}-\Gamma_{d}\rho_{d}$ |  $\bar{B}_{0}$, $\bar{B}_{1}$ |  $\bar{B}_{0}$, $\bar{B}_{1}$ are attractors |  $\bar{B}_{0}$ is a global attractor |  $\bar{B}_{1}:\nu<0$, $\alpha>0$ and $\beta>0$ | if $\beta<2\alpha$; | Figs. 5, 6.
|  $\bar{B}_{2}$, $\bar{B}_{3}$ | if $\nu>0$, $\alpha>0$ and $\beta>0$ | if $\nu<0$, $\alpha>0$ and $\beta>0$ |  $\bar{B}_{2}:$ always ($\alpha>0$, $\beta>0$) |  $q(\bar{B}_{1})\longrightarrow-\infty$ if $\nu>0$ |
|  |  |  |  $\bar{B}_{3}:$ $\nu>0$, $\alpha>0$ and $\beta>0$ | and $q(\bar{B}_{1})\longrightarrow+\infty$ if $\nu<0$; |
|  |  |  |  |  $q(\bar{B}_{2})=1/2$; $q(\bar{B}_{3})=-1$ |
|  |  |  |  |  $q(\bar{C}_{0})=\frac{3\sqrt{5}-7}{4}<0$; |
$Q_{\rm III}=\Gamma\left(\rho_{c}-\rho_{d}-\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}\right)$ |  $\bar{C}_{0}$, $\bar{C}_{1}$ |  $\bar{C}_{0}$, $\bar{C}_{1}$ are attractors |  $\bar{C}_{0}$ is a global attractor |  $\bar{C}_{1}:$ $\nu<0$ and $\gamma>0$ |  $q(\bar{C}_{1})\longrightarrow-\infty$ if $\nu>0$ | Figs. 9, 10.
|  $\bar{C}_{2}$, $\bar{C}_{3}$ | if $\nu>0$ and $\gamma>0$ | if $\nu<0$ and $\gamma>0$ |  $\bar{C}_{2}:$ always ($\gamma>0$) | and $q(\bar{C}_{1})\longrightarrow+\infty$ if $\nu<0$; |
|  |  |  |  $\bar{C}_{3}:$ $\nu>0$ and $\gamma>0$ |  $q(\bar{C}_{2})=1/2$; $q(\bar{C}_{3})=-1$ |
|  |  |  |  $\bar{D}_{00}:\nu\neq 0$ and $\beta>\alpha>0$ | $q(\bar{D}_{0})=\frac{\beta-3\alpha}{2\beta}<0$ |
$Q_{\rm IV}=\Gamma_{c}\rho_{c}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}$ |  $\bar{D}_{0}$, $\bar{D}_{00}$, |  $\bar{D}_{0}$, $\bar{D}_{1}$ are attractors |  $\bar{D}_{0}$ is a global attractor |  $\bar{D}_{1}:$ $\nu<0$ and $\beta>\alpha>0$ | if $\beta<3\alpha$; | Figs. 13, 14, 16.
|  $\bar{D}_{1}$, $\bar{D}_{2}$, | if $\nu>0$ and $\beta>\alpha>0$ | if $\nu<0$ and $\beta>\alpha>0$; |  $\bar{D}_{2}:$ always ($\beta>\alpha>0$) |  $q(\bar{D}_{00})=-1$; |
| $\bar{D}_{3}$ |  |  $\bar{D}_{00}$ is a global attractor |  $\bar{D}_{3}:$ $\nu>0$ and $\beta>\alpha>0$ |  $q(\bar{D}_{1})\longrightarrow-\infty$ if $\nu>0$ |
|  |  | if $\nu<0$ and $\alpha>\beta>0$; |  $\bar{D}_{00}:$ $\nu>0$ and $\alpha>\beta>0$ | and $q(\bar{D}_{1})\longrightarrow+\infty$ if $\nu<0$; |
|  |  |  $\bar{D}_{1}$ is a global attractor |  $\bar{D}_{1}:$ $\nu<0$ and $\alpha>\beta>0$ |  $q(\bar{D}_{2})=1/2$; |
|  |  | if $\nu>0$ and $\alpha>\beta>0$ |  $\bar{D}_{2}:$ always ($\alpha>\beta>0$) | $q(\bar{D}_{3})=-1$ |
|  |  |  |  $\bar{E}_{00}:$ always ($\beta<\alpha<0$) | $q(\bar{E}_{0})=\frac{3\alpha-2\beta}{2\beta}<0$ |
$Q_{\rm V}=\Gamma_{d}\rho_{d}-\Gamma_{cd}\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}}$ |  $\bar{E}_{0}$, $\bar{E}_{00}$ |  $\bar{E}_{0}$, $\bar{E}_{1}$ are attractors |  $\bar{E}_{0}$ is a global attractor |  $\bar{E}_{1}:$ $\nu<0$ with $\beta<\alpha<0$ | if $3\alpha>2\beta$; | Figs. 17, 18.
|  $\bar{E}_{1}$, $\bar{E}_{2}$, | if $\nu>0$ with $\beta<\alpha<0$; | if $\nu<0$ with $\beta<\alpha<0$; |  $\bar{E}_{2}:$ always ($\beta<\alpha<0$) |  $q(\bar{E}_{00})=1/2$; |
| $\bar{E}_{3}$ |  $\bar{E}_{00}$, $\bar{E}_{1}$ are attractors |  $\bar{E}_{00}$ is a global attractor |  $\bar{E}_{3}:$ $\nu>0$ with $\beta<\alpha<0$ |  $q(\bar{E}_{1})\longrightarrow-\infty$ if $\nu>0$ |
|  | if $\nu>0$ with $\alpha<\beta<0$ | if $\nu<0$ with $\alpha<\beta<0$ |  $\bar{E}_{1}:$ $\nu<0$ with $\alpha<\beta<0$ | and $q(\bar{E}_{1})\longrightarrow+\infty$ if $\nu<0$; |
|  |  |  |  $\bar{E}_{2}:$ always ($\alpha<\beta<0$) |  $q(\bar{E}_{2})=1/2$; |
|  |  |  |  $\bar{E}_{3}:$ $\nu>0$ with $\alpha<\beta<0$ | $q(\bar{E}_{3})=-1$ |

Table 6: Summary table describing all the interaction functions, the critical
points and their nature for both constant and dynamical DE equation of state.

####  III.5.2 Dynamical $w_{d}$

We work with the same dynamical $w_{d}$ of eqn. (29) for which the autonomous
system (86) takes the form:

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-\left(\frac{z}{1-z}\right)(1-x)(\alpha-\beta
x)-3x(1-x)\left(1+\nu\frac{(1-z)^{2}(1-x)}{z^{2}}\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)z\left(x-\nu\frac{(1-z)^{2}(1-x)^{2}}{z^{2}}\right),\end{array}\right.$
(92)

where $\nu=\frac{3AH_{0}^{2}}{\kappa^{2}}$. Regularizing the autonomous system
(92), we get

$\displaystyle\left\\{\begin{array}[]{ccc}x^{\prime}&=&-{z^{3}}(1-x)(\alpha-\beta
x)-3(1-z)x(1-x)\left(z^{2}+\nu(1-z)^{2}(1-x)\right),\\\
z^{\prime}&=&\frac{3}{2}(1-z)^{2}z\left(xz^{2}-\nu(1-z)^{2}(1-x)^{2}\right),\end{array}\right.$
(95)

where the regularization amounts to the result that the autonomous systems
(92) and (95) are topologically equivalent. The critical points of the system
(95) are

  * •

$\bar{E}_{0}=\left(\frac{\alpha}{\beta},1\right)$, $\bar{E}_{00}=(1,1)$,
$\bar{E}_{1}=(0,0)$, $\bar{E}_{2}=(1,0)$,
$S=\left\\{\left(x_{c},\frac{3x_{c}}{3x_{c}-(1-x_{c})(\alpha-\beta
x_{c})}\right)\right\\},$

where $S$ represents the set of critical points in which $x_{c}$ is a root of
$\Psi(x)\equiv 9x^{3}-\nu(1-x)^{4}(\alpha-\beta x)^{2}=0$.999Note that
$\Psi(x)$ is obtained from the following two nullclines:  $\displaystyle
xz^{2}-\nu(1-z)^{2}(1-x)^{2}=0,$ (96) $\displaystyle-z^{3}(1-x)[\alpha-\beta
x]-3x(1-x)(1-z)\left[z^{2}+\nu(1-z)^{2}(1-x)\right]=0.$ (97) As $\Psi(x)$
represents a six degree equation in $x$, therefore, $S$ may contain maximum
six critical points. Now we see that $\Psi(0)=-\nu\alpha^{2}<0$ (for $\nu>0$)
and $\Psi(1)>0$, therefore, from the Bolzano’s theorem [169], $\Psi(x)$ has at
least one real root in the interval $(0,1)$. For real and physically
meaningful critical points in our domain $R$, the $z$ component of the
critical point i.e., $z_{c}=\frac{3x_{c}}{3x_{c}-(1-x_{c})(\alpha-\beta
x_{c})}$ must belong to $[0,1]$. Thus, for any $0\leq x_{c}\leq 1$, the
criterion $0\leq\frac{3x_{c}}{3x_{c}-(1-x_{c})(\alpha-\beta x_{c})}\leq 1$
leads to the condition that $(\alpha-\beta x_{c})\leq 0$. Now, since
$\alpha<0$, $\beta<0$, we consider the following cases:

  * •

$\beta<\alpha<0$: In this case, we have $0\leq x_{c}\leq\alpha/\beta$. We also
see that $\Psi(\alpha/\beta)=9\left(\alpha/\beta\right)^{3}>0$. Consequently,
we conclude from the Bolzano’s theorem [169] that there is at least one root
of $\Psi(x)$ in $\left(0,\alpha/\beta\right)$. Now, from the derivative of
$\Psi(x)$,

$\displaystyle\Psi^{\prime}(x)$
$\displaystyle=-6{\beta}^{2}\nu(x-1)^{3}\left(x-\frac{\alpha}{\beta}\right)\left(x-\frac{2\alpha+\beta}{3\beta}\right)+27x^{2},$
(98)

we notice that $\Psi^{\prime}(x)>0$ for $x\in\left(0,\alpha/\beta\right)$.
Hence, the function $\Psi(x)$ is strictly increasing in
$\left(0,\alpha/\beta\right)$ which finally concludes that there is only one
root of $\Psi(x)$ in $\left(0,\alpha/\beta\right)$. Therefore, the set of
critical points $S$ contains only one critical point and for convenience we
label this critical point as $\bar{E}_{3}$. In this case, the point
$\bar{E}_{3}$ qualitatively behaves like the point $E_{3}$ which is described
earlier and correspondingly, the phase portrait looks same as the right plot
of Fig. 18.

On the other hand, for $\nu<0$, following the earlier arguments, we can show
that $\Psi(x)$ has no root in $(0,1)$. This shows that for $\nu<0$, we have
only four critical points: $\bar{E}_{0}$, $\bar{E}_{00}$, $\bar{E}_{1}$, and
$\bar{E}_{2}$ and the phase plot is same as the left plot of Fig. 17.

  * •

$\alpha<\beta<0$: We already know that $x_{c}$ satisfies the inequality
$\alpha-\beta x_{c}\leq 0$, but in contrary to the earlier case, in this
parameter space, we have $\alpha/\beta>1$ and consequently, the critical point
$\bar{E}_{0}$ does not belong to the physical domain $R$. We now investigate
the number of real roots of $\Psi(x)$ in the interval $[0,1]$ for $\nu>0$
because for $\nu<0$, as already commented, there is no root of $\Psi(x)$ in
$[0,1]$. Here, we see that $\frac{2\alpha+\beta}{3\beta}>1$.101010One can
check that
$\frac{2\alpha+\beta}{3\beta}=\frac{2}{3}(\frac{\alpha}{\beta})+\frac{1}{3}>\frac{2}{3}+\frac{1}{3}=1$
(since $\alpha/\beta>1$).  Now, looking at the expression for
$\Psi^{\prime}(x)$ in eqn. (98), we see that

    * –

$(x-1)<0$ as $x\in(0,1)$,

    * –

$\left(x-\frac{\alpha}{\beta}\right)<0$ in $(0,1)$ since
$\frac{\alpha}{\beta}>1$,

    * –

$\left(x-\frac{2\alpha+\beta}{3\beta}\right)<0$ in $(0,1)$ as
$\frac{2\alpha+\beta}{3\beta}>1$,

therefore, $\Psi^{\prime}(x)$ is always positive in $(0,1)$. Thus, $\Psi(x)$
being a strictly increasing function in $(0,1)$ has only one root in $(0,1)$.
Hence, the set $S$ contains only one critical point in the domain $R$. The
phase space stability analysis of this critical point is same as that of the
critical point $E_{3}$ and the phase plot is same as the left plot of Fig. 18.

As before, the polynomial $\Psi(x)$ has no real root in the interval $[0,1]$
for $\nu<0$. Therefore, in this case, we will have only three critical points:
$\bar{E}_{00}$, $\bar{E}_{1}$, and $\bar{E}_{2}$, and the phase plot will be
similar to the right plot of Fig. 17.

##  IV Summary and concluding remarks

Cosmology with non-gravitational interaction between DM and DE is the theme of
this work. This particular theory, according to the existing records, has
occupied a very decent place in the list of alternative cosmological models
beyond $\Lambda$CDM. In this article we have raised some important questions
regarding some not so usual constraints on the interaction functions and
performed a detailed phase space analysis of the interacting scenarios
featuring some novel qualities that distinguish with the existing works in
this direction. According to the existing records in the literature, in almost
every interacting scenarios, some common (but not so natural) assumptions are
considered, such as, the flow of energy should be either from DE to DM or in
the reverse direction, that means either DE will be gainer or DM will be
gainer. The question arises, why such unidirectional property of the
interaction function should be obeyed given the fact that the nature of the
interaction function is still an open question to the astrophysics and the
cosmology community? On the other hand, there are ceaseless debates on the
choice of the interaction functions $-$ whether the interaction functions
should involve the (global) expansion rate explicitly or not. Moreover, in the
context of the DE fluid, should we consider its equation of state to be
dynamical or constant?

Considering these debates, in this work we have considered the following set-
up: i) the assumption of unidirectional interaction functions have been
generalized by means of some sign shifting interaction functions which recover
the unidirectional interaction functions as a special case, and thus, in this
new picture of interacting dynamics, we allow the bidirectional energy flow;
ii) the interaction functions do not depend on the external parameters of the
universe, rather, they depend on the intrinsic nature of the dark components,
and hence, they are expected to offer inherent nature of the dark sector at
the fundamental level, and in addition, iii) along with the constant DE
equation of state, we have considered a dynamical parametrization of the DE
equation of state belonging to the class (28) which adds a new ingredient in
this context, and, so far we are aware of the existing literature, this is the
first time we are reporting such analysis.

We begin the study by considering a very simple but elegant linear interaction
function $Q_{\rm I}=\Gamma(\rho_{c}-\rho_{d})$ of (10) and then considered its
linear and nonlinear extensions in terms of other interaction functions given
in eqns. (11), (12), (13), (14). With the choice of suitable dimensionless
variables (this is very crucial in the analysis, since for some specific
choices of the variables, one may not obtain all the critical points of the
system), we have obtained all the critical points of the interacting
scenarios. The detailed analyses of the interacting scenarios for both
constant and dynamical $w_{d}$ are described in sections III.1, III.2, III.3,
III.4, III.5. Tables 1, 2, 3, 4, 5 summarize the critical points, their
existence, stability and the values of the key cosmological parameters and in
Figs. 1 $-$ 19, we have shown the nature of the critical points and the
evolution of some key cosmological parameters. For the sake of convenience, an
overall summary of results extracted out of all the sign shifting interacting
scenarios is given in Table 6 where mainly the nature of the critical points
for the present sign shifting interacting scenarios has been shown for
different nature of the DE equation of state.

Focusing on the constant DE equation of state, $w_{d}$, we found that, each
sign shifting interacting scenario admits a variety of critical points which
are qualitatively different, namely, the matter dominated critical point which
is unstable in nature; late time stable attractors corresponding to an
accelerating expansion of the universe in which one stable attractor is
completely DE dominated (i.e. $\Omega_{d}=1$) and in one attractor both DE and
DM exist and hence this attractor is physically more interesting according to
the present observational results. Moreover, we found that all the sign
shifting interacting scenarios also admit global attractors for different
regions of the DE equation of state, namely, $w_{d}>-1$ (non-phantom),
$w_{d}=-1$ (cosmological constant) and $w_{d}<-1$ (phantom), see Table 6. At
this point, it is important to mention that the sign shifting nature of the
interaction functions can affect the space of critical points. In
particularly, there is a connection between the late time stable attractors
and the present sign shifting interaction functions because if this sign
shifting nature of the present interaction functions is replaced by the
unidirectional interaction functions, that means when the transfer of energy
between the dark components is restricted to only in one direction111111If the
interaction functions $Q_{\rm I}=\Gamma(\rho_{c}-\rho_{d})$ and $Q_{\rm
III}=\Gamma(\rho_{c}-\rho_{d}-\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}})$,
are replaced by $\widetilde{Q}_{\rm I}=\Gamma(\rho_{c}+\rho_{d})$ and
$\widetilde{Q}_{\rm
III}=\Gamma(\rho_{c}+\rho_{d}+\frac{\rho_{c}\rho_{d}}{\rho_{c}+\rho_{d}})$
respectively (note that $\widetilde{Q}_{\rm I}$ and $\widetilde{Q}_{\rm III}$
represent unidirectional interaction functions), then the critical points
$A_{0}$, $C_{0}$ do not appear in these new unidirectional interacting
scenarios. In a similar fashion, for $Q_{\rm II}$, $Q_{\rm IV}$ and $Q_{\rm
V}$, if we impose that the coupling parameters will have the opposite signs
instead of the same signs which is essential for the sign shifting nature,
then the late time stable attractors $B_{0}$, $D_{0}$ and $E_{0}$ do not
appear in the respective unidirectional interacting scenario.  then the late
time stable attractors, namely, $A_{0}$, $B_{0}$, $C_{0}$, $D_{0}$ and $E_{0}$
do not appear within these unidirectional interaction functions. However, in
the proposed sign shifting interaction models, $A_{0}$, $B_{0}$, $C_{0}$,
$D_{0}$ and $E_{0}$ appear at the transitional point where
$Q(\rho_{c},\rho_{d})$ changes its sign. It is an interesting feature of the
proposed sign shifting interacting models and it is a subject for further
investigations.

The case with dynamical DE equation of state presents a more general
interacting scenario. The present choice of $w_{d}$ (eqn. (29)) covers both
the phantom (for $A>0$, equivalently, $\nu>0$) and non-phantom (for $A<0$,
equivalently, $\nu<0$) regimes. In terms of the number of critical points,
each sign shifting interacting model with dynamical phantom case (i.e.
$\nu>0$) has one extra critical point compared to the corresponding sign
shifting interacting model with dynamical quintessence case ($\nu<0$). Again,
similar to the constant $w_{d}$ case, here too, we observe that if the sign
shifting nature of the interaction functions is replaced by the unidirectional
interaction functions, the late time stable attractors, namely, $\bar{A}_{0}$,
$\bar{B}_{0}$, $\bar{C}_{0}$, $\bar{D}_{0}$ and $\bar{E}_{0}$ do not appear in
this case, but in the context of the proposed sign shifting interaction models
these late time stable attractors arise at the transitional point, that means
where $Q(\rho_{c},\rho_{d})$ changes its sign. Specifically, we have the
following observations:

  1. 1.

Dynamical phantom: In all the sign shifting interacting dynamical phantom
scenarios, we find one matter dominated era (unstable in nature) representing
a decelerating phase, late time stable attractors in which one attractor is
completely DE dominated ($\Omega_{d}=1$) and other attractor allows the
concurrence of DE and DM. Moreover, we noticed that only Model IV in this
series admits one global attractor $\bar{D}_{1}$ provided that the
dimensionless coupling parameters satisfy $\alpha>\beta>0$. According to the
existing literature on the interacting dynamical phantom scenarios [18, 19,
147, 155] concurrent existence of the matter dominated, only DE dominated
($\Omega_{d}=1$), and the co-existence of DE and DM ($\Omega_{d}\neq 0$,
$\Omega_{c}\neq 0$), as we observed within the context of present sign
shifting interacting models, is very rare. For example, even though some
specific interacting models exhibit the matter dominated phase [147, 155] but
the simultaneous occurrence of the DE dominated stable attractor
($\Omega_{d}=1$) and the stable late time scaling attractor corresponding to
an accelerating phase of the universe has not been found.

  2. 2.

Dynamical quintessence: We find that all the scenarios in this category admit
the matter dominated phase which is unstable in nature and it corresponds to a
past decelerating phase. But, unlike in the phantom interacting scenario, here
only one late time stable attractor is allowed which is global in nature (see
Table 6) and this critical point allows the existence of both DE and DM. The
existence of the matter dominated phase within the present sign shifting
models is interesting because such phase is not so common in a variety of
interaction models when $w_{d}$ lies in the quintessence regime, see for
instance [19, 27, 146, 148, 153].

Based on the outcomes of the present article, it is evident that the sign
shifting interaction models are quite appealing. The results further emphasize
that there should not have any particular reason to prefer only the
unidirectional interaction functions in the context of interacting DE, rather,
the bidirectional interaction functions are quite promising and they deserve
further attention. Specifically, the analysis with dynamical $w_{d}$ within
these interaction models is very promising but such analysis is rare in the
literature.

##  V Acknowledgments

The authors thank the referees for their time to read our article and for
giving some useful comments that helped us to improve the article. SH
acknowledges the financial support from the University Grants Commission
(UGC), Govt. of India (NTA Ref. No: 201610019097). JdH is supported by the
Spanish grant PID2021-123903NB-I00 funded by MCIN/AEI/10.13039/501100011033
and by “ERDF A way of making Europe”. TS and SP acknowledge the financial
support from the Department of Science and Technology (DST), Govt. of India
under the Scheme “Fund for Improvement of S&T Infrastructure (FIST)” (File No.
SR/FST/MS-I/2019/41).

## References

  * Riess  _et al._ [1998] A. G. Riess _et al._ (Supernova Search Team), Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ [1999] S. Perlmutter _et al._ (Supernova Cosmology Project), Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae, Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Copeland  _et al._ [2006] E. J. Copeland, M. Sami, and S. Tsujikawa, Dynamics of dark energy, Int. J. Mod. Phys. D 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Frieman  _et al._ [2008] J. Frieman, M. Turner, and D. Huterer, Dark Energy and the Accelerating Universe, Ann. Rev. Astron. Astrophys. 46, 385 (2008), arXiv:0803.0982 [astro-ph] .
  * Bamba  _et al._ [2012] K. Bamba, S. Capozziello, S. Nojiri, and S. D. Odintsov, Dark energy cosmology: the equivalent description via different theoretical models and cosmography tests, Astrophys. Space Sci. 342, 155 (2012), arXiv:1205.3421 [gr-qc] .
  * Nojiri and Odintsov [2006] S. Nojiri and S. D. Odintsov, Introduction to modified gravity and gravitational alternative for dark energy, eConf C0602061, 06 (2006), arXiv:hep-th/0601213 .
  * Nojiri and Odintsov [2011] S. Nojiri and S. D. Odintsov, Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models, Phys. Rept. 505, 59 (2011), arXiv:1011.0544 [gr-qc] .
  * De Felice and Tsujikawa [2010] A. De Felice and S. Tsujikawa, f(R) theories, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Capozziello and De Laurentis [2011] S. Capozziello and M. De Laurentis, Extended Theories of Gravity, Phys. Rept. 509, 167 (2011), arXiv:1108.6266 [gr-qc] .
  * Clifton  _et al._ [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified Gravity and Cosmology, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Koyama [2016] K. Koyama, Cosmological Tests of Modified Gravity, Rept. Prog. Phys. 79, 046902 (2016), arXiv:1504.04623 [astro-ph.CO] .
  * Cai  _et al._ [2016] Y.-F. Cai, S. Capozziello, M. De Laurentis, and E. N. Saridakis, f(T) teleparallel gravity and cosmology, Rept. Prog. Phys. 79, 106901 (2016), arXiv:1511.07586 [gr-qc] .
  * Nojiri  _et al._ [2017] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Akrami  _et al._ [2021] Y. Akrami _et al._ (CANTATA),  _Modified Gravity and Cosmology: An Update by the CANTATA Network_, edited by E. N. Saridakis, R. Lazkoz, V. Salzano, P. Vargas Moniz, S. Capozziello, J. Beltrán Jiménez, M. De Laurentis, and G. J. Olmo (Springer, 2021) arXiv:2105.12582 [gr-qc] .
  * Bahamonde  _et al._ [2023] S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. Levi Said, J. Mifsud, and E. Di Valentino, Teleparallel gravity: from theory to cosmology, Rept. Prog. Phys. 86, 026901 (2023), arXiv:2106.13793 [gr-qc] .
  * Amendola [2000] L. Amendola, Coupled quintessence, Phys. Rev. D 62, 043511 (2000), arXiv:astro-ph/9908023 .
  * Amendola and Quercellini [2003] L. Amendola and C. Quercellini, Tracking and coupled dark energy as seen by WMAP, Phys. Rev. D 68, 023514 (2003), arXiv:astro-ph/0303228 .
  * Guo  _et al._ [2005] Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Cosmological evolution of interacting phantom energy with dark matter, JCAP 05, 002, arXiv:astro-ph/0412624 .
  * Gumjudpai  _et al._ [2005] B. Gumjudpai, T. Naskar, M. Sami, and S. Tsujikawa, Coupled dark energy: Towards a general description of the dynamics, JCAP 06, 007, arXiv:hep-th/0502191 .
  * Cai and Wang [2005] R.-G. Cai and A. Wang, Cosmology with interaction between phantom dark energy and dark matter and the coincidence problem, JCAP 03, 002, arXiv:hep-th/0411025 .
  * Wang  _et al._ [2006] B. Wang, C.-Y. Lin, and E. Abdalla, Constraints on the interacting holographic dark energy model, Phys. Lett. B 637, 357 (2006), arXiv:hep-th/0509107 .
  * Barrow and Clifton [2006] J. D. Barrow and T. Clifton, Cosmologies with energy exchange, Phys. Rev. D 73, 103520 (2006), arXiv:gr-qc/0604063 .
  * Berger and Shojaei [2006] M. S. Berger and H. Shojaei, An Interacting Dark Energy Model for the Expansion History of the Universe, Phys. Rev. D 74, 043530 (2006), arXiv:astro-ph/0606408 .
  * Setare [2007] M. R. Setare, Interacting holographic dark energy model and generalized second law of thermodynamics in non-flat universe, JCAP 01, 023, arXiv:hep-th/0701242 .
  * Rosenfeld [2007] R. Rosenfeld, Reconstruction of interacting dark energy models from parameterizations, Phys. Rev. D 75, 083509 (2007), arXiv:astro-ph/0701213 .
  * Feng  _et al._ [2007] C. Feng, B. Wang, Y. Gong, and R.-K. Su, Testing the viability of the interacting holographic dark energy model by using combined observational constraints, JCAP 09, 005, arXiv:0706.4033 [astro-ph] .
  * Boehmer  _et al._ [2008] C. G. Boehmer, G. Caldera-Cabral, R. Lazkoz, and R. Maartens, Dynamics of dark energy with a coupling to dark matter, Phys. Rev. D 78, 023505 (2008), arXiv:0801.1565 [gr-qc] .
  * Olivares  _et al._ [2008] G. Olivares, F. Atrio-Barandela, and D. Pavon, The Integrated Sachs-Wolfe Effect in Interacting Dark Energy Models, Phys. Rev. D 77, 103520 (2008), arXiv:0801.4517 [astro-ph] .
  * Valiviita  _et al._ [2008] J. Valiviita, E. Majerotto, and R. Maartens, Instability in interacting dark energy and dark matter fluids, JCAP 07, 020, arXiv:0804.0232 [astro-ph] .
  * Mohseni Sadjadi and Vadood [2008] H. Mohseni Sadjadi and N. Vadood, Notes on interacting holographic dark energy model in a closed universe, JCAP 08, 036, arXiv:0806.2767 [gr-qc] .
  * Cruz  _et al._ [2008] N. Cruz, S. Lepe, F. Pena, and J. Saavedra, Holographic dark energy interacting with dark matter in a Closed Universe, Phys. Lett. B 669, 271 (2008), arXiv:0807.3854 [gr-qc] .
  * Wang and Zhang [2008] S. Wang and Y. Zhang, Alleviation of Cosmic Age Problem In Interacting Dark Energy Model, Phys. Lett. B 669, 201 (2008), arXiv:0809.3627 [astro-ph] .
  * Caldera-Cabral  _et al._ [2009a] G. Caldera-Cabral, R. Maartens, and L. A. Urena-Lopez, Dynamics of interacting dark energy, Phys. Rev. D 79, 063518 (2009a), arXiv:0812.1827 [gr-qc] .
  * Gavela  _et al._ [2009] M. B. Gavela, D. Hernandez, L. Lopez Honorez, O. Mena, and S. Rigolin, Dark coupling, JCAP 07, 034, [Erratum: JCAP 05, E01 (2010)], arXiv:0901.1611 [astro-ph.CO] .
  * Caldera-Cabral  _et al._ [2009b] G. Caldera-Cabral, R. Maartens, and B. M. Schaefer, The Growth of Structure in Interacting Dark Energy Models, JCAP 07, 027, arXiv:0905.0492 [astro-ph.CO] .
  * Valiviita  _et al._ [2010] J. Valiviita, R. Maartens, and E. Majerotto, Observational constraints on an interacting dark energy model, Mon. Not. Roy. Astron. Soc. 402, 2355 (2010), arXiv:0907.4987 [astro-ph.CO] .
  * Lopez Honorez  _et al._ [2010] L. Lopez Honorez, B. A. Reid, O. Mena, L. Verde, and R. Jimenez, Coupled dark matter-dark energy in light of near Universe observations, JCAP 09, 029, arXiv:1006.0877 [astro-ph.CO] .
  * Yu  _et al._ [2010] F. Yu, J. Zhang, J. Lu, W. Wang, and Y. Gui, A more general interacting model of holographic dark energy, Phys. Lett. B 688, 263 (2010), arXiv:1004.2092 [astro-ph.CO] .
  * Wei [2010] H. Wei, Revisiting the Cosmological Constraints on the Interacting Dark Energy Models, Phys. Lett. B 691, 173 (2010), arXiv:1004.0492 [gr-qc] .
  * Duran  _et al._ [2010] I. Duran, D. Pavon, and W. Zimdahl, Observational constraints on a holographic, interacting dark energy model, JCAP 07, 018, arXiv:1007.0390 [astro-ph.CO] .
  * Chen  _et al._ [2011] X. Chen, B. Wang, N. Pan, and Y. Gong, Constraining the interacting dark energy models from weak gravity conjecture and recent observations, Phys. Lett. B 695, 30 (2011), arXiv:1008.3455 [astro-ph.CO] .
  * Cao  _et al._ [2011] S. Cao, N. Liang, and Z.-H. Zhu, Testing the phenomenological interacting dark energy with observational $H(z)$ data, Mon. Not. Roy. Astron. Soc. 416, 1099 (2011), arXiv:1012.4879 [astro-ph.CO] .
  * Clemson  _et al._ [2012] T. Clemson, K. Koyama, G.-B. Zhao, R. Maartens, and J. Valiviita, Interacting Dark Energy – constraints and degeneracies, Phys. Rev. D 85, 043007 (2012), arXiv:1109.6234 [astro-ph.CO] .
  * Marulli  _et al._ [2012] F. Marulli, M. Baldi, and L. Moscardini, Clustering and redshift-space distortions in interacting dark energy cosmologies, Mon. Not. Roy. Astron. Soc. 420, 2377 (2012), arXiv:1110.3045 [astro-ph.CO] .
  * Baldi and Salucci [2012] M. Baldi and P. Salucci, Constraints on interacting dark energy models from galaxy Rotation Curves, JCAP 02, 014, arXiv:1111.3953 [astro-ph.CO] .
  * He  _et al._ [2011] J.-H. He, B. Wang, and E. Abdalla, Deep connection between f(R) gravity and the interacting dark sector model, Phys. Rev. D 84, 123526 (2011), arXiv:1109.1730 [gr-qc] .
  * Avelino and da Silva [2012] P. P. Avelino and H. M. R. da Silva, Effective dark energy equation of state in interacting dark energy models, Phys. Lett. B 714, 6 (2012), arXiv:1201.0550 [astro-ph.CO] .
  * Harko and Lobo [2013] T. Harko and F. S. N. Lobo, Irreversible thermodynamic description of interacting dark energy-dark matter cosmological models, Phys. Rev. D 87, 044018 (2013), arXiv:1210.3617 [gr-qc] .
  * Sun and Yue [2013] C.-Y. Sun and R.-H. Yue, Stable large-scale perturbations in interacting dark-energy model, JCAP 08, 018, arXiv:1303.0684 [astro-ph.CO] .
  * Chimento  _et al._ [2013] L. P. Chimento, M. G. Richarte, and I. E. Sánchez García, Interacting dark sector with variable vacuum energy, Phys. Rev. D 88, 087301 (2013), arXiv:1310.5335 [gr-qc] .
  * Li and Zhang [2014] Y.-H. Li and X. Zhang, Large-scale stable interacting dark energy model: Cosmological perturbations and observational constraints, Phys. Rev. D 89, 083009 (2014), arXiv:1312.6328 [astro-ph.CO] .
  * Yang and Xu [2014a] W. Yang and L. Xu, Cosmological constraints on interacting dark energy with redshift-space distortion after Planck data, Phys. Rev. D 89, 083517 (2014a), arXiv:1401.1286 [astro-ph.CO] .
  * Yang and Xu [2014b] W. Yang and L. Xu, Testing coupled dark energy with large scale structure observation, JCAP 08, 034, arXiv:1401.5177 [astro-ph.CO] .
  * Yang and Xu [2014c] W. Yang and L. Xu, Coupled dark energy with perturbed Hubble expansion rate, Phys. Rev. D 90, 083532 (2014c), arXiv:1409.5533 [astro-ph.CO] .
  * Li  _et al._ [2014] Y.-H. Li, J.-F. Zhang, and X. Zhang, Parametrized Post-Friedmann Framework for Interacting Dark Energy, Phys. Rev. D 90, 063005 (2014), arXiv:1404.5220 [astro-ph.CO] .
  * Tamanini [2015] N. Tamanini, Phenomenological models of dark energy interacting with dark matter, Phys. Rev. D 92, 043524 (2015), arXiv:1504.07397 [gr-qc] .
  * Goncalves  _et al._ [2015] R. S. Goncalves, G. C. Carvalho, and J. S. Alcaniz, Low-z test for interacting dark energy, Phys. Rev. D 92, 123504 (2015), arXiv:1507.01921 [astro-ph.CO] .
  * Pan  _et al._ [2015] S. Pan, S. Bhattacharya, and S. Chakraborty, An analytic model for interacting dark energy and its observational constraints, Mon. Not. Roy. Astron. Soc. 452, 3038 (2015), arXiv:1210.0396 [gr-qc] .
  * Yang  _et al._ [2016] W. Yang, H. Li, Y. Wu, and J. Lu, Cosmological constraints on coupled dark energy, JCAP 10, 007, arXiv:1608.07039 [astro-ph.CO] .
  * Nunes  _et al._ [2016] R. C. Nunes, S. Pan, and E. N. Saridakis, New constraints on interacting dark energy from cosmic chronometers, Phys. Rev. D 94, 023508 (2016), arXiv:1605.01712 [astro-ph.CO] .
  * Yang  _et al._ [2017a] W. Yang, N. Banerjee, and S. Pan, Constraining a dark matter and dark energy interaction scenario with a dynamical equation of state, Phys. Rev. D 95, 123527 (2017a), [Addendum: Phys.Rev.D 96, 089903 (2017)], arXiv:1705.09278 [astro-ph.CO] .
  * Dutta  _et al._ [2017] J. Dutta, W. Khyllep, and N. Tamanini, Scalar-Fluid interacting dark energy: cosmological dynamics beyond the exponential potential, Phys. Rev. D 95, 023515 (2017), arXiv:1701.00744 [gr-qc] .
  * Di Valentino  _et al._ [2017] E. Di Valentino, A. Melchiorri, and O. Mena, Can interacting dark energy solve the $H_{0}$ tension?, Phys. Rev. D 96, 043503 (2017), arXiv:1704.08342 [astro-ph.CO] .
  * Santos  _et al._ [2017] L. Santos, W. Zhao, E. G. M. Ferreira, and J. Quintin, Constraining interacting dark energy with CMB and BAO future surveys, Phys. Rev. D 96, 103529 (2017), arXiv:1707.06827 [astro-ph.CO] .
  * Mifsud and Van De Bruck [2017] J. Mifsud and C. Van De Bruck, Probing the imprints of generalized interacting dark energy on the growth of perturbations, JCAP 11, 001, arXiv:1707.07667 [astro-ph.CO] .
  * Yang  _et al._ [2017b] W. Yang, S. Pan, and D. F. Mota, Novel approach toward the large-scale stable interacting dark-energy models and their astronomical bounds, Phys. Rev. D 96, 123508 (2017b), arXiv:1709.00006 [astro-ph.CO] .
  * Linton  _et al._ [2018] M. S. Linton, A. Pourtsidou, R. Crittenden, and R. Maartens, Variable sound speed in interacting dark energy models, JCAP 04, 043, arXiv:1711.05196 [astro-ph.CO] .
  * Yang  _et al._ [2018a] W. Yang, S. Pan, and J. D. Barrow, Large-scale Stability and Astronomical Constraints for Coupled Dark-Energy Models, Phys. Rev. D 97, 043529 (2018a), arXiv:1706.04953 [astro-ph.CO] .
  * Yang  _et al._ [2017c] W. Yang, L. Xu, H. Li, Y. Wu, and J. Lu, Testing the Interacting Dark Energy Model with Cosmic Microwave Background Anisotropy and Observational Hubble Data, Entropy 19, 327 (2017c).
  * Yang  _et al._ [2018b] W. Yang, S. Pan, E. Di Valentino, R. C. Nunes, S. Vagnozzi, and D. F. Mota, Tale of stable interacting dark energy, observational signatures, and the $H_{0}$ tension, JCAP 09, 019, arXiv:1805.08252 [astro-ph.CO] .
  * von Marttens  _et al._ [2019] R. von Marttens, L. Casarini, D. F. Mota, and W. Zimdahl, Cosmological constraints on parametrized interacting dark energy, Phys. Dark Univ. 23, 100248 (2019), arXiv:1807.11380 [astro-ph.CO] .
  * Yang  _et al._ [2018c] W. Yang, A. Mukherjee, E. Di Valentino, and S. Pan, Interacting dark energy with time varying equation of state and the $H_{0}$ tension, Phys. Rev. D 98, 123527 (2018c), arXiv:1809.06883 [astro-ph.CO] .
  * Li  _et al._ [2020] C. Li, X. Ren, M. Khurshudyan, and Y.-F. Cai, Implications of the possible 21-cm line excess at cosmic dawn on dynamics of interacting dark energy, Phys. Lett. B 801, 135141 (2020), arXiv:1904.02458 [astro-ph.CO] .
  * Mishra and Chakraborty [2019] S. Mishra and S. Chakraborty, Stability analysis of an interacting holographic dark energy model, Mod. Phys. Lett. A 34, 19 (2019), arXiv:2005.06296 [gr-qc] .
  * Pan  _et al._ [2019a] S. Pan, W. Yang, E. Di Valentino, E. N. Saridakis, and S. Chakraborty, Interacting scenarios with dynamical dark energy: Observational constraints and alleviation of the $H_{0}$ tension, Phys. Rev. D 100, 103520 (2019a), arXiv:1907.07540 [astro-ph.CO] .
  * Di Valentino  _et al._ [2020] E. Di Valentino, A. Melchiorri, O. Mena, and S. Vagnozzi, Interacting dark energy in the early 2020s: A promising solution to the $H_{0}$ and cosmic shear tensions, Phys. Dark Univ. 30, 100666 (2020), arXiv:1908.04281 [astro-ph.CO] .
  * von Marttens  _et al._ [2020] R. von Marttens, L. Lombriser, M. Kunz, V. Marra, L. Casarini, and J. Alcaniz, Dark degeneracy I: Dynamical or interacting dark energy?, Phys. Dark Univ. 28, 100490 (2020), arXiv:1911.02618 [astro-ph.CO] .
  * Cheng  _et al._ [2020] G. Cheng, Y.-Z. Ma, F. Wu, J. Zhang, and X. Chen, Testing interacting dark matter and dark energy model with cosmological data, Phys. Rev. D 102, 043517 (2020), arXiv:1911.04520 [astro-ph.CO] .
  * Pan  _et al._ [2020a] S. Pan, G. S. Sharov, and W. Yang, Field theoretic interpretations of interacting dark energy scenarios and recent observations, Phys. Rev. D 101, 103533 (2020a), arXiv:2001.03120 [astro-ph.CO] .
  * Pan  _et al._ [2020b] S. Pan, J. de Haro, W. Yang, and J. Amorós, Understanding the phenomenology of interacting dark energy scenarios and their theoretical bounds, Phys. Rev. D 101, 123506 (2020b), arXiv:2001.09885 [gr-qc] .