# Interacting dark sector: a dynamical system perspective

Chonticha Kritpetch<EMAIL_ADDRESS>High Energy Physics and Cosmology
Research Group, School of Science, University of Phayao, Phayao 56000,
Thailand  Nandan Roy<EMAIL_ADDRESS>(Corresponding Author) Centre
for Theoretical Physics and Natural Philosophy,
Nakhon Sawan Studiorum for Advanced Studies, Mahidol University, Nakhonsawan,
60130, Thailand  Narayan Banerjee<EMAIL_ADDRESS>Department of
Physical Sciences, Indian Institute of Science Education and Research Kolkata,
Mohanpur 741246,WB, India

###### Abstract

We investigate the interaction between the dark sectors from a dynamical
systems perspective. A general setup for interacting dark energy models that
incorporates both quintessence and phantom fields through a switch parameter,
allowing an interaction in the dark sectors has been considered. In the first
part of our analysis, we have not assumed any specific form of the interaction
and in the second part, we invoked examples in general framework of the
interaction. The potentials of the scalar field are classified into two broad
classes of potentials: exponential and non-exponential. We identify the
potential late-time attractors of the system which have a complete dark energy
domination. From our analysis, it is evident there could be an interaction
between the dark sector. The interaction, if any, weakens over time. We find
for the quintessence field the transfer of energy from dark matter to dark
energy can flip the direction and on the contrary, for the phantom field, it
is only from dark matter to dark energy.

##  I Introduction

For the past two decades, various cosmological observations have provided
substantial evidence of a universe expanding at an accelerated rate[1, 2, 3,
4, 5]; however, the explanation of this behavior still remains a challenge.
Although the cosmological constant is the simplest and most successful
candidate for dark energy, the driver of the acceleration, it faces
significant challenges such as that of the huge discrepancy of its
observationally required value with the theoretically predicted one and the
coincidence problem[6]. Recently from high-precision cosmological data a
statistically significant discrepancy in the estimation of the current value
of the Hubble parameter ($H_{0}$) between early-time and late-time
observations has been reported. Explanation of this tension is another open
problem in cosmology that challenges the cosmological constant. Data from the
early universe measurement estimate $H_{0}\sim(67.0-68.5)$ km/s/Mpc [7, 8, 9],
while the measurement of the $H_{0}$ observing the local universe using the
distance ladder measurements reported $H_{0}=(74.03\pm 1.42)$ km/s/Mpc [10,
11, 12, 13]. This tension in the measurement of the $H_{0}$ indicates the
possibility of involvement of new physics while explaining the accelerated
expansion of the universe.

Dynamical dark energy models are considered to be the alternatives to the
cosmological constant which can alleviate these problems. Different such
models such as quintessence, k-essence, phantom dark energy, etc., have been
proposed [14, 15]. Generally in these models, the accelerated expansion of the
universe is caused by a scalar field rolling down the potential and generating
sufficient negative pressure[16, 17, 18, 19, 20, 21, 22].

In dynamical dark energy models the possibility of interaction between the
dark matter and dark energy was introduced in order to alleviate the cosmic
coincidence problem[23, 24, 25, 26, 27]. Recently these models have shown the
potential to alleviate the $H_{0}$ and $\sigma_{8}$ tensions as well[28, 29,
30, 31, 32, 33, 34, 35]. With the interaction, the dark matter and the dark
energy are not conserved independently rather they are conserved together.

The presence of interaction between the dark sectors can affect the overall
evolution of the universe[26]. A wide variety of investigations are already
there in the literature that look at the effect of the interaction in the dark
sector on different cosmological observables [36, 37, 24, 38, 39, 40, 41, 42].
The techniques of dynamical systems analysis have been extensively used to
study the stability and late-time behavior of different interacting dark
energy models[43, 44, 45, 46, 47, 48, 49].

In this work, we have performed the dynamical system analysis of the
interacting scalar field dark energy model from a general perspective without
considering any particular form of the interaction as far as possible. Our
setup considers both the quintessence and the phantom scalar field in a single
analysis. Also, the present work considers a general class of potentials,
broadly, an exponential form or a non-exponential one. To test our approach,
we have proposed a general class of interaction to start with that can
encompass a wide range of forms of interactions suggested in the literature.
We have conducted analytical and numerical investigations to assess the
stability of the system, track its evolution, and compare the model with the
observed data set.

The broad result of this work is that the models normally settle with a final
dark energy dominated scenario; a non-gravitational interaction, if there is
any, dies down with the evolution and the direction of the transfer of energy
due to the interaction may have a flip for the quintessence case but not such
for the phantom models.

The paper is organized in the following way, in Sec: II the mathematical setup
of the system is discussed. Section III deals with the investigation of the
system with exponential potential and Section IV deals with the investigation
of the non-exponential potential. In Section V we have summarized our results
and findings.

##  II Mathematical Setup

Let us consider a universe that is spatially flat, where all the matter
components of the universe follow the barotropic relation given by
$p_{i}=w_{i}\rho_{i}$, where $p_{i}$ represents the pressure, $\rho_{i}$
represents the density and $w_{i}$ represents the equation of state (EoS) of a
component.

In such a universe, if the dark energy is considered to be a scalar field, the
Einstein field equations can be written as follows:

$3H^{2}=\kappa^{2}\left(\sum_{i}\rho_{i}+\frac{1}{2}\epsilon\dot{\phi}^{2}+V(\phi)\right),$
(1)

$\dot{H}=-\frac{\kappa^{2}}{2}\left(\sum_{i}\left(\rho_{i}+P_{i}\right)+\epsilon\dot{\phi}^{2}\right).$
(2)

We introduce the switch parameter $\epsilon$ to incorporate both quintessence
and phantom field dynamics within a single framework. The equation corresponds
to a quintessence field when $\epsilon=+1$ and a phantom field when
$\epsilon=-1$. In a scenario where dark matter and dark energy interact with
each other such that the density of each component is conserved together but
not individually, the continuity equation for each component can be written as
follows,

$\displaystyle\dot{\rho}_{m}+3H\rho_{m}=-Q,$ (3)
$\displaystyle\dot{\rho}_{\phi}+3H\left(\rho_{\phi}+P_{\phi}\right)=Q,$ (4)

where the components of dark matter and dark energy are identified by
subscripts $m$ and $\phi$, respectively. The interaction term is denoted by
$Q$. If $Q$ is positive the transfer of energy happens from dark matter to
dark energy and vice versa.

In this interacting dark sector scenario Klein-Gordon equation of the scalar
field can be written as;

$\ddot{\phi}+3H\dot{\phi}+\epsilon\frac{dV}{d\phi}=\epsilon\frac{Q}{\dot{\phi}}.$
(5)

To perform a dynamical system analysis, we introduce the following set of
dimensionless variables,

$x^{2}=\frac{\kappa^{2}\dot{\phi}^{2}}{6H^{2}},y^{2}=\frac{\kappa^{2}V(\phi)}{3H^{2}},\lambda=-\frac{1}{V(\phi)}\frac{dV(\phi)}{d\phi}.$
(6)

With these, the system reduces to the following set of autonomous equations,

$\displaystyle x^{\prime}=-3x+\sqrt{3/2}\epsilon\lambda
y^{2}+\frac{3}{2}x\left(1+\epsilon x^{2}-y^{2}\right)+\epsilon f(x,y),$ (7a)
$\displaystyle y^{\prime}=-\sqrt{3/2}\lambda xy+\frac{3}{2}y\left(1+\epsilon
x^{2}-y^{2}\right),$ (7b)
$\displaystyle\lambda^{\prime}=-\sqrt{6}\lambda^{2}(\Gamma-1)x,$ (7c)

where
$\Gamma=\frac{V(\phi)\frac{\partial^{2}V(\phi)}{\partial\phi^{2}}}{\left(\frac{\partial
V(\phi)}{\partial\phi}\right)^{2}}$. We consider a general form of the
interaction term as $Q=\sqrt{6}\dot{\phi}H^{2}f(x,y)$, where $f(x,y)$
represents a wide range of functions involving the dynamical variables $x$ and
$y$. This formulation allows for the incorporation of different types of
interactions that have been studied in cosmology. For a comprehensive list of
various choices of $f(x,y)$, see [49] and the references therein. From the
physical point of view in a purely dark matter or dark energy-dominated state
the interaction and its derivatives should vanish, that is
$f(x,y)=0=\frac{\partial f(x,y)}{\partial x}=\frac{\partial f(x,y)}{\partial
y}$. The dynamical system variables can be used to express various
cosmological parameters in a concise manner as,

$\displaystyle\Omega_{\phi}=\epsilon x^{2}+y^{2},$ (8) $\displaystyle
w_{\phi}=\frac{\epsilon x^{2}-y^{2}}{\epsilon x^{2}+y^{2}},$ (9)
$\displaystyle q=-1+\frac{3}{2}(1+\epsilon x^{2}-y^{2}).$ (10)

In order to close the autonomous system given in equation (7), it is necessary
to specify a particular form of the $\Gamma$ function, which can be
essentially related to choosing a specific form of the potential. Depending on
the choice of the $\Gamma$ we can classify the system into two classes. The
first class arises when $\Gamma=1$, resulting in an exponential potential that
effectively reduces the system to a 2-dimensional form. On the other hand, the
second class corresponds to $\Gamma\neq 1$, which corresponds to any potential
except the exponential potential. This classification has already been used in
[50, 51].

In the next section, we discuss the fixed points of the system and the
corresponding stability of those fixed points. Fixed points are obtained by
simultaneously solving the autonomous equations given in eq.(7) with
$x^{\prime}=y^{\prime}=\lambda^{\prime}=0$.

##  III Exponential Potential

In this case, we assume that the potential has an exponential form and set
$\Gamma=1$. This reduces the dimension of the system from 3D to 2D since
$\lambda$ is constant.

###  III.1 Fixed point and stability

The fixed points of this system are given in Table.1 and in the last column we
tabulated the form of the interaction $f(x,y)$ at the fixed points. There are
only three classes of fixed points.

Fixed Points | $x$ | $y$ | $f(x,y)$
---|---|---|---
$Q_{1}$ | 0 | 0 | 0
$Q_{2i}$ |  | $\frac{3}{2}x(1+\epsilon x^{2})-3x+\epsilon f(x,0)=0$
---
0 |  |  $f(x,y)=0$ for $x=0,x\pm 1$
---
$f(x,y)\neq 0$ for $0<x<1$
$Q_{3i}$ | $\frac{1}{2}\left(\epsilon\left(2f(x,y)+\sqrt{6}\lambda\right)+2\sqrt{6}\lambda x^{2}-2x\left(\lambda^{2}\epsilon+3\right)\right)=0$ | $\pm\frac{\sqrt{3x^{2}\epsilon-\sqrt{6}\lambda x+3}}{\sqrt{3}}$ | $f(x,y)$

Table 1: List of the fixed points.

####  III.1.1 Fixed Point $Q_{1}$

The fixed point $Q_{1}$ corresponds to a completely dark matter-dominated
($\Omega_{\phi}=0$) regime and the interaction term $f(x,y)$ vanishes at this
fixed point. The eigenvalues associated with $Q_{1}$ are given by
($\frac{3}{2},-\frac{3}{2}$). This fixed point is a saddle in nature.

####  III.1.2 Fixed Point $Q_{2}$

Depending on the choice of the function $f(x,y)$, the class of the fixed
points $Q_{2}$ can be a single or multiple fixed point upon solving the
equation given in Table.1,

$\frac{3}{2}x(1+\epsilon x^{2})-3x+\epsilon f(x,0)=0$ (11)

The eigenvalues associated with this fixed point are given by
$\frac{1}{2}\left(3x^{2}\epsilon-\sqrt{6}\lambda
x+3\right),\frac{1}{2}\left(2\epsilon\partial_{x}f+9x^{2}\epsilon-3\right)$.
For a completely matter-dominated fixed point ($\Omega_{m}=1$) for which
$x=0$, one gets back the $Q_{1}$ as the fixed point. For the complete
quintessence field domination, $(x=\pm 1)$ these fixed points are unstable
fixed points since the eigenvalues reduce to
$[\frac{1}{2}(6\pm\sqrt{6}\lambda),3]$. From the complete phantom field
domination one of the fixed points is stable and another is unstable depending
on the choice of $\lambda$ as the eigenvalues for this case are
$[\frac{1}{2}(\pm\sqrt{6}\lambda),-6]$. Here we have used the fact that
$\partial_{x}f=0$ for a complete dark energy domination.

For cases where both dark matter and dark energy contribute, the stability
conditions are as follows. For the quintessence field ($\epsilon=+1$),

$\displaystyle\left(-1\leq
x<0;\lambda<\frac{3x^{2}+3}{\sqrt{6}x};\partial_{x}f<\frac{1}{2}\left(3-9x^{2}\right)\right),$
(12a) $\displaystyle\left(0<x\leq
1;\lambda>\frac{3x^{2}+3}{\sqrt{6}x};\partial_{x}f<\frac{1}{2}\left(3-9x^{2}\right)\right),$
(12b)

and for the phantom field ($\epsilon=-1$),

$\displaystyle\left(x<0;\lambda<\frac{3-3x^{2}}{\sqrt{6}x};\partial_{x}f>\frac{1}{2}\left(-9x^{2}-3\right)\right),$
(13a)
$\displaystyle\left(x>0;\lambda>\frac{3-3x^{2}}{\sqrt{6}x};\partial_{x}f>\frac{1}{2}\left(-9x^{2}-3\right)\right).$
(13b)

Contrary to the canonical field ($0\leq\Omega_{\phi}\leq 1$), there is no
strict positivity condition on the energy density of the phantom field. It can
be negative ($\Omega_{\phi}<0$) too. The existence of these fixed points
previously has been reported in [52, 53]

####  III.1.3 Fixed Point $Q_{3}$

The general form of the class of fixed points $Q_{3}$ is given in the form of
an equation.

$\frac{1}{2}\left(\epsilon\left(2f(x,y)+\sqrt{6}\lambda\right)+2\sqrt{6}\lambda
x^{2}-2x\left(\lambda^{2}\epsilon+3\right)\right)=0.$

One can solve the above algebraic equation for a given form of the interaction
$f(x,y)$ to find all associated fixed points. Although the eigenvalues
corresponding to this point can be quite complicated in form, for a 2D system,
one can use the trace and determinant of the Jacobian matrix to investigate
the stability. The trace $T_{Q3}$ and determinant $D_{Q3}$ of the Jacobian
matrix at these fixed points are given in appendix A. The condition for
stability of this fixed point is $T_{Q3}<0$ and $D_{Q3}>0$.

As an example, let us consider a special case where $x=0$ and $y=\pm 1$
depicting a completely dark energy-dominated universe and hence $f(x,y)=0$ and
$\partial_{x}f=\partial_{y}f=0$, the trace and determinant of the fixed point
reduce to the following simple form:

$\displaystyle T_{Q3}$ $\displaystyle=-6,$ $\displaystyle D_{Q3}$
$\displaystyle=3\epsilon\lambda^{2}+9.$

For the quintessence field($\epsilon=+1$), the fixed point in this example is
stable. On the contrary for the phantom field ($\epsilon=-1$), it depends on
the choice of $\lambda$. Since we have chosen $x=0$ and $y=\pm 1$, the field
has a slow roll, and $\lambda<<1$. Therefore, even for the phantom field, this
fixed point is an attractor.

###  III.2 Example

Inorder to investigate the general setup further, here, we propose a general
form of interaction as an example to test our approach:

$f(x,y)=\alpha(1-\epsilon x^{2}-y^{2})^{m}x^{\gamma}$ (14)

This particular form of interaction allows us to examine a wide range of
interactions [49, 53, 46, 54, 55]. In table.3 of appendix:B we have given a
list of interactions that are used in literature and can be incorporated into
the above general form. It is important to note that this form is not limited
to only those specific interactions mentioned in the table.3.

For this choice of the $f(x,y)$ the fixed point $Q_{1}$ exists only when
$\gamma>0$, as the condition $f(x,y)=0$ is necessary for the existence of this
fixed point. Any choice of interaction that violates this criterion would miss
the possibility of a pure matter-dominated universe as an initial condition.
Furthermore, this fixed point is inherently unstable for both the quintessence
and phantom field, regardless of the specific form of the interaction chosen.

The particular fixed points included in $Q_{2}$ class can be found by solving
the quadratic equation;

$\frac{3}{2}x(1+\epsilon x^{2})+\epsilon\alpha(1-\epsilon
x^{2})^{m}x^{\gamma}-3x=0.$ (15)

For both the quintessence and phantom case, there are multiple solutions to
the above equations leading to multiple fixed points, with $x=0$ being the
trivial solution. The fixed point $x=0$ is indistinguishable from the matter-
dominated case $Q_{1}$. A complete dark energy domination is represented by
the fixed point $x^{2}=1$. It has been already shown in the general analysis
in the previous subsection that for the complete quintessence field
domination, these fixed points are stable, and for the complete phantom field
domination one of the fixed points is stable another is unstable.

Depending on the choice of $\gamma,m$, there could be fixed points that can
represent a state of the universe where there are both dark matter and dark
energy contributions. From our general analysis in the previous subsection,
these points could be unstable or saddle in nature. For the quintessence
field, this particular interaction renders the stability conditions in the
expressions 12 to the following,

$\displaystyle\left(-1\leq
x<0;\lambda<\frac{3x^{2}+3}{\sqrt{6}x};f<\frac{3x(1-x^{2})(1-3x^{2})}{\gamma(1-x^{2})-2x^{2}}\right)$
$\displaystyle\left(0<x\leq
1;\lambda>\frac{3x^{2}+3}{\sqrt{6}x};f<\frac{3x(1-x^{2})(1-3x^{2})}{\gamma(1-x^{2})-2x^{2}}\right).$

Similarly, for the phantom field, conditions in the expressions 13 reduce to;

$\displaystyle\left(x<0;\lambda<\frac{3-3x^{2}}{\sqrt{6}x};f>-\frac{3x(1+x^{2})(1+3x^{2})}{\gamma(1+x^{2})+2x^{2}}\right)$
$\displaystyle\left(x>0;\lambda>\frac{3-3x^{2}}{\sqrt{6}x};f>-\frac{3x(1+x^{2})(1+3x^{2})}{\gamma(1+x^{2})+2x^{2}}\right)$

.

Here we have used $\partial_{x}f(x,y)=f(x,y)(\frac{\gamma}{x}-\frac{2\epsilon
x}{1-\epsilon x^{2}-y^{2}})$ and $y=0$ at this fixed point.

To obtain all the fixed points belonging to this $Q_{3}$ class one needs to
solve the equation,

$\begin{split}\frac{1}{2}\Biggl{(}\epsilon\left(\sqrt{6}\lambda+2\alpha
x^{\gamma}\left(1-x^{2}\epsilon\right)^{m}\right)\\\ +2\sqrt{6}\lambda
x^{2}-2x\left(\lambda^{2}\epsilon+3\right)\Biggr{)}=0\end{split}$ (16)

There could be multiple fixed points depending on the choice of $\gamma,m$.
One can easily compute the trace and determinant given in the appendix A and
find the stability of these fixed points. We then use numerical techniques to
find the phase space behaviour and evolution of the system for different
choices of the model parameters.

In Fig.1 we have shown the phase plot of the system with the exponential
potential for different choices of the $\lambda,\gamma$, and $m$ parameters
for the quintessence field. Here we have considered $\alpha=-0.2$. This choice
is motivated by the posterior obtained for the interaction parameter $\alpha$
using a similar mathematical setup in [56]. The first row (blue background)
and the third row (yellow background) represent $\lambda=-1$ and $\lambda=+1$
respectively. The middle row (orange background) represents the $\lambda=0$
case. All the cases indicate a completely dark energy-dominated universe as a
late-time attractor. Notice that with the change in the sign of $\lambda$ the
sign of the value of the $x$ at the late time attractor also changes. Positive
$\lambda$ corresponds to a positive value of $x$ and vice versa.

Also, these plots agree with the analytical finding that for $\gamma=0$ case
there would not be a fixed point that is completely dominated by dark matter
component. The first two columns correspond to the $\gamma=0$ case where there
is no fixed point at $x=0,y=0$ on the other hand the last two columns where
$\gamma=1$, $x=0,y=0$ clearly seen to have a saddle fixed point. The circle
represents $0\leq\Omega_{\phi}\leq 1$ where all the physical dynamics happen.

For the phantom field in Fig.2 we have shown the phase plot of the system for
the same choice of the parameters as in the quintessence case for the region
$0\leq\Omega_{\phi}\leq 1$. In this region, the phantom cases have only two
classes of fixed points as those belonging to the $Q_{2i}$ coincide with the
$Q_{1}$. Similar to the quintessence case phantom field also has complete dark
energy domination as the late-time attractor. The trajectories outside the
region $0\leq\Omega_{\phi}\leq 1$ enter the region as there are no attractors
outside the region [52].

For further investigation of the evolution of the system, we have numerically
evolved the autonomous system given in Eq.(7). For the numerical solution of
the system, either one needs to supply the initial condition or the current
condition of the $x,y$. Here we consider the second approach where we have
estimated the current values of the $x,y$ from the observation by solving,
$\Omega_{\phi 0}=\epsilon x_{0}^{2}+y_{0}^{2}$ and
$q_{0}=-1+\frac{3}{2}(1+\epsilon x_{0}^{2}-y_{0}^{2})$ where $\Omega_{\phi 0}$
and $q_{0}$ are the current density parameter of the scalar field and the
current value of the deceleration parameter respectively. Here we consider
$\Omega_{\phi 0}=0.68$[7] and $q_{0}=-0.51$[57]. For quintessence field
($\epsilon=+1$) we estimate $x_{0}=0.09,y_{0}=0.825$ and for the phantom field
$x_{0}=0.01,y=0.824$. In Fig.3 we have shown the evolution of different
cosmological variables. We have chosen $\lambda=0.5,\alpha=-0.3$ for the
quintessence field and $\lambda=0.2,\alpha=-0.3$ for the phantom field. The
choice of $\lambda$ is arbitrary but the particular choices we have here is to
avoid difficulties in the numerical integration and also to fit data. For the
choice of $\alpha$ parameter, we have used the constraints obtained from
recent cosmological observation using a similar mathematical setup[56]. In
Fig.3(a) we have shown the evolution of the density parameters
$\Omega_{\phi}$,$\Omega_{m}$ and $f(x,y)$. Evolution of the density parameters
$\Omega_{\phi}$ and $\Omega_{m}$ have the expected behavior whereas $f(x,y)$
(in dashed line) shows some intriguing nature. In the distant past, the
magnitude of the interaction $f(x,y)$ was much larger when compared to the
present epoch and approaches zero asymptotically for the future. For all the
cases the maximum value of the interaction is during the matter domination and
starts to decrease as the dark energy gradually dominates. To understand the
evolution of the interaction better in the figure.3(d) we have presented the
phase plot of the interaction term $f(x,y)$ vs $f^{\prime}(x,y)$ using the
eq.(14) and the eq.(7). It is interesting to note that the evolution of the
$f(x,y)$ started from a non-zero value and evolved to zero as the universe
will be completely dark energy-dominated. Another interesting fact to notice
here for all the cases the interaction term $f(x,y)$ has a flip in signature.

In Fig.3(b) we have plotted the evolution of the EoS of the scalar field
together with the deceleration parameter $q$ and in Fig.3(c) we have plotted
the evolution of the Hubble parameter $H(N)$ with respect to $N$. For
comparison with the observational data, we have also shown the measurement
from various observations. The data set used in this plot is given in
Appendix.C. As it is evident these models of interaction can fit the data
quite well, particularly at the late time.

In fig.4 we present the phantom case for the same choice of $\lambda$ and
$\alpha$ parameters. $\Omega_{m}$ and $\Omega_{\phi}$ show the expected
behavior. The evolution of $f(x,y)$ has a behavior similar to the quintessence
case. The interaction is found to have the maximum in magnitude during the
matter domination and decrease gradually to zero in the future which is
completely dominated by the phantom field. But if we notice the evolution of
$f(x,y)$ from the phase plot given in fig.4(d), there is a qualitative
difference with the quintessence case. For all the cases we have considered
the evolution of the interaction is unidirectional for the phantom case since
there is no flip in the signature of the $f(x,y)$, and the interaction
vanishes faster than the quintessence case.

The evolution of the $w_{\phi}$ and the deceleration parameter $q$ for the
phantom field is shown in fig.3(d). As it is expected $w_{\phi}$ evolves from
$w_{\phi}<-1$ to $w_{\phi}=-1$ at present. The comparison with the
observational data is shown in the fig.4(c) by plotting $H(N)$ vs $N$ and it
can be seen the data can be fitted very well even with the phantom field as
dark energy.

##  IV Non-Exponential Potential

Here we present the phase space behavior of the class of potentials that are
non-exponential, characterized by the condition $\Gamma\neq 1$. The fixed
points associated with this class of potentials are listed in Table 2. In
total, there are four classes of fixed points.

Fixed Points | $x$ | $y$ | $\lambda$ | $f(x,y)$
---|---|---|---|---
$P_{1}$ | 0 | 0 | $\lambda$ | 0
$P_{2i}$ | $\frac{3}{2}x(1-\epsilon x^{2})-\epsilon f(x,0)=0$ | 0 | 0 |  | $f(x,y=0$
---
for $x=0,\pm 1$
$P_{3\pm}$ | 0 | $\pm 1$ | $\lambda=-\sqrt{\frac{2}{3}}f(x,y)$ | $f(x,y)$
$P_{4\pm}$ |  | $3x-\epsilon f(x,y)=0$
---
$y^{2}=(1+\epsilon x^{2})$ | 0 | $f(x,y)$

Table 2: List of the fixed points.

####  IV.0.1 Fixed Point $P_{1}$

Fixed Point $P_{1}$ represents a completely matter-dominated situation for
both the quintessence and the phantom fields ($\epsilon=\pm 1$). The
corresponding eigenvalues of this fixed point are
($\frac{3}{2},0,-\frac{3}{2}$). Regardless of the specific form of the
interaction function $f(x,y)$, this fixed point is inherently saddle in
nature.

####  IV.0.2 Fixed Point $P_{2}$

The cosmological behaviour of the fixed point $P_{2}$ is similar to the fixed
point $Q_{2}$ in the exponential case except that $\lambda=0$ for $P_{2}$. The
particular fixed points corresponding to this class can be obtained from
Eq.(11). The eigenvalues corresponding to these fixed points are
($0,\frac{3}{2}\left(x^{2}\epsilon+1\right),\frac{1}{2}\left(2\epsilon\partial_{x}f+9x^{2}\epsilon-3\right)$).
These fixed points are nonhyperbolic fixed points but can be easily checked
for the quintessence field ($\epsilon=+1$) the second eigenvalue can not be
negative hence for the quintessence field these fixed points are unstable. For
the phantom field, the second eigenvalue is $\frac{3}{2}(1-x^{2})$ which can
be written as $\frac{3}{2}(1+\Omega_{\phi_{P_{2}}})$ as
$\Omega_{\phi_{P_{2}}}=-x^{2}$ at this fixed point. This fixed point is
unstable for $-1<\Omega_{\phi}\leq 1$. For complete dark matter or dark energy
domination the interaction is zero as it is expected but during the epoch when
there is contribution from both the dark matter and dark energy the
interaction may not vanish.

####  IV.0.3 Fixed Point $P_{3}$

The fixed point $P_{3}$ represents a completely dark energy-dominated
situation where the value of the $\lambda$ depends on the choice of the form
of the interaction. Since this is a completely dark energy-dominated fixed
point so from a physical point of view there could not be any interaction
between the dark energy and the dark matter because of the absence of the
latter. Hence we consider the interaction term $f(x,y)=\partial_{x}f(x,y)=0$,
therefore $\lambda=0$ at this fixed point.

The eigenvalues corresponding to this fixed point are $(-3,0,-3)$. This is a
nonhyperbolic fixed point so one can not use the linear stability analysis. A
more complex analytical tool like the central manifold theorem or the
numerical tools has to be used to analyse the stability of this fixed point.

####  IV.0.4 Fixed Point $P_{4}$

The fixed point $P_{4\pm}$ are intermediate fixed points for which there can
be contributions from both the dark energy and dark matter. The corresponding
eigenvalues at this fixed point are the following;

$\begin{aligned}
&\left\\{0,\frac{1}{6}\left(-A+3\epsilon\partial_{x}f-18\right),\frac{1}{6}\left(A+3\epsilon\partial_{x}f-18\right)\right\\}\end{aligned}$,
where

$A=\sqrt{3}\sqrt{4f^{2}(\partial_{x}f-3\epsilon)+4\epsilon f\sqrt{9+\epsilon
f^{2}}\partial_{y}f+3(\partial_{x}f)^{2}}$.

This fixed point is also nonhyperbolic and more sophisticated numerical or
mathematical methods should be used to study its stability with a particular
choice of interaction.

###  IV.1 Example

For further investigation of the non-exponential potential case here we
consider the same form of the interaction given in Eq.(14).

Like the exponential case, fixed point $P_{1}$ only exists for $\gamma>0$
models and it is unstable, independent of the form of the interaction.

All the fixed points belonging to $P_{2}$ class can be found by solving the
quadratic equation given in Eq.(15). Irrespective of any particular form of
the interaction this fixed point is unstable for both the quintessence and
phantom field.

The fixed point $P_{3}$ exists for any choice of $\gamma$ and $m$ and it is a
completely dark energy-dominated state. The eigenvalues at this fixed point
are ($-3,0,3$) as $\partial_{x}f=0$ at this fixed point. This fixed point is
non-hyperbolic so its stability can not be understood using linear stability
analysis. We shall investigate it numerically for some specific choices of the
model parameters $\gamma,m$, and $\Gamma$.

To get the particular fixed point belonging to the $P_{4}$ class one needs to
solve the following equation,

$3x-\epsilon\alpha x^{\gamma}\left(1-\epsilon x^{2}-y^{2}\right)^{m}=0.$

One can notice from another equation corresponding to this fixed point
$y^{2}=(1+\epsilon x^{2})$, for the quintessence for any $x\neq 0,y^{2}>1$.
Hence the only physical solution to this equation is $x=0$ which makes this
fixed point indistinguishable from $P_{3}$. For the phantom field, the
behaviour is richer as this fixed point can represent complete matter
domination to dark energy domination, and a combination of both depending on
the choice of $\gamma,m$. This fixed point is also nonhyperbolic. In the next,
we shall numerically investigate the phase space behaviour and the evolution
of this system for this particular choice of interaction.

For the numerical investigation, we chose $\Gamma=1/2$ for which the potential
becomes $V(\phi)=(A+B\phi)^{2}$. In Fig.5 and Fig.6 we have plotted the phase
diagram for $\Gamma=1/2$ considering the same combination of the $\gamma,m$ as
in the exponential case for the quintessence and phantom field respectively.

In Fig.5, we showed the 3D phase plot for the quintessence model (for
$\Gamma=\frac{1}{2}$). The physically allowed region confined by Friedmann
constraint takes the shape of a cylinder described with $\lambda$ axis as the
axis of the cylinder which is shown as the shaded cylindrical region in the
plot. The surface of the cylinder represents the complete dark energy
domination whereas the axis of the cylinder represents the complete matter
domination. For the quintessence case, it can be noticed for all the models we
have considered the solutions originated from the points $x=\pm 1,y=0$
111Though the solutions originated from the dark energy domination is not
expected from the physical point of view it is due to not considering the
contribution of radiation in our analysis since we are interested in late-time
dynamics. As in the early times, the contribution from the matter is
negligible the scalar field works as a proxy to the radiation to fulfill the
Friedmann constraint equation $\Omega_{m}+\Omega_{\phi}=1$. This argument can
be validated from the evolution of the EoS of the scalar field in Fig:7(b)
where the EoS of the quintessence field remains positive during the early dark
energy domination hence unable to drive the acceleration of the universe. This
has been reported before in [51] . and are attracted towards complete dark
energy domination where there could be a contribution from the both potential
and kinetic part of the field. The whole $\lambda$ axis represents the matter
domination as expected from our analytical finding remains saddle for all
cases.

In Fig.6 we have shown the 3D phase plot of the phantom case for $\Gamma=1/2$.
The blue-shade hyperbolic regions show the Friedmann constraint region for the
phantom field. Compared to the quintessence case the phase space for the
phantom case is more involved. The solutions can start from the region where
$\Omega_{\phi}<0$ but enter the Friedmann constraint region and are attracted
towards completely phantom field dominated late time attractors at the
$\lambda=0$ plane.

For further investigation of the system’s evolution, we have considered the
same strategy as the exponential case. The numerical simulation of the system
has been considered with the estimated current values of $x_{0},y_{0}$ as in
the exponential case. In Fig.7 and Fig.8, we have shown the evolution of the
cosmological parameters for the quintessence and phantom fields for the choice
of $\Gamma=1/2$. For both the quintessence and the phantom case the evolution
of the energy density $\Omega_{\phi}$ and $\Omega_{m}$ shows expected behavior
and the the interaction starts from a non-zero value and vanishes with time as
the late time attractor is always dark energy dominated. The phase plot
behavior of the interaction for the quintessence and the phantom field in
Fig.7(d) and Fig.8(d) show a behavior similar to the exponential potential
counterpart. For the quintessence case, there is a flip in the sign of the
$f(x,y)$ indicating a change in the direction of the flow of energy from dark
matter to dark energy. However, for the phantom case, it is unidirectional
from dark matter to dark energy only. For the quintessence field $w_{\phi}>-1$
at an early time and approaches to $w_{\phi}\simeq-1$ at late times and for
the phantom field $w_{\phi}<-1$ at an early time and approaches to
$w_{\phi}\simeq-1$. It can be also seen from the plots Fig.7(c) and Fig.8(c)
both the quintessence and the phantom field with the nonexponential potentials
can fit the $H(z)$ vs $z$ data set given in Appendix C, particularly at late
time.

##  V Conclusions

This study uses dynamical system analysis to investigate interacting dark
energy models, including both quintessence and phantom scalar field in a
single setup through a switch parameter $\epsilon$. At first, the equations of
motion of the scalar fields are recast to a set of autonomous systems by
considering suitable variable transformations. Also, a general approach has
been adopted for the choice of the potential. The choice of potential is
classified into two general classes exponential potential and non exponential
potential. Stability analysis has been performed without considering any
particular form of the interaction both for the exponential and non-
exponential classes of the potentials.

A general form of interaction has been proposed as an example that can
incorporate a wide class of popular forms of interaction. The numerical
evolution of the system considering this form of interaction has been studied
for both the exponential and non-exponential potentials. We have compared the
evolution of the Hubble parameter in these models against the observed data
and it can be seen that these models can fit the data very well, particularly
at the late time.

From our analysis, it is evident that the late-time attractor is a completely
dark energy-dominated state of the universe. The numerical evolution of the
universe suggests that there can be some interaction between the dark sectors,
although this interaction becomes weaker with time and becomes negligible as
the universe becomes more and more dark energy dominated. Our finding shows
that for the quintessence field, the interaction can be from dark matter to
dark energy and vice versa. During the evolution of the universe, the
interaction started from dark matter to dark energy and it reversed its
direction since there is a flip in the sign of the interaction term $f(x,y)$
and vanishes in the future. For the phantom field, the interaction is from
dark matter to dark energy and vanishes over time. This generalized method of
studying the interaction can be applied to study other dark energy models,
which will be presented in future work.

###### Acknowledgements.

This work (Grant No. RGNS 65-132) was supported by Office of the Permanent
Secretary, Ministry of Higher Education, Science, Research and Innovation (OPS
MHESI), Thailand Science Research and Innovation (TSRI), and university of
Phayao. CK would like to acknowledge that this work was partially supported by
the University of Phayao’s International Sabbatical Leave Fund 2022. In
addition, the authors would like to thank Burin Gumjudpai and Manabendra
Sharma for their useful discussion on the work.

##  Appendix A Eigenvalues of the $Q_{3}$ fixed point.

The trace and determinant of this fixed point are;

$\begin{split}T_{Q3}&=\epsilon\partial_{x}f+3\sqrt{\frac{3}{2}}\lambda x-6,\\\
D_{Q3}&=-3x^{2}\left(-5\lambda^{2}+M-3\epsilon\right)\\\ &+\sqrt{6}\lambda
x\left(M\epsilon+\lambda^{2}(-\epsilon)-9\right)\\\ &-3M\epsilon-
Nx\sqrt{9x^{2}\epsilon-3\sqrt{6}\lambda x+9}\\\ &+\frac{\lambda
N\epsilon\sqrt{3x^{2}\epsilon-\sqrt{6}\lambda x+3}}{\sqrt{2}}\\\
&-6\sqrt{6}\lambda x^{3}\epsilon+3\lambda^{2}\epsilon+9\end{split}$ (17)

The condition to have a stable region $T_{Q3}<0$ and $D_{Q3}>0$.

##  Appendix B Form of interactions

In the following, we give a list of popular forms of interactions already
studied in the literature that can be incorporated into our general
parametrization of the interaction given in Eq.14,

$Q$ | $q$ | References
---|---|---
$\beta\rho\dot{\phi}$ | $\frac{\sqrt{3}}{\sqrt{2}}\beta\left(1-x^{2}-y^{2}\right)$ | [58, 59, 46]
$\beta H\rho$ | $\frac{\beta}{2}\left(1-x^{2}-y^{2}\right)/x$ | [60, 54]
$\beta H\dot{\phi}^{2}$ | $\beta x$ | [55, 61]
$\beta\rho\dot{\phi}^{2}/H$ | $3\beta x\left(1-x^{2}-y^{2}\right)$ | [54]
$\alpha\left(\rho+\rho_{\phi}\right)\dot{\phi}$ | $\alpha/\sqrt{6}$ | [62]
Table 3: A list of example interactions that can be incorporated in the
general form of the interaction considered in this work.

##  Appendix C  Hubble’s parameter vs. redshift & scale factor.

$z$ | $N=\ln(\frac{1}{1+z})$ | $H(z)\left(\frac{\text{km/s}}{\text{Mpc}}\right)$ | Ref.
---|---|---|---
0.07 | -0.067 | $69\pm 19.6$ | [63]
0.09 | -0.086 | $69\pm 12$ | [64]
0.100 | -0.095 | $69\pm 12$ | [64]
0.120 | -0.113 | $68.6\pm 26.2$ | [63]
0.170 | -0.157 | $83\pm 8$ | [64]
0.179 | -0.164 | $75\pm 4$ | [65]
0.199 | -0.181 | $75\pm 5$ | [65]
0.200 | -0.182 | $72.9\pm 29.6$ | [63]
0.270 | -0.239 | $77\pm 14$ | [64]
0.280 | -0.246 | $88.8\pm 36.6$ | [63]
0.320 | -0.277 | $79.2\pm 5.6$ | [66]
0.352 | -0.301 | $83\pm 14$ | [65]
0.3802 | -0.322 | $83\pm 13.5$ | [65]
0.400 | -0.336 | $95\pm 17$ | [64]
0.4004 | -0.336 | $77\pm 10.2$ | [65]
0.4247 | -0.353 | $87.1\pm 11.2$ | [65]
0.440 | -0.364 | $82.6\pm 7.8$ | $[82]$
0.4497 | -0.371 | $92.8\pm 12.9$ | [65]
0.470 | -0.385 | $89\pm 50$ | [83]

$z$ | $N=\ln(\frac{1}{1+z})$ | $H(z)\left(\frac{\text{km/s}}{\text{Mpc}}\right)$ | Ref.
---|---|---|---
0.4783 | -0.390 | $80.9\pm 9$ | [65]
0.480 | -0.392 | $97\pm 62$ | [84]
0.570 | -0.451 | $100.3\pm 3.7$ | [66]
0.593 | -0.465 | $104\pm 13$ | [65]
0.600 | -0.470 | $87.9\pm 6.1$ | $[82]$
0.680 | -0.518 | $92\pm 8$ | [65]
0.730 | -0.548 | $97.3\pm 7$ | $[82]$
0.781 | -0.577 | $105\pm 12$ | [65]
0.875 | -0.628 | $125\pm 17$ | [65]
0.880 | -0.631 | $90\pm 40$ | [84]
0.900 | -0.641 | $117\pm 23$ | [64]
1.037 | -0.711 | $154\pm 20$ | [65]
1.300 | -0.832 | $168\pm 17$ | [64]
1.363 | -0.859 | $160\pm 33.6$ | [85]
1.430 | -0.887 | $177\pm 18$ | [64]
1.530 | -0.928 | $140\pm 14$ | [64]
1.750 | -1.011 | $202\pm 40$ | [64]
1.965 | -1.086 | $186.5\pm 50.4$ | $[85]$
2.340 | -1.205 | $222\pm 7$ | $[86]$
2.360 | -1.211 | $226\pm 8$ | 887

Table 4: Observastional $H(z)$ vs $z$ with references which has been used in
this work.  Figure 1: Phase plot of $x$ vs $y$ for the exponential potential
($\Gamma=1$) for different choices of the $\lambda,\gamma,m$ for the
quintessence field. Figure 2: Phase plot of $x$ vs $y$ for the exponential
potential for different choices of the $\lambda,\gamma,m$ for the phantom
field.

(a)

(b)

(c)

(d)

Figure 3: Plot of the different cosmological parameters for the quintessence
field for the exponential potential ($\Gamma=1$).

(a)

(b)

(c)

(d)

Figure 4: Plot of the different cosmological parameters for the phantom field
for exponential potential $(\Gamma=1)$. Figure 5: Phase plot of $x$ vs $y$ for
the non-exponential potential($\Gamma=\frac{1}{2}$) for different choices of
the $\gamma,m$ for the Quintessence field. Figure 6: Phase plot of $x$ vs $y$
for the non-exponential potential($\Gamma=2$) for different choices of the
$\gamma,m$ for the Quintessence field.

(a)

(b)

(c)

(d)

Figure 7: Plot of the different cosmological parameters for the quintessence
field for the $\Gamma=1/2$.

(a)

(b)

(c)

(d)

Figure 8: Plot of the different cosmological parameters for the phantom field
for the $\Gamma=1/2$.

## References

  * [1] Adam G. Riess et al.  Observational evidence from supernovae for an accelerating universe and a cosmological constant.  Astron. J., 116:1009–1038, 1998.
  * [2] S. Perlmutter et al.  Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae.  Astrophys. J., 517:565–586, 1999.
  * [3] Attila Meszaros.  On the Reality of the accelerating universe.  Astrophys. J., 580:12–15, 2002.
  * [4] M. Arnaud et al.  Planck intermediate results. XXXI. Microwave survey of Galactic supernova remnants.  Astron. Astrophys., 586:A134, 2016.
  * [5] Christopher P Ahn, Rachael Alexandroff, Carlos Allende Prieto, Scott F Anderson, Timothy Anderton, Brett H Andrews, Éric Aubourg, Stephen Bailey, Eduardo Balbinot, Rory Barnes, et al.  The ninth data release of the sloan digital sky survey: first spectroscopic data from the sdss-iii baryon oscillation spectroscopic survey.  The Astrophysical Journal Supplement Series, 203(2):21, 2012.
  * [6] T Padmanabhan.  Dark energy: mystery of the millennium.  In AIP Conference Proceedings, volume 861, pages 179–196. American Institute of Physics, 2006.
  * [7] N. Aghanim et al.  Planck 2018 results. VI. Cosmological parameters.  2018\.
  * [8] Shadab Alam, Metin Ata, Stephen Bailey, Florian Beutler, Dmitry Bizyaev, Jonathan A. Blazek, Adam S. Bolton, Joel R. Brownstein, Angela Burden, Chia-Hsun Chuang, and et al.  The clustering of galaxies in the completed sdss-iii baryon oscillation spectroscopic survey: cosmological analysis of the dr12 galaxy sample.  Monthly Notices of the Royal Astronomical Society, 470(3):2617–2652, Mar 2017.
  * [9] S. Joudaki et al.  KiDS+VIKING-450 and DES-Y1 combined: Cosmology with cosmic shear.  Astron. Astrophys., 638:L1, 2020.
  * [10] Adam G. Riess, Louise Breuval, Wenlong Yuan, Stefano Casertano, Lucas M. Macri, J. Bradley Bowers, Dan Scolnic, Tristan Cantat-Gaudin, Richard I. Anderson, and Mauricio Cruz Reyes.  Cluster cepheids with high precision gaia parallaxes, low zero-point uncertainties, and hubble space telescope photometry.  The Astrophysical Journal, 938(1):36, oct 2022.
  * [11] Adam G. Riess, Stefano Casertano, Wenlong Yuan, Lucas M. Macri, and Dan Scolnic.  Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM.  Astrophys. J., 876(1):85, 2019.
  * [12] Kenneth C. Wong et al.  H0LiCOW – XIII. A 2.4 percent measurement of H0 from lensed quasars: 5.3$\sigma$ tension between early- and late-Universe probes.  Mon. Not. Roy. Astron. Soc., 498(1):1420–1439, 2020.
  * [13] Wendy L. Freedman et al.  The Carnegie-Chicago Hubble Program. VIII. An Independent Determination of the Hubble Constant Based on the Tip of the Red Giant Branch.  Astrophys. J., 882:34, 2019.
  * [14] Luca Amendola and Shinji Tsujikawa.  Dark energy: theory and observations.  Cambridge University Press, 2010.
  * [15] Kazuharu Bamba, Salvatore Capozziello, Shin’ichi Nojiri, and Sergei D. Odintsov.  Dark energy cosmology: the equivalent description via different theoretical models and cosmography tests.  Astrophys. Space Sci., 342:155–228, 2012.
  * [16] Edmund J Copeland, Mohammad Sami, and Shinji Tsujikawa.  Dynamics of dark energy.  International Journal of Modern Physics D, 15(11):1753–1935, 2006\.
  * [17] P. J. E. Peebles and B. Ratra.  Quintessence: A review.  Reviews of Modern Physics, 75(2):559–606, 2003.
  * [18] C. Armendariz-Picon, V. Mukhanov, and P. J. Steinhardt.  k-essence as a model for dark energy.  Physical Review Letters, 85(15):4438–4441, 2001.
  * [19] Nandan Roy, Sangita Goswami, and Sudipta Das.  Quintessence or phantom: study of scalar field dark energy models through a general parametrization of the hubble parameter.  Physics of the Dark Universe, 36:101037, 2022.
  * [20] Aritra Banerjee, Haiying Cai, Lavinia Heisenberg, Eoin Ó. Colgáin, M. M. Sheikh-Jabbari, and Tao Yang.  Hubble sinks in the low-redshift swampland.  Phys. Rev. D, 103(8):L081305, 2021.
  * [21] Bum-Hoon Lee, Wonwoo Lee, Eoin Ó. Colgáin, M. M. Sheikh-Jabbari, and Somyadip Thakur.  Is local H 0 at odds with dark energy EFT?  JCAP, 04(04):004, 2022.
  * [22] C. Krishnan, E. Ó. Colgáin, M. M. Sheikh-Jabbari, and Tao Yang.  Running Hubble Tension and a H0 Diagnostic.  Phys. Rev. D, 103(10):103509, 2021.
  * [23] Rong-Gen Cai and Anzhong Wang.  Cosmology with interaction between phantom dark energy and dark matter and the coincidence problem.  JCAP, 03:002, 2005.
  * [24] G Mangano, Gennaro Miele, and V Pettorino.  Coupled quintessence and the coincidence problem.  Modern Physics Letters A, 18(12):831–842, 2003.
  * [25] H. Mohseni Sadjadi and M. Alimohammadi.  Cosmological coincidence problem in interactive dark energy models.  Phys. Rev. D, 74:103007, 2006.
  * [26] B. Wang, E. Abdalla, F. Atrio-Barandela, and D. Pavon.  Dark Matter and Dark Energy Interactions: Theoretical Challenges, Cosmological Implications and Observational Signatures.  Rept. Prog. Phys., 79(9):096901, 2016.
  * [27] J. F. Jesus, A. A. Escobal, D. Benndorf, and S. H. Pereira.  Can dark matter–dark energy interaction alleviate the cosmic coincidence problem?  Eur. Phys. J. C, 82(3):273, 2022.
  * [28] V. Salvatelli, A. Marchini, L. Pogosian, N. Vittorio, Y.-C. Wu, and J. Zavala.  Indications of a late-time interaction in the dark sector.  Physical Review Letters, 113(18):181301, 2014.
  * [29] A. Costa and P. G. Ferreira.  Hubble tension and interacting dark energy.  Journal of Cosmology and Astroparticle Physics, 2017(12):013, 2017\.
  * [30] E. Di Valentino, A. Melchiorri, and J. Silk.  Cosmological constraints from the combination of latest data sets: the role of dark energy interactions.  The European Physical Journal C, 79(2):139, 2019.
  * [31] S. Kumar, S. Kumar, K. Liao, and Y. Wang.  Interacting dark energy models with a logarithmic interaction term and their implications on the hubble tension.  Astrophysics and Space Science, 365(6):207, 2020.
  * [32] Eleonora Di Valentino, Alessandro Melchiorri, Olga Mena, and Sunny 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.
  * [33] Eleonora Di Valentino, Alessandro Melchiorri, Olga Mena, and Sunny Vagnozzi.  Nonminimal dark sector physics and cosmological tensions.  Phys. Rev. D, 101(6):063502, 2020.
  * [34] Weiqiang Yang, Supriya Pan, Eleonora Di Valentino, Rafael C. Nunes, Sunny Vagnozzi, and David F. Mota.  Tale of stable interacting dark energy, observational signatures, and the $H_{0}$ tension.  JCAP, 1809:019, 2018.
  * [35] Deng Wang.  The multi-feature universe: Large parameter space cosmology and the swampland.  Phys. Dark Univ., 28:100545, 2020.
  * [36] Luca Amendola.  Coupled quintessence.  Physical Review D, 62(4):043511, 2000.
  * [37] Glennys R Farrar and P James E Peebles.  Interacting dark matter and dark energy.  The Astrophysical Journal, 604(1):1, 2004.
  * [38] Nicola Tamanini.  Phenomenological models of dark energy interacting with dark matter.  Physical Review D, 92(4):043524, 2015.
  * [39] Luis P Chimento.  Linear and nonlinear interactions in the dark sector.  Physical Review D, 81(4):043525, 2010.
  * [40] Supriya Pan, Subhra Bhattacharya, and Subenoy Chakraborty.  An analytic model for interacting dark energy and its observational constraints.  Monthly Notices of the Royal Astronomical Society, 452(3):3038–3046, 2015.
  * [41] Valeria Pettorino, Carlo Baccigalupi, and Gianpiero Mangano.  Extended quintessence with an exponential coupling.  Journal of Cosmology and Astroparticle Physics, 2005(01):014, 2005\.
  * [42] Valeria Pettorino and Carlo Baccigalupi.  Coupled and extended quintessence: theoretical differences and structure formation.  Physical Review D, 77(10):103003, 2008.
  * [43] Wompherdeiki Khyllep, Jibitesh Dutta, Spyros Basilakos, and Emmanuel N. Saridakis.  Background evolution and growth of structures in interacting dark energy scenarios through dynamical system analysis.  Phys. Rev. D, 105(4):043511, 2022.
  * [44] Gabriela Caldera-Cabral, Roy Maartens, and L. Arturo Urena-Lopez.  Dynamics of interacting dark energy.  Phys. Rev. D, 79:063518, 2009.
  * [45] Luca Amendola.  Coupled quintessence.  Phys. Rev. D, 62:043511, 2000.
  * [46] Christian G. Boehmer, Gabriela Caldera-Cabral, Ruth Lazkoz, and Roy Maartens.  Dynamics of dark energy with a coupling to dark matter.  Phys. Rev. D, 78:023505, 2008.
  * [47] Hmar Zonunmawia, Wompherdeiki Khyllep, Nandan Roy, Jibitesh Dutta, and Nicola Tamanini.  Extended Phase Space Analysis of Interacting Dark Energy Models in Loop Quantum Cosmology.  Phys. Rev. D, 96(8):083527, 2017.
  * [48] Saddam Hussain, Saikat Chakraborty, Nandan Roy, and Kaushik Bhattacharya.  Dynamical systems analysis of tachyon-dark-energy models from a new perspective.  Phys. Rev. D, 107(6):063515, 2023.
  * [49] Sebastian Bahamonde, Christian G. Böhmer, Sante Carloni, Edmund J. Copeland, Wei Fang, and Nicola Tamanini.  Dynamical systems applied to cosmology: dark energy and modified gravity.  Phys. Rept., 775-777:1–122, 2018.
  * [50] Nandan Roy and Narayan Banerjee.  Generalized Brans-Dicke Theory: A Dynamical Systems Analysis.  Phys. Rev. D, 95(6):064048, 2017.
  * [51] Nandan Roy and Narayan Banerjee.  Dynamical systems study of Chameleon scalar field.  Annals Phys., 356:452–466, 2015.
  * [52] L. Arturo Urena-Lopez.  Scalar phantom energy as a cosmological dynamical system.  JCAP, 09:013, 2005.
  * [53] Burin Gumjudpai, Tapan Naskar, M. Sami, and Shinji Tsujikawa.  Coupled dark energy: Towards a general description of the dynamics.  JCAP, 06:007, 2005.
  * [54] Xi-ming Chen and Yungui Gong.  Fixed points in interacting dark energy models.  Phys. Lett. B, 675:9–13, 2009.
  * [55] Jose P. Mimoso, Ana Nunes, and Diego Pavon.  Asymptotic behavior of the warm inflation scenario with viscous pressure.  Phys. Rev. D, 73:023502, 2006.
  * [56] Nandan Roy.  Exploring the possibility of interacting quintessence model as an alternative to the $\Lambda$CDM model.  Gen. Rel. Grav., 55(10):115, 2023.
  * [57] Adam G Riess, Wenlong Yuan, Lucas M Macri, Dan Scolnic, Dillon Brout, Stefano Casertano, David O Jones, Yukei Murakami, Gagandeep S Anand, Louise Breuval, et al.  A comprehensive measurement of the local value of the hubble constant with 1 km $s^{-1}~{}mpc^{-1}$ uncertainty from the hubble space telescope and the sh0es team.  The Astrophysical journal letters, 934(1):L7, 2022.
  * [58] Luca Amendola.  Scaling solutions in general nonminimal coupling theories.  Phys. Rev. D, 60:043501, 1999.
  * [59] Tame Gonzalez, G. Leon, and I. Quiros.  Dynamics of quintessence models of dark energy with exponential coupling to dark matter.  Class. Quant. Grav., 23:3165–3179, 2006.
  * [60] Andrew P. Billyard and Alan A. Coley.  Interactions in scalar field cosmology.  Phys. Rev. D, 61:083503, 2000.
  * [61] Nandan Roy and Kazuharu Bamba.  Arbitrariness of potentials in interacting quintessence models.  Phys. Rev. D, 99(12):123520, 2019.
  * [62] Fabrizio F. Bernardi and Ricardo G. Landim.  Coupled quintessence and the impossibility of an interaction: a dynamical analysis study.  Eur. Phys. J. C, 77(5):290, 2017.
  * [63] Zhang Cong, Zhang Han, Yuan Shuo, Liu Siqi, Zhang Tong-Jie, and Sun Yan-Chun.  Four new observational h(z) data from luminous red galaxies in the sloan digital sky survey data release seven.  Research in Astronomy and Astrophysics, 14(10):1221, oct 2014.
  * [64] Joan Simon, Licia Verde, and Raul Jimenez.  Constraints on the redshift dependence of the dark energy potential.  Phys. Rev. D, 71:123001, Jun 2005.
  * [65] M Moresco, A Cimatti, R Jimenez, L Pozzetti, G Zamorani, M Bolzonella, J Dunlop, F Lamareille, M Mignoli, H Pearce, P Rosati, D Stern, L Verde, E Zucca, C.M Carollo, T Contini, J.-P Kneib, O. Le Fè vre, S.J Lilly, V Mainieri, A Renzini, M Scodeggio, I Balestra, R Gobat, R McLure, S Bardelli, A Bongiorno, K Caputi, O Cucciati, S. de la Torre, L. de Ravel, P Franzetti, B Garilli, A Iovino, P Kampczyk, C Knobel, K Kovač, J.-F. Le Borgne, V. Le Brun, C Maier, R Pelló, Y Peng, E Perez-Montero, V Presotto, J.D Silverman, M Tanaka, L.A.M Tasca, L Tresse, D Vergani, O Almaini, L Barnes, R Bordoloi, E Bradshaw, A Cappi, R Chuter, M Cirasuolo, G Coppa, C Diener, S Foucaud, W Hartley, M Kamionkowski, A.M Koekemoer, C López-Sanjuan, H.J McCracken, P Nair, P Oesch, A Stanford, and N Welikala.  Improved constraints on the expansion rate of the universe up to z 1.1 from the spectroscopic evolution of cosmic chronometers.  Journal of Cosmology and Astroparticle Physics, 2012(08):006–006, aug 2012.
  * [66] Antonio J. Cuesta et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations in the correlation function of LOWZ and CMASS galaxies in Data Release 12.  Mon. Not. Roy. Astron. Soc., 457(2):1770–1785, 2016.