Attractor Solutions in Interacting Dark Energy Models
${\bf Yousef~{}Bisabr}$***e-mail<EMAIL_ADDRESS>
Department of Physics, Shahid Rajaee Teacher Training University, Lavizan,
Tehran 16788, Iran

###### Abstract

We investigate a cosmological model in which dark energy, represented by a
quintessential scalar field, is coupled to a dark-matter perfect fluid in the
spatially flat Friedmann-Robertson-Walker Universe. This allows an energy
exchange in the dark sector which could happen both at early times before
recombination era or at late times. We use the coupling function
$Q=\gamma\rho_{dm}\dot{\varphi}$ which is induced by conformal transforming
scalar-tensor and $f(R)$ gravity theories to Einstein frame. It is argued that
there is a connection between this coupling function and $Q\propto\rho_{dm}H$.
A dynamical analysis is used to show that there are early- and late-time
attracting solutions for which the system evolves for a wide range of initial
conditions. These attractors generalize the scaling solutions which have been
already found in the non-interacting case.

##  1 Introduction

There are a large number of evidences which confirm that $\Lambda$CDM is a
successful model in describing cosmological observations. The model allows us
to describe the expansion history of the Universe from the beginning until
now. One important quantity in this description is $H_{0}$, the present value
of the Hubble parameter. There are two sets of measurements of $H_{0}$ which
are operative at high- and low-redshift patches. The early-time measurements
are based on CMB radiation or observations of baryonic acoustic oscillations
(BAO) which give $H_{0}=67.36\pm 0.54~{}kms^{-1}Mpc^{-1}$ [1]. The other set
is based on late-time (or local) measurements such as Cepheids-Calibrated
Supernovae which gives $H_{0}=74.03\pm 1.42~{}kms^{-1}Mpc^{-1}$ [2]. These two
values are discrepant at $>4\sigma$. In absence of systematic errors [3], new
physics or changing the expansion history seems to be necessary to remove the
tension between the two sets of observations.
Some extensions of $\Lambda$CDM (either at early times or late times) have
been proposed to resolve the problem [4]. An early-time proposal consists of
increasing the expansion rate around matter/radiation equality due to a new
energy component known as early dark energy (EDE) [5]. EDE behaves as a
cosmological constant before some critical redshift $z_{c}$ while after that
it dilutes faster than radiation to leave later expansion history unchanged.
Another proposal relies on possible interactions between dark energy (DE) and
dark matter (DM). In fact, different authors have already argued that due to
current lack of knowledge of the nature of DM and DE it is quite possible that
they interact and the interaction plays roles in expansion history of the
Universe. For instance, interacting dark energy (IDE) models can offer
explanations to the cosmic coincidence problem [7] and the cosmological
constant problem [8]. It is shown [6] that these interactions can also
alleviate the current tension on the value of the Hubble constant $H_{0}$.
Just like quintessence, dynamical dark energy which gives accelerating
expansion at late times, EDE is usually taken as a dynamical scalar field
which can generate accelerating expansion near the recombination era. In order
for EDE to have successful dynamics, it should arise without fine tuning and
imposing special initial conditions. To avoid fine tuning, it is necessary to
find attractor solutions in dynamics of EDE. If such attractor solutions exist
the system evolves for a wide range of initial conditions towards a unique
asymptotic behavior. In the present work, we show that such attracting
solutions actually exist when there is an interaction between DE and DM.

##  2 The Model

We will consider a two-component interacting system consisting of DM and a
real scalar field $\varphi$ as DE. Such an interaction is characterized by the
following conservation equations

$\nabla^{\mu}T_{\mu\nu}^{dm}=-\nabla^{\mu}T_{\mu\nu}^{\varphi}=Q_{\nu}$ (1)

where $T_{\mu\nu}^{dm}$ and $T_{\mu\nu}^{de}\equiv T_{\mu\nu}^{\varphi}$ are
energy-momentum tensors of DM and DE, respectively. The interaction current
$Q_{\nu}$ characterizes the energy transfer between the two energy components.
Assuming a flat Friedmann-Robertson-Walker metric
$ds^{2}=-dt^{2}+a^{2}(t)(\delta_{ij}dx^{i}dx^{j})$ with $"a"$ being the scale
factor, the energy balance equations (1) take the form

$\dot{\rho}_{\varphi}+3H(\omega_{\varphi}+1)\rho_{\varphi}=-Q$ (2)
$\dot{\rho}_{dm}+3H\rho_{dm}=Q$ (3)

where $\rho_{dm}$ is energy density of DM,
$\rho_{\varphi}\equiv\frac{1}{2}\dot{\varphi}^{2}+V(\varphi)$,
$p_{\varphi}\equiv\frac{1}{2}\dot{\varphi}^{2}-V(\varphi)$,
$\omega_{\varphi}\equiv p_{\varphi}/\rho_{\varphi}$ and $V(\varphi)$ is the
potential function. The Hubble parameter $H\equiv\frac{\dot{a}}{a}$ satisfies
the Friedmann equation†††Hereafter we use the units $\kappa^{2}\equiv 8\pi
G=1$ where $G$ is the gravitational constant.

$3H^{2}=\rho_{dm}+\rho_{\varphi}$ (4)

The exact form of the coupling function $Q$ is unknown and a wide variety of
theoretical and phenomenological interacting scenarios have been proposed and
investigated in the literature [9] [10] [11]. Here we consider the
consequences of the interaction function‡‡‡The covariant form of this function
is $Q_{\nu}=\gamma T^{dm}\nabla_{\nu}\varphi$ where $T^{dm}$ is the trace of
the energy-momentum tensor of dark matter. The trace $T^{dm}$ is solely given
by $\rho_{dm}$ for the pressureless dark matter. The trace dependence of
$Q_{\nu}$ means that any traceless energy component such as radiation does not
contribute to the interacting system (2) and (3). [12]

$Q=\gamma\rho_{dm}\dot{\varphi}$ (5)

This arises in different cases such as string theory, Einstein frame
representation of scalar-tensor and $f(R)$ gravity theories. Although $\gamma$
is generally an evolving function and depends on several variables, we will
take it as a constant parameter. This is a good approximation if the coupling
function $\gamma(\varphi,...)$ is slowly varying or the process of interaction
takes place in a short period of time with respect to the whole expansion
history of the Universe. The latter is particularly true for EDE scenario.
The equation (3) can be immediately solved to give

$\rho_{dm}=\rho_{0dm}a^{-3}e^{\gamma\varphi}$ (6)

with $\rho_{0dm}$ being an integration constant. This equation is equivalent
to $\rho_{dm}=\rho_{0dm}~{}a^{(-3+\varepsilon)}$ with
$\varepsilon=\frac{\gamma\varphi}{\ln a}$. While the coupling function $Q$
represents the interaction of DM and DE, $\varepsilon$ indicates how the
interaction affects the expansion rate of the Universe. Note that if
$\varepsilon$ is taken to be a constant parameter then
$\dot{\varphi}=\frac{\varepsilon}{\gamma}H$. In this case, (5) would be
equivalent to $Q=\varepsilon\rho_{dm}H$. This is the coupling function of
interacting models discussed by some authors [13].
In the case of null interaction ($Q=0$), DM energy density follows
$\rho_{dm}=\rho_{0dm}~{}a^{-3}$, the standard matter evolution in a matter-
dominated Universe. On the other hand, $Q\neq 0$ corresponds to an energy
exchange between DM and DE. In this case when $\varepsilon>0$ ($\gamma>0$),
energy is injecting into DM so that the latter will dilute more slowly
compared to the standard evolution. When $\varepsilon<0$ ($\gamma<0$) the
reverse is true, namely that DM is annihilated and the rate of its dilution is
faster than the case $\rho_{dm}\propto a^{-3}$.

##  3 Phase space analysis

To study dynamics of the equations (2), (3) and (4) one should write them as a
system of autonomous ordinary differential equations by defining a set of
appropriate variables. In general, there are many possible ways to define
these variables. We will use [14],

$X=\frac{1}{H}\frac{\dot{\varphi}}{\sqrt{6}}~{}~{}~{}~{},~{}~{}~{}Y=\frac{1}{H}\sqrt{\frac{V}{3}}$
(7)

We also consider an exponential potential $V(\varphi)=e^{-\lambda\varphi}$
[15] where $\lambda$ is a constant parameter. In terms of (7), the constraint
(4) takes the form

$1=\Omega_{dm}+\Omega_{\varphi}=\Omega_{dm}+X^{2}+Y^{2}$ (8)

where $\Omega_{dm}=\rho_{dm}/3H^{2}$ and
$\Omega_{\varphi}=\rho_{\varphi}/3H^{2}=X^{2}+Y^{2}$ are the density
parameters of DM and DE, respectively. Hence the two-dimensional phase space
is restricted to the region $0\leq X^{2}+Y^{2}\leq 1$ and $Y\geq 0$. This
means that the physical phase space of the model is contained within the unit
circle. The deceleration parameter and the effective equation of state of the
two interacting fluids are given by

$q=-1-\frac{\dot{H}}{H^{2}}=\frac{1}{2}+\frac{3}{2}(X^{2}-Y^{2})$ (9)
$\omega_{eff}=-1-\frac{2\dot{H}}{3H^{2}}=X^{2}-Y^{2}$ (10)

Now the set of equations (2), (3) and (4) with variables (7) gives rise to the
autonomous equations

$X^{\prime}=-3X+\sqrt{\frac{3}{2}}\lambda
Y^{2}-\sqrt{\frac{3}{2}}\gamma(1-X^{2}-Y^{2})-X\frac{H^{\prime}}{H}$ (11)
$Y^{\prime}=-\sqrt{\frac{3}{2}}\lambda XY-Y\frac{H^{\prime}}{H}$ (12)

where the prime indicates derivative with respect to $\ln a$. The expression
for $\frac{H^{\prime}}{H}$ can be straightforwardly obtained by taking
derivative of (4) and using (2) and (3). It gives

$\frac{H^{\prime}}{H}=-\frac{3}{2}(1+X^{2}-Y^{2})$ (13)

Combining the latter with (11) and (12) leads finally to

$X^{\prime}=f_{X}(X,Y)\equiv-3X+\sqrt{\frac{3}{2}}\lambda
Y^{2}-\sqrt{\frac{3}{2}}\gamma(1-X^{2}-Y^{2})+\frac{3X}{2}(1+X^{2}-Y^{2})$
(14) $Y^{\prime}=f_{Y}(X,Y)\equiv-\sqrt{\frac{3}{2}}\lambda
XY+\frac{3Y}{2}(1+X^{2}-Y^{2})$ (15)

The fixed points are those that satisfy $X^{\prime}=0$ and $Y^{\prime}=0$.
These points together with $\omega_{eff}$, $\Omega_{\varphi}$ and the
deceleration parameter are given in Tab.1 for the system (11) and (12).

Table 1: The fixed points and the corresponding effective equation of state for the system (14) and (15)  Name | X | Y | $\omega_{eff}$ | $\Omega_{\varphi}$ | $q$ | $q<0$
---|---|---|---|---|---|---
$A_{+}$ | $1$ | $0$ | $1$ | $1$ | $2$ | $no$
$A_{-}$ | $-1$ | $0$ | $1$ | $1$ | $2$ | $no$
B | $-\sqrt{\frac{2}{3}}\gamma$ | $0$ | $\frac{2}{3}\gamma^{2}$ | $\frac{2}{3}\gamma^{2}$ | $\frac{1}{2}+\gamma^{2}$ | $no$
$C_{+}$ | $\frac{\lambda}{\sqrt{6}}$ | $\sqrt{1-\frac{\lambda^{2}}{6}}$ | $-1+\frac{\lambda^{2}}{3}$ | $1$ | $-1+\frac{\lambda^{2}}{3}$ | $\lambda<\sqrt{2}$
$C_{-}$ | $\frac{\lambda}{\sqrt{6}}$ | $-\sqrt{1-\frac{\lambda^{2}}{6}}$ | $-1+\frac{\lambda^{2}}{3}$ | $1$ | $-1+\frac{\lambda^{2}}{3}$ | $\lambda<\sqrt{2}$
$D_{+}$ | $\frac{\sqrt{\frac{3}{2}}}{(\gamma+\lambda)}$ | $\sqrt{\frac{2\gamma(\gamma+\lambda)+3}{2(\gamma+\lambda)^{2}}}$ | $-\frac{\gamma}{(\gamma+\lambda)}$ | $\frac{\gamma(\gamma+\lambda)+3}{(\gamma+\lambda)^{2}}$ | $\frac{1}{2}-\frac{3\gamma}{2(\gamma+\lambda)}$ | $\frac{\gamma}{(\gamma+\lambda)}>\frac{1}{3}$
$D_{-}$ | $\frac{\sqrt{\frac{3}{2}}}{(\gamma+\lambda)}$ | $-\sqrt{\frac{2\gamma(\gamma+\lambda)+3}{2(\gamma+\lambda)^{2}}}$ | $-\frac{\gamma}{(\gamma+\lambda)}$ | $\frac{\gamma(\gamma+\lambda)+3}{(\gamma+\lambda)^{2}}$ | $\frac{1}{2}-\frac{3\gamma}{2(\gamma+\lambda)}$ | $\frac{\gamma}{(\gamma+\lambda)}>\frac{1}{3}$
Table 2: The fixed points and the corresponding effective equation of state for the system (14) and (15)  point | $\mu_{1}$ | $\mu_{2}$ | $stability$
---|---|---|---
$A_{+}$ | $3+2\sqrt{\frac{3}{2}}\gamma$ | $3-\sqrt{\frac{3}{2}}\lambda$ |  $\gamma<-\frac{3}{2}$, $\lambda>\sqrt{6}$
$A_{-}$ | $3-2\sqrt{\frac{3}{2}}\gamma$ | $3+\sqrt{\frac{3}{2}}\lambda$ |  $\gamma>\frac{3}{2}$, $\lambda<-\sqrt{6}$
B | $-\frac{3}{2}+\gamma^{2}$ | $\frac{3}{2}+\gamma(\lambda+\gamma)$ |  $\gamma<\pm\sqrt{\frac{3}{2}}$ , $\gamma(\lambda+\gamma)<-\frac{3}{2}$
$C_{\pm}$ | $\frac{\lambda}{2}(\lambda^{2}-6)$ | $-3+\lambda(\lambda+\gamma)$ |  $\\{\lambda<-\sqrt{6}\\}\bigcup\\{0<\lambda<\sqrt{6}\\}$, $\lambda(\lambda+\gamma)<3$
$D_{+}$ |  |  | $Fig.1$
$D_{-}$ |  |  | $Fig.1$

In order to study stability of the fixed points, one should consider
perturbations around them. In the linearized case, $\partial f_{i}/\partial
x_{j}$ is of particular importance. In our two-dimensional case, the
eigenvalues (denoted in the following by $\mu_{1}$ and $\mu_{2}$) of the
Jacobian matrix

$J=\left[{\begin{array}[]{cc}\partial f_{X}/X&\partial f_{X}/Y\\\ \partial
f_{Y}/X&\partial f_{Y}/Y\\\ \end{array}}\right]$

determine stability of the fixed points. The fixed points are stable if both
eigenvalues are negative, unstable if both are positive and saddle if one is
positive and the other is negative. For the dynamical system (11) and (12),
eigenvalues of $J$ and the stability status of the fixed points are given in
Tabs.2. The points $D_{\pm}$ have more complicated stability conditions. The
stability status of $D_{+}$ is depicted in Fig.1. The figure shows that
$D_{+}$ is stable ($\mu_{1,2}<0$) or saddle ($\mu_{1}<0,~{}\mu_{1}>0$ and vice
versa) in most regions of the parameters space.
The fixed points $A_{\pm}$ represent solutions for which the kinetic energy of
DE is dominated ($X=\pm 1$, $Y=0$). For these solutions, the right hand side
of (4) corresponds to an effective fluid with a stiff equation of state
$\omega_{eff}=1$. The point $B$ is generalization of the unstable solution
$x=0,~{}y=0$ of [14] in the non-interacting case $\gamma=0$. In our case, this
point is DM-dominated (for $\gamma<1$) which is affected by the interaction
and can now be stable for some values of $\gamma$ and $\lambda$ (see Tab. 2).
The stream plot of the system (14) and (15) is indicated in Fig. 2. In
Fig.2(a), $A_{\pm}$ are unstable nodes for which all the trajectories evolve
away from them and then are attracted towards the late-time attractor $C_{+}$.
In Fig.2(b),(c), all trajectories emerge from the early-time attractor $A_{-}$
before being attracted towards the saddle point $A_{+}$. These trajectories
ultimately end to the late-time attractors $B$ and $D_{+}$.

##  4 Conclusions

Dynamical dark energy models or quintessence are based on scalar fields with
appropriate potentials which can produce a late-time accelerating expansion.
The same scenario is recently used to alleviate the tension between local
determination of $H_{0}$ and the values inferred by CMB experiments at
recombination era. In fact, EDE is a scalar field which is responsible for the
required accelerating expansion at early times to alleviate the Hubble
tension. It behaves as a cosmological constant just before recombination and
then dilutes away like radiation or faster at late times. Due to the unknown
nature of the dark sector it is quite possible that there is an interaction
between DM and DE both at early and late times. At late times, it is shown
that interacting dark energy models play roles for addressing some
cosmological puzzles such as the coincidence and the cosmological constant
problems.
In this Letter, the main focus is on the problem of initial conditions. We
consider a cosmological model in which DM and DE are allowed to have mutual
interaction via the coupling function $Q=\gamma\rho_{dm}\dot{\varphi}$. The
advantage of this choice is that it arises from Einstein frame representation
of scalar-tensor theories and $f(R)$ gravity models. We first connect this
model to the class of models characterized by the coupling function
$Q\propto\rho_{dm}H$. This connection is established by taking $\varepsilon$,
denoting how the rate of dilution of $\rho_{dm}$ differs from
$\rho_{dm}\propto a^{-3}$, as a constant parameter. Then we have used a
dynamical analysis approach and an exponential potential for the DE scalar
field to show that there are early- and late-time attractor solutions in such
a dynamical system. This ensures that evolution of trajectories of the system
converges to a unique asymptotic behavior for a much wider range of initial
conditions. Moreover, we have shown that a DE-dominated ($\Omega_{\varphi}=1$)
solution exists which is a late-time attractor (Fig. 1(a)). We have also found
different late-time attractors where neither DE nor DM dominates the evolution
of the Universe.

## References

  * [1] N. Aghanim et al. (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), DOI:10.1051/0004-6361/201833910
Erratum: Astron.Astrophys. 652, C4 (2021)], DOI: 10.1051/0004-6361/201833910e

  * [2] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri and D. Scolnic, Astrophys. J. 876, 85 (2019), DOI:10.3847/1538-4357/ab1422
  * [3] G. Efstathiou, Mon. Not. Roy. Astron. Soc. 440, 1138 (2014), DOI: 10.1093/mnras/stu278
G. E. Addison et al., Astrophys. J. 818, 132 (2016), arXiv:1511.00055, DOI:
10.3847/0004-637X/818/2/132

  * [4] L. Knox and M. Millea, Phys. Rev. D101, 043533 (2020), DOI: 10.1103/PhysRevD.101.043533
  * [5] C. Wetterich, Phys. Lett. B 594, 17-22 (2004), DOI: 10.1016/j.physletb.2004.05.008
V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122,
221301 (2019), DOI: 10.1103/PhysRevLett.122.221301

  * [6] S. Kumar and R. C. Nunes, Phys. Rev. D94, 123511 (2016), DOI: 10.1103/PhysRevD.94.123511
E. Valentino, A. Melchiorri and O. Mena, Phys. Rev. D96, 043503 (2017), DOI:
10.1103/PhysRevD.96.043503

  * [7] D. Pavon and W. Zimdahl, Phys. Lett. B 628, 206 (2005), DOI: 10.1016/j.physletb.2005.08.134
S. del Campo, R. Herrera and D. Pavon, JCAP 0901, 020 (2009), DOI:
10.1088/1475-7516/2009/01/020
Y. Bisabr, Phys. Rev.D 82, 124041 (2010), DOI: 10.1103/PhysRevD.82.124041

  * [8] C. Wetterich, Astron. Astrophys. 301, 321 (1995), DOI: 10.48550/arXiv.hep-th/9408025
  * [9] Patrocinio Pérez, Ulises Nucamendi and Roberto De Arcia, Eur. Phys. J. C 81, 1063 (2021), DOI: 10.1140/epjc/s10052-021-09857-4
  * [10] M.B. Gavela et al., JCAP 0907, 034 (2009) (Erratum-ibid.1005:E01 (2010), (DOI: 10.1088/1475-7516/2010/05/E01), DOI:10.1088/1475-7516/2009/07/034
Bo-Yu Pu, Xiao-Dong Xu, Bin Wang and Elcio Abdalla, Phys. Rev. D 92, 123537
(2015), DOI: 10.1103/PhysRevD.92.123537
C. G. Boehmer, et al., Phys. Rev. D 78, 023505 (2008), DOI:
10.1103/PhysRevD.78.023505

  * [11] Ole Eggers Bjaelde, Subinoy Das, Adam Moss, JCAP 10, 017 (2012), DOI: 10.1088/1475-7516/2012/10/017
  * [12] Luca Amendola, Phys. Rev. D 62, 043511 (2000), DOI: 10.1103/PhysRevD.62.043511
  * [13] A. Paliathanasis, G. Leon, W. Khyllep, J. Dutta and S. Pan, Phys. J. C 81, 607 (2021), DOI: 10.1140/epjc/s10052-021-09362-8
G. Olivares, F. Atrio-Barandela and D. Pavon, Phys. Rev. D 77, 063513 (2008),
DOI: 10.1103/PhysRevD.77.063513
E. Valentino, A. Melchiorri and O. Mena, Phys. Rev. D 96, 043503 (2017), DOI:
10.1103/PhysRevD.96.043503

  * [14] E. J. Copeland, A. R. Liddle and D. Wands, Phys. Rev. D 57, 4686 (1998), DOI: 10.1103/PhysRevD.57.4686
  * [15] J. J. Halliwell, Phys. Lett. B185 341 (1987), DOI: 10.1016/0370-2693(87)91011-2
J. Barrow, Phys. Lett. B187 12 (1987), DOI: 10.1016/0370-2693(87)90063-3
B. Ratra, P. Peebles, Phys. Rev. D 37 3406 (1988), DOI:
10.1103/PhysRevD.37.3406

Figure 1: The parameters spaces corresponding to the eigenvalues (a) $\mu_{1}$
and (b) $\mu_{2}$ of the point $D_{+}$.

Figure 2: The phase plots of the dynamical system (14) and (15) for a)
$\lambda=1,~{}\gamma=0.1$ , b) $\lambda=4,~{}\gamma=-0.5$ and c)
$\lambda=10,~{}\gamma=0.5$ .