Phase Space Analysis of the Accelerating Multi-fluid Universe
S. D. Odintsov V. K. Oikonomou Petr V. Tretyakov
gr-qc astro-ph.CO hep-th





odintsov@ieec.uab.es ^1)Institut de
Ciencies de lEspai (IEEC-CSIC),  Carrer de Can Magrans, s/n, 08193
Barcelona, Spain 
 ^2)ICREA, Passeig LluAs Companys, 23,
    08010 Barcelona, Spain



v.k.oikonomou1979@gmail.com ^3)Laboratory
for Theoretical Cosmology, Tomsk State University of Control Systems
and Radioelectronics (TUSUR), 634050 Tomsk, Russia
 ^4)Tomsk
State Pedagogical University, 634061 Tomsk, Russia




tpv@theor.jinr.ru ^5)Joint Institute for
Nuclear Research, Joliot-Curie 6, 141980 Dubna, Moscow region,
Russia
 ^6)Institute of Physics, Kazan Federal University,
Kremlevskaya street 18, 420008 Kazan, Russia


We study in detail the phase space of a Friedmann-Robertson-Walker
Universe filled with various cosmological fluids which may or may
not interact. We use various expressions for the equation of state,
and we analyze the physical significance of the resulting fixed
points. In addition we discuss the effects of the stability or an
instability of some fixed points. Moreover we study an interesting
phenomenological scenario for which there is an oscillating
interaction between the dark energy and dark matter fluid. As we
demonstrate, in the context of the model we use, at early times the
interaction is negligible and it starts to grow as the cosmic time
approaches the late-time era. Also the cosmological dynamical system
is split into two distinct dynamical systems which have two distinct
de Sitter fixed points, with the early-time de Sitter point being
unstable. This framework gives an explicit example of the
unification of the early-time with late-time acceleration. Finally,
we discuss in some detail the physical interpretation of the various
models we present in this work.


04.50.Kd, 95.36.+x, 98.80.-k, 98.80.Cq, 11.25.-w 

Phase Space Analysis of the Accelerating Multi-fluid Universe 
    Petr V. Tretyakov^5,6
    December 30, 2023
==============================================================






§ INTRODUCTION



In modern theoretical cosmology, the most striking event was the
observation of the late-time acceleration <cit.> that our
Universe undergoes at present time. Admittedly this observation has
utterly changed the way of thinking of modern cosmologists, since
this late-time acceleration is a feature of our Universe that was
never thought it would actually occur. Consequently, the focus for
the last nearly 20 years is to model in a successful way this
late-time acceleration and also to harbor the late and early-time
acceleration era in a unified theoretical framework. Towards this
unified description, many proposals, especially those suggesting to
modify the gravitational sector, have been introduced ever since,
see the reviews <cit.> for
details. From the first moment that the late-time acceleration has
been observed, it was realized that no perfect matter fluid known at
that time was able to realize the late-time acceleration era, and
therefore the need for alternative generalized cosmological fluids
was compelling. By using generalized cosmological fluids, both the
late and early-time acceleration era can be realized, and up to date
there are many theoretical proposals that use generalized fluids,
for example in Refs.
<cit.>
imperfect fluids are used in order to describe the cosmological
evolution of our Universe, and in some cases certain particular
examples are used, called viscous fluids are (see
<cit.> for reviews). It is notable that the imperfect
fluids may describe even phantom evolution of our Universe, without
using phantom scalar fields, which violate the energy conditions,
see Refs. <cit.> for details. Furthermore,
other cosmological evolutions like bouncing cosmology, in the
context of both classical and loop quantum cosmology imperfect
fluids were studied in <cit.>, and also singular
cosmology can be realized by imperfect fluids
<cit.>. Furthermore, several models
which take into account bulk viscosity were discussed in Refs.
<cit.>, and also
an important class of models which assume an interaction between
dark matter and dark energy fluids, can be found in Refs.
<cit.>.

In this paper we shall perform a detailed phase-space analysis of a
Friedmann-Robertson-Walker Universe, filled with different fluid
components, which may or may not interact between them. The
dynamical evolution of such kind of model is described by the
Friedmann equations,

    3/κ^2H^2=∑ρ_i,


    -1/κ^2(2Ḣ + 3H^2)=∑ p_i,

or equivalently,

    -2/κ^2Ḣ = ∑(ρ_i+ p_i),

and also the energy conservation equations hold true,

    ρ̇_i=-3H(ρ_i+p_i).

In the above, κ^2=8π G and also the equation of state (EoS)
p_i=p_i(ρ_i) may be highly non-trivial for some fluid
components. We shall appropriately choose the variables in order to
capture the phase space dynamics in the most optimal way, and we
shall analyze the structure of the phase space by providing an
analytic treatment of the cosmological dynamical equations. Seeing
the cosmological equations as a dynamical system is a particularly
appealing way to investigate the phase space in many cosmological
contexts, see for example
<cit.>,
but also in modified gravity too, see for example
<cit.>.
In most cases, the choice of the dynamical system variables plays a
crucial choice, and in some cases the resulting cosmological
dynamical system may be rendered autonomous <cit.>.
Also it is possible to choose dimensionless variables, see for
example Refs. <cit.> for an F(R)
gravity cosmological dynamical system, and also see Ref. <cit.>,
for a cosmological theory with higher derivatives of the scalar
curvature. The dynamical systems approach for cosmological systems
has many attributes, with the most important being the fact that the
fixed points of the dynamical system actually provide new insights
with regards to the behavior of the attractor solutions and also
reveals the stability structure of the dynamical system near the
attractors. It is conceivable that the choice of the variables plays
an important role, as we also demonstrate by this work.


This paper is organized as follows: In section <ref> we
present some well known features of dynamical systems approach in
cosmological context with generalized fluids, in section <ref>
we discuss in brief how interactions between dark matter and dark
energy may be introduced and by using appropriately chosen variables
we present how the dynamical systems analysis can be performed in
this case. In section <ref> we generalize the formalism we
developed in the previous sections and by using dimensionless
variables, we investigate the physical consequences of having
various equations of state for the fluid components of the
cosmological system. In section <ref> by using appropriately
chosen dimensionless variables, we investigate in detail how the
interaction of dark matter and dark energy fluids may affects the
phase space structure. We study the stability and behavior of the
fixed points of the dynamical system and we also discuss how the
early and late-time acceleration eras are affected by the various
functional forms of the interaction coupling between dark matter and
dark energy. Finally the conclusions follow in the end of the paper.







§ STANDARD APPROACH ON DYNAMICAL SYSTEMS AND COSMOLOGICAL DYNAMICS


In this section we present the simplest case of the dynamical
systems approach in cosmological dynamics. We consider the simplest
case in which the Universe is filled with radiation ρ_r and a
perfect fluid with a non-trivial EoS of the form
p=-ρ+f(ρ)+G(H) <cit.>. Such non-trivial EoS can be
considered as some sort of viscous fluid or generalized EoS fluid,
see <cit.> for a review on this topic. In addition,
such an EoS can be considered as an effective fluid presentation of
some modified gravity theory <cit.>. In the case
at hand, the full dynamical system takes the following form,

    3/κ^2H^2      =ρ_r+ρ=ρ_tot,
    ρ̇_r      =-4Hρ_r,
    ρ̇      =-3H[f(ρ)+G(H)].

The appearance of the term G(H) might seem unconventional, from a
thermodynamic point of view, and we need to briefly describe the
motivation for using such a term. This term encompasses the viscous
part of the cosmological fluid, so it mainly quantifies the
viscosity of the fluid. In the Universe, nd especially in the very
early stages of it's evolution, the effects of a viscous
cosmological component are most likely expected to occur during the
neutrino decoupling process, which occurs at the end of the lepton
era <cit.>. Hence, viscosity is encompassed in the
very own fabric of the Universe. In addition, a strong motivation
for using viscous fluid components comes from the fact that the
perfect fluid approach among cosmologists-hydrodynamicists is just
an ideal approach, and does not describe the real world. Finally,
due to the fact that early and late-time acceleration may be
described by an unknown form of a cosmological fluid, it is natural
to assume that the fluid has the most general form, which means that
a viscous component is needed[Note that terms of the form
G(H,Ḣ) in the cosmological relation between the effective
pressure and the energy density, naturally occur in modified gravity
thermodynamics <cit.>].

Having described the motivation for the use of viscous fluid
components, we can rewrite the above cosmological equations in terms
of dimensionless variables. In the case at hand, there is only one
independent variable due to the constraint equation (<ref>). By
using the e-foldings number[In this case time derivatives
for some variable X
  transform as X'=dX/dN=Ẋ/H] N=ln a and by introducing a new dimensionless
  variable defined as,

    x=κ^2/3H^2ρ=ρ/ρ_tot,

we obtain the next dynamical equation,

    x'=3/ρ_tot [ f(ρ) + G(H) ](x-1) -4x^2 +4x,

where the variables ρ, H and ρ_tot must be expressed
in terms of x depending on the choices of the functions f and
G.


Let us here discuss the simplest choice of EoS, which is p=w_0ρ
+w_1H^2, which implies that f(ρ)=ρ(1+w_0) and
G(H)=w_1H^2. It is easy to see that in this case, Eq. (<ref>)
takes the following form,

    x'= [ 3(1+w_0)x +κ^2w_1 ](x-1) -4x^2 +4x≡ m(x),

where the “prime” denotes differentiation with respect to the
e-foldings number. The equation that determines the stationary
points for the above dynamical equation, is quadratic with
discriminant 𝔇=(κ^2w_1 -4 + 3(1+w_0)
)^2, so the existing solutions are always real. For the case
𝔇=0 there is the only one solution, which is,
x=1/2-κ^2w_1/2(3(1+w_0)-4). For the case
𝔇>0 there are two solutions, which are:

    x_1      =1,
    
    x_2      =-κ^2w_1/3(1+w_0)-4,

and this case is the most interesting, since it provides us plenty
dynamical solutions. Now note that the physical values that the
variable x can take, are 0⩽ x⩽ 1, but for x_2
we have 0<x_2<1, which give us the following restrictions for the
free parameters:

    w_1<0,  w_0>-1/4,  -κ^2w_1<-4+3(1+w_0),
          w_1>0,  w_0<-1/4,  -κ^2w_1>-4+3(1+w_0).

Let us study the stability conditions of the existing stationary
points. It is clear that for every point x_i there is only one
eigenvalue μ_i=m'(x=x_i), and for the fixed points at hand, we
have,

    μ_1(x_1)      =-4+3(1+w_0)+κ^2w_1,
    μ_2(x_2)      =4-3(1+w_0)-κ^2w_1=-μ_1 .

Therefore, we find that for any values of parameters w_0, w_1
one of the fixed points is stable and the other one is unstable.
Note that for the case 𝔇=0 we have μ=0 and
therefore it is compelling to investigate the corresponding center
manifold, but fortunately this is not interesting case from a
physical point of view. The physical significance of the stationary
point x_1 is that it corresponds to a Universe with
ρ_tot=ρ and ρ_r=0. With regard to the fixed point
x_2, it corresponds to a Universe with some fixed relation between
ρ and ρ_r.





§ MODELS WITH DARK MATTER INTERACTING TO DARK ENERGY




Let us discuss at this point some models which describe interactions
of dark matter with dark energy. The dynamical system in this case
can be written in the following form,

    3/κ^2H^2      =ρ_r+ρ+ρ_dm=ρ_tot,
    ρ_r'      =-4ρ_r,
    ρ'      =-3[f(ρ)+G(H)]-Q/H,
    ρ'_dm      =-3ρ_dm+Q/H-3(-3Hζ),

where the prime denotes as previously differentiation with respect
to the e-foldings number N, and Q quantifies the interaction
between dark energy and dark matter, and the term -3Hζ
corresponds to the bulk viscous pressure of the dark matter fluid.

We will assume that the bulk viscous coefficient ζ has the
following form <cit.>,

    ζ=ζ_0/√(3κ^2)ρ_tot^1/2=1/κ^2Hζ_0.

The motivation for using this kind of ansatz, comes from
astronomical estimates on the relation between ζ and
ρ_tot, which was studied in Ref. <cit.>,
where a general form of the ζ-ρ_tot relation was assumed,
and it was of the form ζ∼ρ^λ_tot. The choice
ζ∼ρ^1/2_tot we used in Eq. (<ref>) yields a
good fit between astronomical constraints and the fluid approach.

In this case we introduce the following set of dimensionless
variables,

    x      =κ^2/3H^2ρ=ρ/ρ_tot,
    
    y      =κ^2/3H^2ρ_dm=ρ_dm/ρ_tot,
    
    q      =κ^2/3H^3Q=Q/Hρ_tot,

where only x and y are independent dynamical variables. By using
the set of dimensionless variables, the dynamical equations
(<ref>)-(<ref>) take the form[Note that the
constraint (<ref>) and the equation for ρ_r are already
took into account in this system.]:

    x'      =-4x^2+4x-yx-3xζ_0-q+3(f+G)/ρ_tot(x-1),
    
    y'      =-y^2-4xy+y(1-3ζ_0)+3ζ_0+q+3(f+G)/ρ_toty,

and by using the EoS we used in the previous section, namely,
f(ρ)=ρ(1+w_0), G(H)=w_1H^2 we obtain,





    x'      =x^2(3w_0-1)+x(w_1κ^2-3ζ_0-3w_0+1)-xy-(q+w_1κ^2),
    
    y'      =-y^2+xy(3w_0-1)+y(w_1κ^2+1-3ζ_0)+3ζ_0+q.

We can easily find the stationary points for the above dynamical
system, by multiplying (<ref>) by y, (<ref>) by (1-x)
and by combining the resulting equations we find,

    (1-x_0-y_0)(y+3ζ_0+q)=0.


The above equation indicates that there exist at most three
stationary points, and it is clear that even in the most general
case, the stationary points may be found analytically, and these are
equal to,

    x_0=q+w_1κ^2/1-3w_0,  
    y_0=-3ζ_0-q,


    x_0=1-y_0,  
    y_0=1/6w_0 [ 3w_0+w_1κ^2-3ζ_0 ±√((3w_0+w_1κ^2-3ζ_0)^2+12w_0q+36w_0ζ_0) ].

Since there are a many free parameters, we investigate some
particular cases, which are interesting from a physical point of
view. We start off with the case w=0, in which case the first
fixed point takes the form,

    y_0=-3ζ_0-q,   x_0=q+w_1κ^2,

and the corresponding eigenvalues are,

    μ_1      =1,
    μ_2      =1-w_1κ^2+3ζ_0,

Accordingly, the second fixed point is,

    y_0=3ζ_0+q/3ζ_0-w_1κ^2,   x_0=1-y_0,

and the corresponding eigenvalues are,

    μ_1      =w_1κ^2-3ζ_0,
    μ_2      =w_1κ^2-3ζ_0-1,

so the first point always unstable, whereas the second one may be
stable, depending on the choice of the free parameters. However, if
the first point lies in the physically allowed region, the second
fixed point is rendered always unstable. Typical phase portrait with
the two fixed points in the physical region, are presented in
Fig.<ref>.

Now let us consider another physically interesting case, for which
q=-3ζ_0. In this case, there exist three fixed points which
are,

    y_0      =0,  x_0=3ζ_0-w_1κ^2/3w_0,  μ_1=w_1κ^2-3ζ_0+1+x_0(3w_0-1),  μ_2=2μ_1-w_1κ^2+3ζ_0-3w_0-1,
    
    y_0      =0,  x_0=1,  μ_1=w_1κ^2-3ζ_0+3w_0,  μ_2=w_1κ^2-3ζ_0+3w_0-1,
    
    y_0      =3w_0+w_1κ^2-3ζ_0/3w_0,  x_0=1-y_0,  μ_1=-1,  μ_2=  -3w_0-w_1κ^2+3ζ_0.

A particularly interesting subcase of the above is if we further
choose, 3ζ_0=w_1κ^2, then the fixed points
become,

    x_0      =0,  y_0=1,  μ_1,2=-1,  -3w_0,
    
    x_0      =1,  y_0=0,  μ_1,2=3w_0,  3w_0-1,
    
    x_0      =0,  y_0=0,  μ_1,2=1,  1-3w_0.

We can see that two fixed points are always unstable, whereas one of
the three, the first or the second one depending on the sign of the
parameter w_0, is stable.





In Fig. <ref> we plotted the phase portrait for the case
w_0=0.5, while in Fig. <ref> we plotted the phase
portrait for w_0=-0.5, and finally in Fig. <ref> we can
see the phase portrait for w_0=0. As it can be seen in all
figures, two of the three fixed points are unstable and one of the
three is stable. Furthermore, it can be seen that by using Eq.
(<ref>), it is possible to find the case (by appropriately
choosing the parameters) for which there will be some stationary
point with fixed relation x/y≠0,1. This task may be solved
numerically but we refrain from going into details on this.








§ GENERALIZED FORM OF COSMOLOGICAL FLUIDS



In this section we extend the cases we presented in section
<ref> to include generalized form of the EoS. We consider the
simplest scenario for which the Universe is filled with ρ_r and
a perfect fluid with non-trivial EoS p=-ρ+f(ρ)+G(H)
<cit.>. In this case, the dynamical system takes the following
form,

    3/κ^2H^2      =ρ_r+ρ=ρ_tot,
    -2/κ^2Ḣ      = 4/3ρ_r + [f(ρ)+G(H)],
    ρ̇_r      =-4Hρ_r,
    ρ̇      =-3H[f(ρ)+G(H)].

By using the e-foldings number as independent variable and also by
introducing the dimensionless variables,

    x=κ^2/3H^2ρ=ρ/ρ_tot,


    z=1/κ^2H^2=3/ρ_totκ^4

the dynamical system can be cast in the following form,

    x'      =κ^4z [ f(ρ) + G(H) ](x-1) -4x^2 +4x,
    
    z'      =4z(1-x)+z^2κ^4 [ f(ρ) + G(H) ],

where according to the definition, ρ=3x/κ^4z and
H^2=1/κ^2z. Moreover we can see from the system
(<ref>)-(<ref>) that the fixed point x_0=1,  z_0=0
always exists, except for some very special choices of functions f
and G.


It is interesting to note that actually due to the constraint
equation (<ref>), in the case at hand, there is only one
independent variable, but by introducing the variable z and by
taking into account the equation (<ref>), allows us to obtain
more information about the evolution of the dynamical system. For
instance, one of the most hard tasks in such kind of dynamical
systems analysis, is to interpret correctly the physical meaning of
the stationary points. This approach was firstly proposed in Ref.
<cit.>. Let us introduce an additional parameter, which will be
very helpful for this interpretation, which is the effective
equation of state, which we denote as w_eff, and it is defined
as follows:

    w_eff≡ -1-2Ḣ/3 H^2.

The parameter w_eff can be expressed in terms of the
dimensionless variables (<ref>)-(<ref>) as follows,
 
    w_eff=1/3-4/3x+1/ρ_tot(f+G).

By specifying the EoS, it is possible to obtain various physically
interesting evolution scenarios, so in the rest of this section we
shall specify the EoS and we study in detail the dynamical evolution
stemming from the choice of EoS.





 §.§ A Simplified Form of EoS


Let us discuss the simplest case of EoS, which is p=w_0ρ
+w_1H^2, which in turn implies f(ρ)=ρ(1+w_0),
G(H)=w_1H^2. It is easy to verify that in this case, the equations
(<ref>)-(<ref>) take the following form,

    x'      = x^2(3w_0-1) +x(1-3w_0+w_1κ^2)-w_1κ^2,
    
    z'      = z[ x(3w_0-1) +4+w_1κ^2].

Thus we have the following stationary points for the dynamical
system above,


  * First fixed point  x_0=w_1κ^2/1-3w_0,   z_0=0.

  * Second fixed point x_0=1,   z_0=0.


For the first fixed point we need to note that, this point lies in
the physical region only if 0<w_1κ^2/1-3w_0<1.
The corresponding eigenvalues are μ_1=1-3w_0-w_1κ^2, and
μ_2=4. Using Eq. (<ref>) we find that at this point we
have w_eff=1/3.

With regard to the second fixed point,
the eigenvalues are μ_1=3w_0+w_1κ^2-1, and μ_2=μ_1+4.
Correspondingly we find that at this point
w_eff=w_0+1/3w_1κ^2. Thus we can see that the first
fixed point always unstable and also that the second fixed point may
be stable only if one (or both) of parameters w_0, w_1 are
strictly negative. For instance if we require that the first point
lies in the physical region, we find that the second point is stable
for w_0<-1. In Figs.  (<ref>-<ref>) we plot the
typical behavior of the phase trajectories.

Particularly, Fig. <ref> corresponds to w_0=-2, and
w_1κ^2=0.5. The left fixed point corresponds to the effective
EoS parameter w_eff=1/3 and for the right fixed point,
we have w_eff=-11/6. Clearly the left fixed point
represents radiation, while the right one corresponds to some
phantom evolution. In Fig. <ref>, the phase portrait
corresponds to the following choices for the parameters, w_0=0.1,
w_1κ^2=0.3. The left point corresponds to an effective EoS
parameter w_eff=1/3, while the right point corresponds
to w_eff=0.2, which describes a form of collisional matter
<cit.>.




 §.§ More Complicated Forms of the EoS


Now let us study more complicated forms of the EoS, and we choose it
to be of the form f(ρ)+G(H)=Aρ^α+BH^2β. This
EoS is known to lead the cosmological system to finite-time
singularities, as this was demonstrated in Ref. <cit.>. For this
EoS, the equations (<ref>)-(<ref>) take the form,

    x'      = -4x^2+4x+κ^4z[ A(3x/κ^4z)^α+B(1/κ^2z)^β](x-1),
    
    z'      =4z(1-x)+κ^4z^2[ A(3x/κ^4z)^α+B(1/κ^2z)^β].

Correspondingly, the effective EoS in this case reads,
 
    w_eff=1/3 [1-4x + Aκ^4z(3x/κ^4 z)^α+Bκ^4 z1/(κ^2 z)^β ].




It is clear that in the most general case, this system may be solved
only numerically, so let study some appropriately chosen cases which
admit analytical solutions. Consider first the case for which
α=1, β=2, in which case, the dynamical system
(<ref>)-(<ref>) takes the following form,

    x'      = ( 3Ax-4x+B/z)(x-1),
    
    z'      =4z(1-x)+3Axz+B.

In this case there are two stationary points: the first is x_0=1,
z_0=0[Note here that in all these three cases existence
of the point x_0=1, z_0=0 is not an obvious solution, but it can
be confirmed by numerical investigations.] and the second is
x_0=1, z_0=-B/3A. For the first point, we have
w_eff=sign(B)∞[Infinite values of w_eff looks
like the Ruzmaikin solution at t→ 0.] and for the second
one we have w_eff=-1.

We can see that the second point corresponds to some new non-trivial
de Sitter state with ρ_tot=ρ and H^2=H_0^2≠ 0. The
eigenvalues of the second fixed point are μ_1=-4, and
μ_2=3A. The typical behavior of the phase trajectories
corresponding to this case, can be found in Fig.<ref>, for
A=-1, and B=1. As it can be seen, there exist trajectories which
start from the first fixed point  x_0=1, z_0=0 (note that zero
values of z correspond to infinite values of H, which indicates
a singularity) and end up to the second fixed point, with
non-singular and non-zero values of H. In effect, the second fixed
point may be interpreted as a late-time acceleration de Sitter point
of the cosmological dynamical system.


Now let us consider the case for which α=2, and β=1, in
which case, the dynamical system (<ref>)-(<ref>) takes the
form,

    x'      = ( 9A/κ^4x^2/z-4x+κ^2 B)(x-1),
    
    z'      =4z(1-x)+9A/κ^4x^2+κ^2 B z.

In this case there are two stationary points, namely, x_0=1,
z_0=0 and also x_0=1, z_0=-9A/κ^6B and the
situation is very similar to the previous case. For the first fixed
point we have w_eff=sign(A)∞ and for the second one
w_eff=-1. The eigenvalues of the second point are μ_1=-4,
and μ_2=κ^2B. Note also that for both these cases, A and
B must have opposite sign, in order for the second fixed point to
lie in the physical region. We also need to note that for this case
there are additional fixed points, which are difficult to find
analytically, but the most physically interesting cases of fixed
points are the ones we just presented. The phase space behavior
corresponding to this case can be found in Fig. <ref>, for
A=1/2, B=-9, κ^2=1, and as it can be seen, the
behavior of the trajectories is similar to the previous case.

Concluding this section, let us briefly discuss another choice of
parameters for which it is possible to find analytically the fixed
points, and this occurs for the choice α=2, β=2, in
which case, the dynamical system (<ref>)-(<ref>) takes the
form,

    x'      = ( 9A/κ^4x^2/z-4x+B/z)(x-1),
    
    z'      =4z(1-x)+9A/κ^4x^2+B.

In this case, the stationary points are the following, x_0=1,
z=0 and x_0^2=-Bκ^4/9A, z=0. For the first fixed
point we have w_eff=sign(Bκ^4+9A)∞ and for the second
one w_eff=1/3-4/9κ^2√(-B/A).
The behavior of the trajectories corresponding to this case can be
found in Fig. <ref>, for the choice A=-1, B=1, and in
Fig. <ref>, for the choice A=1, B=-1.






§ DARK MATTER INTERACTING WITH DARK ENERGY: SOME NON-TRIVIAL MODELS




Let us now discuss some non-trivial models which describe
interactions between the dark matter and dark energy fluids. The
dynamical system in this case may be written in the following form,

    3/κ^2H^2      =ρ_r+ρ+ρ_dm=ρ_tot,
    -2/κ^2Ḣ      = 4/3ρ_r + [f(ρ)+G(H)]+ρ_dm+(-3Hζ),
    ρ̇_r      =-4Hρ_r,
    ρ̇      =-3H[f(ρ)+G(H)+Q/3HH^2k],
    ρ̇_dm      =-3H[ρ_dm-Q/3HH^2k+(-3Hζ)],

where we modified the interaction between dark energy and dark
matter by using the multiplier H^2k[Note that case
k=0 corresponds to the standard interaction we presented in a
previous section.], and also the term -3Hζ, which corresponds
to the bulk viscous pressure of the dark matter fluid. We will
assume that the bulk viscous coefficient ζ has the following
form,

    ζ=ζ_0/√(3κ^2)ρ_tot^1/2=1/κ^2Hζ_0.

We define the set of dimensionless variables as follows,

    x      =κ^2/3H^2ρ=ρ/ρ_tot,
    
    y      =κ^2/3H^2ρ_dm=ρ_dm/ρ_tot,
    
    z      =1/κ^2H^2=3/ρ_totκ^4,
    
    q      =κ^2/3H^3Q=Q/Hρ_tot,

where only x and y are dynamical independent variables. In the
new variables system, Eqs. (<ref>)-(<ref>) take the
form[Note that the constraint (<ref>) and the equation
for ρ_r have already been taken into account for this system.]:

    x'      =-4x^2+4x-yx-3xζ_0-q/(κ^2z)^k+κ^4z(f+G)(x-1),
    
    y'      =-y^2-4xy+y(1-3ζ_0)+3ζ_0+q/(κ^2z)^k+κ^4zy(f+G),
    
    z'      =4z(1-x-y)+z^2κ^4(f+G)+3yz-3zζ_0,

and using that the EoS has the form we used in the previous section,
namely, f(ρ)=ρ(1+w_0), G(H)=w_1H^2, we obtain,

    x'      =(x-1)(w_1κ^2+3xw_0-x)-x(y+3ζ_0)-q/(κ^2z)^k,
    
    y'      =-y^2+y(1-3ζ_0+3xw_0-x+w_1κ^2)+3ζ_0+q/(κ^2z)^k,
    
    z'      =z(4-x-y+3xw_0-3ζ_0+w_1κ^2).

The effective EoS in this case (for arbitrary functions f and G)
reads,
 
    w_eff=1/3 [1-4x -y + κ^4z(f+G) -3ζ_0  ],

which in the case that f(ρ)=ρ(1+w_0), G(H)=w_1H^2,
becomes,
 
    w_eff=1/3 [1-x -y + 3xw_0 +w_1κ^2 -3ζ_0  ].

Now we consider some special cases, and we start with the case
k=0. This case corresponds to the usual interaction term. Since
the first two  equations are identical to (<ref>)-(<ref>),
we have the solutions (<ref>), (<ref>), where we
need to add z_0=0.

Also, since the first two equations do not depend on z, the first
and second eigenvalues will be totally identical to the ones we
obtained in section <ref>, and one additional eigenvalue
appears for every stationary point, which is,

    μ_3=4-x_0-y_0+3x_0w_0-3ζ_0+w_1κ^2.

Also note that the stability or instability of the point z_0=0
with respect to this additional coordinate z, implies stability in
the past or in the future correspondingly. It mean that a stable
with respect to coordinates x, y point which have z_0=0 will
be stable in the past if it is stable with respect to coordinate
z, or equivalently will be stable for infinite values of H. In
addition, it will stable in the future, if it is unstable with
respect to coordinate z, or equivalently it will be stable for
small (zero) values of H.

The effective EoS (<ref>) for the fixed point
(<ref>) is equal to w_eff=1/3 for any values of
parameters, which clearly describes a radiation dominated era.
However, the effective EoS for the fixed points (<ref>),
has a more complicate structure, which is given below,

    w_eff=1/2w_0+1/6w_1κ^2-1/2ζ_0 ±1/6√((3w_0+w_1κ^2-3ζ_0)^2+12w_0q+36w_0ζ_0).

Let us now further analyze the case at hand, by specifying the
values of the free parameters, so we start with the parameter w_0,
and assume for the moment that w_0=0. In this case the fixed point
(<ref>) has an additional eigenvalue μ_3=4 and the
corresponding effective EoS becomes w_eff=1/3, which
describes radiation. Moreover, the fixed point (<ref>) has the
additional eigenvalue μ_3=3-3ζ_0+w_1κ^2 and the
corresponding effective EoS for this point is
w_eff=1/3w_1κ^2-ζ_0. Thus in Fig.
<ref> the left fixed point corresponds to radiation with
w_eff=1/3, and the right point has
w_eff=11/3, and for both the fixed points, H has
infinite values. In Table <ref> we have gathered all the fixed
points which correspond to the case k=0, w_0=0.

As it can be seen in Table <ref>, the fixed point
P_1b may describe late-time acceleration. Indeed, if we
set w_eff=-1+α, with 0<α≪ 1, we obtain
μ_1=-3+3α, μ_2=-4+3α and μ_3=α and as we
already noted, this means that this point is stable in the future
(for small values of H). Moreover, by changing the parameter q,
we can provide any interesting relation between ρ_dm and
ρ. For instance if put q=3/4-3ζ_0 we obtain for
this point y_0≡ρ_dm/ρ_tot=1/4.



Let us discuss some alternative choices for the parameters, so
consider the case q=-3ζ_0, in which case the fixed point
(<ref>) has an additional eigenvalue, which we denote
μ_3, and it is equal to,
μ_3=4-3ζ_0-w_1κ^2/3w_0 and the corresponding
effective EoS is w_eff=1/3. In addition, the fixed point
(<ref>) has the additional eigenvalue
μ_3=3+3w_0-3ζ_0+w_1κ^2 with w_eff=0, which
describes a matter dominated state. Finally, the fixed point
(<ref>) has μ_3=3 with
w_eff=w_0+1/3w_1κ^2-ζ_0. In Table <ref> we
have gathered all the fixed points for the case q=-3ζ_0.

As it can be seen in Table <ref>, the fixed point P_2a is
always unstable, if it lies in the physical region, that is, when
x_0⩽ 1. The fixed point P_2b may be stable in the past
and in the future, depending on the values of the parameters.
Finally, the fixed point P_2c may be stable only if w_eff>0.
So in this case, no fixed point describes late-time acceleration,
however, the fixed points P_2a and P_2b may describe a
radiation dominated era and matter dominated era respectively.




Consider now the case q=-3ζ_0, q=-w_1κ^2, in which
case, the fixed point (<ref>) has the additional eigenvalue
μ_3=3 and the corresponding effective EoS is w_eff=w_0.
Accordingly, the fixed point (<ref>) has μ_3=3+3w_0 with
w_eff=0 and finally the fixed point (<ref>) has μ_3=4
with w_eff=1/3. Thus in Fig.<ref> we have
w_eff=1/3 for the fixed point (x_0=0,y_0=0),
w_eff=0 for the fixed point (x_0=1,y_0=0) and
w_eff=1/2 for the fixed point (x_0=0,y_0=1).
Correspondingly, in Fig. <ref> we have
w_eff=1/3 for the fixed point (x_0=0,y_0=0),
w_eff=0 for the fixed point (x_0=1,y_0=0) and
w_eff=-1/2 for the fixed point (x_0=0,y_0=1). Note
that all the aforementioned fixed points correspond to states which
have infinite values of H.

In Table <ref> we have gathered all the fixed points for the
case q=-3ζ_0, q=-w_1κ^2.



Now let us consider some alternative choices for the parameter k,
and we start with the case k=1, in which case the only stationary
point is the following,

    x_0=3ζ_0-3-w_1κ^2/3w_0,   y_0=1-x_0,   z_0=qw_0/κ^2(3+3w_0+w_1κ^2-3ζ_0-3w_0ζ_0),

and the expression (<ref>) yields for this fixed point,

    w_eff=-1.

The effective EoS parameter above describes a de Sitter evolution,
with a finite value for the Hubble parameter, namely H=H_0, which
may be arbitrarily small, depending on the values of the other
parameters. The eigenvalues for this fixed point in the most general
case, have a quite complicate form, but there is one special case
for which we can compute these analytically, and this occurs when
w_0=1/3. In this case, the eigenvalues are,

    μ_1      =-4,
    μ_2,3      =1/2[ 3ζ_0-7-w_1κ^2±√((3ζ_0-7-w_1κ^2)^2+12(4-4ζ_0+w_1κ^2)) ],

where the first eigenvalue corresponds to variable x. We can see
that the stability with respect to all dimensions, requires the
following conditions to hold true,

    ζ_0      <7/3+1/3w_1κ^2,
    ζ_0      >1+1/4w_1κ^2,

which is quite compatible with 0⩽ x_0⩽ 1, but
incompatible with z_0>0, which reads
ζ_0<1+1/4w_1κ^2. This means that by varying the
parameters, we can make this point stable with respect to the
coordinates x and y, and unstable with respect z, so this
point may be used for the construction of the late-time acceleration
phase. Note also that by changing value of w_0, it is quite
possible to make this point stable with respect to all the
coordinates.


Consider now the case k=-1, in which case there are four
stationary points, which can be found in Table <ref>. Note that
the parameter a appearing in Table <ref> is equal to
a=(3w_0-w_1κ^2+3ζ_0)^2+12κ^2w_0w_1. The
corresponding eigenvalues can be found in Table <ref>.

By looking the Table <ref> it can be seen that the most
phenomenologically interesting fixed point is P_4b. First of all
it is easy to make this fixed point a stable attractor. Moreover
w_eff is equal to -1 exactly, so this point corresponds to
some de Sitter solution. Finally, for this point we have z_0≠ 0
which means that this point corresponds to some state with non-zero
(and non-infinite) value of the Hubble parameter H=H_0. So by
changing the values of the parameters it is easy to make z_0
sufficiently large, which implies sufficiently small values of
H_0. For instance, by choosing w_0>0, w_1>0 and ζ_0<0 we
can make this fixed point stable. By assigning sufficiently large
values of w_1 for sufficiently small (positive) values of q, we
get large values for z_0. And finally we can see that by changing
the values of the parameters, it is easy to construct some fixed
relation y_0/x_0≡ρ_dm/ρ. In
conclusion, this de Sitter attractor may be viewed as a late-time
attractor.

It is interesting to note that there are five different regions in
the parameter space, for which all four stationary points lie in the
physical region, and these regions are the following,[i.e.
0⩽ x_0⩽ 1, 0⩽ y_0⩽ 1,
z⩾0 for all points simultaneously.]

    w_1κ^2⩽-4,  q>0,   -1/3⩽ζ_0<a,  w_0⩾ b,
    
    w_1κ^2⩽-4,  q>0,   a⩽ζ_0⩽0,  w_0⩾ c,
    
    -4<w_1κ^2<-2/3√(13)-2/3,  q>0,   -c+4/3⩽ζ_0<a,  w_0⩾ b,
    
    -4<w_1κ^2<-2/3√(13)-2/3,  q>0,   a⩽ζ_0<0,  w_0⩾ c,
    
    -2/3√(13)-2/3⩽ w_1κ^2⩽-3,  q>0,   -c+4/3⩽ζ_0⩽ 0,  w_0⩾ c,

where

    a=1/3(2w_1κ^2-1)+2/3√(w_1^2κ^4-w_1κ^2),



    b=-1/3w_1κ^2-ζ_0+2/3√(3w_1κ^2ζ_0),



    c=1/3(1-w_1κ^2).

However, the above cases are quite hard to be tackled analytically,
so a numerical study is needed, which exceeds the purposes of this
article.




 §.§ Oscillating Dark Energy-Dark Matter Interaction



In the previous sections we considered cases of dark energy and dark
matter interactions, by specifying the free parameters, and in this
section we follow a different approach: we will directly modify the
functional form of the dark energy-dark matter interaction by making
it oscillating. So assume that in Eq. (<ref>) we make the
following replacement, -QH^2k→ -Q cos(h_0H^2). In
effect, instead of Eqs. (<ref>)-(<ref>) we have

    x'      =(x-1)(w_1κ^2+3xw_0-x)-x(y+3ζ_0)-qcos(h_0H^2),
    
    y'      =-y^2+y(1-3ζ_0+3xw_0-x+w_1κ^2)+3ζ_0+qcos(h_0H^2),
    
    z'      =z(4-x-y+3xw_0-3ζ_0+w_1κ^2),

where cos(h_0H^2)=cos (h_0/κ^2 z ). Let
us suppose that the argument of the cosine function changes
monotonically, from π/2 during the inflationary era,
when H^2 M_Pl, to 0 which corresponds at the present
time. In practice we may realize this kind of behavior by assuming
sufficiently small values of the parameter h_0. Then, the physical
interpretation of the situation at hand is the following, at early
times (inflationary era), there is no interaction between the
various matter fluid components, and as the time grows, this
interaction starts to develop, and it grows until it reaches its
maximum  value at H=0. From a mathematical point of view, this
would mean that there are two distinct asymptotic states at
t→ 0 and t→∞, which are described by
two distinct dynamical systems. In the following we discuss these
two asymptotic dynamical systems separately.







  §.§.§ Asymptotic State I: The Inflationary Epoch


Let us now focus on the inflationary asymptotic state, in which case
instead of (<ref>)-(<ref>), we have the following
dynamical system,

    x'      =(x-1)(w_1κ^2+3xw_0-x)-x(y+3ζ_0),
    
    y'      =-y^2+y(1-3ζ_0+3xw_0-x+w_1κ^2)+3ζ_0,
    
    z'      =z(4-x-y+3xw_0-3ζ_0+w_1κ^2),

where for the stationary points in the last equation above, we must
put z_0=0. We can see that the last equation is decoupled from the
system, and mainly governs the stability with respect to the
variable z only. The stationary points of the system
(<ref>)-(<ref>) are already found in section <ref>
and these are presented in Eqs. (<ref>)-(<ref>),
where we need to put q=0. The first fixed point should be excluded
from the future investigation, since it has w_eff=1/3.
Thus we have one possible fixed point candidate that may describe
inflation (let us denote it P_5a):

    x_0=1-y_0,  
    y_0=1/6w_0 [ 3w_0+w_1κ^2-3ζ_0 ±√((3w_0+w_1κ^2-3ζ_0)^2+36w_0ζ_0) ],

with,
 
    w_eff=1/3 [ 3x_0w_0 +w_1κ^2 -3ζ_0  ].

Now in order to describe inflation, we need to construct a de Sitter
solution, which should be stable with respect to the coordinates
x, y and unstable with respect to z, in order to provide an
exit from inflation. Thus let us put for (<ref>)
w_eff=-1+γ, where 0<γ≪ 1. In this case, the third
eigenvalue associated with the coordinate z will exactly be
μ_3=γ, so this provides an instability of de Sitter point.
Moreover deriving x_0 from (<ref>) and equating it to
(<ref>), we find an additional relation between the free
parameters,
 
    ζ_0=(1-γ)(3-3γ+3w_0+w_1κ^2)/3(1+w_0-γ),

or if we require that all the values of parameters are not very
small in comparison with γ, we have,
 
    ζ_0=3+3w_0+w_1κ^2/3(1+w_0)-6+3w_0+w_1κ^2/3(1+w_0)γ.





  §.§.§ Asymptotic State II: Present Time and Late-time Acceleration


Now let us consider the present-time asymptotic case, in which case
instead of the dynamical system of Eqs. (<ref>)-(<ref>) we
have the following system,

    x'      =(x-1)(w_1κ^2+3xw_0-x)-x(y+3ζ_0)-q,
    
    y'      =-y^2+y(1-3ζ_0+3xw_0-x+w_1κ^2)+3ζ_0+q,
    
    z'      =z(4-x-y+3xw_0-3ζ_0+w_1κ^2),

where for stationary points in the last equation we must require
z_0≠ 0. Thus, from the last equation we have
 
    y_0=4-3ζ_0+w_1κ^2+x_0(3w_0-1).

Substituting the above in Eq. (<ref>) we find,

    x_0=-q-w_1κ^2/3(w_0+1),

and by substituting in Eq. (<ref>) we find,

    x_0=12ζ_0+q-3w_1κ^2-12/3(3w_0-1).

In effect, the system will compatible only if the following
condition holds true,

    -q-w_1κ^2/3(w_0+1)=12ζ_0+q-3w_1κ^2-12/3(3w_0-1),

which gives us,
 
    q=3w_0+3-3ζ_0+w_1κ^2-3ζ_0w_0/w_0,

and the corresponding stationary point is (let us denote it
P_5b)
 
    x_0=3ζ_0-w_1κ^2-3/3w_0,  y_0=3w_0+3-3ζ_0+w_1κ^2/3w_0,  x_0+y_0=1.

Combining (<ref>) with (<ref>) we find,
 
    q=(6+3w_0+w_1κ^2)γ,

so we can see that q is very small but does not
vanish.[Putting w_eff=-1 exactly for the first point,
we get q=0.] Note also that the effective EoS for this point,
calculated for (<ref>) yields exactly w_eff=-1 for any
values of the parameters.






  §.§.§ combining and stability analysis.



By looking the differential equations in Eqs.
(<ref>)-(<ref>) and (<ref>)-(<ref>), it can be
seen that these are identical apart for some additive constants, and
in effect the eigenvalues will be identical for the two systems. Let
us denote x'=f and y'=g for notational simplicity. In effect,
the equation that determines the eigenvalues, takes the following
form,

    |
    [   (f_x)_0-μ (f_y)_0;                    ; (g_x)_0   (g_y)_0-μ ] | =0,

where (f_x)_0, (f_y)_0, (g_x)_0 and (g_y)_0, are equal to,

    (f_x)_0=w_1κ^2+6w_0x_0-x_0-3ζ_0-3w_0,



    (f_y)_0=-x_0,



    (g_x)_0=(3w_0-1)(1-x_0),



    (g_y)_0=-1+x_0-3ζ_0+3w_0x_0+w_1κ^2,

and in the above relations we took into account that for both points
we have x_0+y_0=1.

The solution of Eq. (<ref>) for the fixed point P_5a with
positive sign in (<ref>) reads,
 
    μ_1=-√((3w_0+w_1κ^2-3ζ_0)^2+36w_0ζ_0),  μ_2=1/2[ -2+ 3w_0+w_1κ^2-3ζ_0 +μ_1 ],

and for the point P_5b
 
    μ_1=-4,  μ_2= -6 - 3w_0 - w_1κ^2+3ζ_0,

and by substituting the expression (<ref>) for ζ_0, for
the fixed point P_5a, we have,
 
    μ_1=-3 + 6w_0 + 3w_0^2 + w_0w_1κ^2/w_0+1,  μ_2=-4,

and for the point P_5b we get,
 
    μ_1=-4,  μ_2= -3 + 6w_0 + 3w_0^2 +
    w_0w_1κ^2/w_0+1.

Note that we have set γ=0, since all the eigenvalues, even in
this case, do not vanish. Thus we can see that both points P_5a
and P_5b are stable or unstable with respect to coordinates x
and y simultaneously, and stability condition reads,
 
    (3 + 6w_0 + 3w_0^2 + w_0w_1κ^2)(w_0+1)>0.

Moreover, we need to require that both fixed points lie in the
physical region, namely at 0⩽ x_0⩽ 1 and
0⩽ y_0⩽ 1, so the following conditions must be
satisfied,

    w_1κ^2<-3,  w_0⩾ c,      
    
    w_1κ^2=-3,  w_0> 0,      
    
    -3<w_1κ^2<0,  c⩽ w_0 <0  w_0>0,      
    
    w_1κ^2>0,  w_0⩽ c,

where,

    c=-1/3(w_1κ^2+3).

In Table <ref> we have gathered all the fixed points and the
corresponding eigenvalues for the case of an oscillating form of the
interaction term Q.

It is worth discussing some interesting scenarios, so assume that
w_0=1, w_1κ^2=-1, then the fixed points and the eigenvalues
become,

    P_5a:   x_0=1/6,  y_0=5/6,   μ_1=-11/2,  μ_2=-4,



    P_5b:   x_0=5/6,  y_0=1/6,   μ_1=-4,  μ_2=-11/2.

So by changing the values of the parameters w_0 and w_1, it is
possible to appropriately fix the fixed point P_5b, which recall
that it corresponds to late-time acceleration, in order some fixed
relation between ρ_dm and ρ_de holds true. Note that
in this section we studied only case for which the choices of the
functions f(ρ) and G(H^2) were the simplest choices, but in
principle a more involved functional form for these functions may
lead to more interesting phenomenology. In such a case though, the
analytical treatment will be possibly insufficient, so a concrete
numerical analysis will be needed.


In conclusion, we constructed a cosmological model that describes
the inflationary and the late-time acceleration era, due to an
oscillating interaction term between the dark energy and dark matter
fluids. We found two fixed points, with the first being P_5a,
which describes the initial de Sitter solution, namely the
inflationary era, and this fixed point was unstable, a feature which
indicates that the graceful exit is triggered. The second fixed
point was also found to be a de Sitter fixed point, which describes
the present day acceleration. As we demonstrated, in the oscillating
model we discussed, during the early Universe there was no
interaction between dark energy and dark matter, however during the
late-time the interaction was present. Finally, as we showed, by
appropriately choosing the parameters, it is possible to produce
some fixed relation for ρ_dm/ρ_de.



Before closing this section, let us briefly discuss an interesting
issue with regards the early-time era. It would be interesting to
calculate the slow-roll indices and the corresponding observational
indices for the inflationary era we presented in section <ref>.
For example by assuming the perfect fluid approach
<cit.>, we can express the spectral index of primordial
curvature perturbations n_s and the scalar-to-tensor ratio r in
the usual way these are given in the case of a canonical scalar
field,

    n_s=1-6ϵ+2η,



    r=16ϵ.

with the slow-roll indices being defined in terms of the Hubble rate
as follows,

    ϵ= -Ḣ/H^2,



    η= ϵ -Ḧ/2HḢ.

By combining Eqs. (<ref>) and (<ref>) it is easy to
calculate the slow-roll parameter ϵ, which reads,

    ϵ=1/2 [ 4- y +x(3w_0-1) +w_1κ^2 -3ζ_0  ].

Differentiating (<ref>) and by combining with equations
(<ref>)-(<ref>) we can find a similar expression for the
slow-roll parameter η. Note here that the resulting expression
for η is much more complicated so we omit it. There is a major
obstacle in calculating the slow-roll indices however, since in the
expressions for ϵ and η, we must use values for x
and y not on the stationary point, but from some point near the
stationary point. The resulting slow-roll indices must be expressed
in terms of the e-foldings number N, but doing this analytically
is a rather formidable task, that exceeds the purposes of this work.
A numerical approach though might be less difficult to perform, so
we hope to address this task in a future work.




§ CONCLUSIONS
 



In this work we analyzed in detail the phase space of a cosmological
system that contains cosmological fluids that have various forms of
equation of state. We firstly discussed the simplest forms of EoS,
and we found the fixed points of the cosmological dynamical system
and we discussed the physical significance of these fixed points. In
addition, we discussed theories that admit interactions between the
dark energy and dark matter fluids. In addition we introduced a new
class of interaction between dark energy and dark matter, in which
theories the interaction term is oscillating, allowing different
form of interactions for various eras during the cosmological
evolution. As we demonstrated, it is possible to have almost
negligible interactions at early times, that is, during the
inflationary era, and for the same model the interaction is turned
on at late times. The cosmological dynamical system of the
oscillating interaction term is in turn decomposed into two distinct
dynamical systems at early and late times. Interestingly enough the
two dynamical systems predict two de Sitter fixed points
corresponding to early and late times, with the early-time de Sitter
point being unstable in one of the coordinates, a feature which
indicates the possible exit from the inflationary era. This
framework gives us the interesting fluid-filled Universe evolution,
unifying the early-time acceleration with late-time acceleration.

The fluid description offers an alternative viewpoint in modern
cosmology, which may describe successfully many eras of our
Universe's evolution. The dynamical system approach offers many new
insights since the fixed points of the dynamical system reveal the
attractors of the whole theory and their stability indicates whether
these attractors are final attractors of the system. A compelling
extension of this work is to include Loop Quantum Cosmology effects
in the EoS, as was performed in Ref. <cit.>, so
we defer this task to a future work. Another strong motivation to
adopt the fluid approach in cosmology, is the late-time acceleration
era, since this can be modeled by fluid cosmology
<cit.>. Furthermore, as modified gravity maybe easily
presented in the effective fluid representation, the dynamical
systems approach turns out to be useful also for the study of the
cosmology in modified gravity. More interestingly, the dynamical
system approach can be applied to the study of bouncing cosmology,
described by a multi-component fluid. We aim to address this latter
issue in a future work. In conclusion, we believe that the
multi-component fluid cosmology may have some valuable insights to
offer to the cosmologists community.





This work was supported by the Russian Science Foundation (RSF)
grant 16-12-10401 (P.V.T.), by the Russian Foundation for Basic Research (RFBR) grant 17-02-01008  (P.V.T.),
by MINECO (Spain), project
FIS2013-44881, FIS2016-76363-P and by CSIC I-LINK1019 Project
(S.D.O) and by the Russian Ministry of Education and Science,
Project No. 3.1386.2017 (S.D.O and V.K.O).



99

riess A.G. Riess et al. (High-z Supernova Search Team), Astronom. J. 116, 1009 (1998) [arXiv:astro-ph/9805201]


reviews1

S. Nojiri, S. D. Odintsov and V. K. Oikonomou,
  
  arXiv:1705.11098 [gr-qc].
  




reviews2


S. Nojiri and S. D. Odintsov,
  
  Phys. Rept.  505 (2011) 59
  doi:10.1016/j.physrep.2011.04.001
  [arXiv:1011.0544 [gr-qc]]



reviews3

S. Nojiri and S. D. Odintsov,
  
  eConf C 0602061 (2006) 06
   [Int. J. Geom. Meth. Mod. Phys.  4 (2007) 115]
  doi:10.1142/S0219887807001928
  [hep-th/0601213].



reviews4
  V. Faraoni and S. Capozziello,
  
  Fundam. Theor. Phys.  170 (2010).
  doi:10.1007/978-94-007-0165-6




inhomogen1

S. Nojiri and S. D. Odintsov, Phys. Rev. D 72 (2005) 023003
[hep-th/0505215]


nojodineos1

K. Bamba, S. Nojiri, S. D. Odintsov and D. Sáez-Gómez,
Phys. Rev. D 90 (2014) 124061 [arXiv:1410.3993 [hep-th]].







inhomogen2

S. Capozziello, V. F. Cardone, E. Elizalde, S. Nojiri and
S. D. Odintsov, Phys. Rev. D 73 (2006) 043512
[astro-ph/0508350].

brevik1


I. Brevik and Ø. Grøn,

arXiv:1409.8561 [gr-qc].





brevik2


M. Szydlowski, W. Godlowski and R. Wojtak,

Gen. Rel. Grav. 38 (2006) 795 [astro-ph/0505202].








fluid1
X. Chen, H. Firouzjahi, M. H. Namjoo and M. Sasaki,

JCAP 1309 (2013) 012 [arXiv:1306.2901 [hep-th]]

fluid2 R. Myrzakulov and L. Sebastiani,

Astrophys. Space Sci. 356 (2015) 1, 205 [arXiv:1410.3573
[gr-qc]]

Brevik:2017msy
  I. Brevik, O. Gron, J. de Haro, S. D. Odintsov and E. N. Saridakis,
  
  arXiv:1706.02543 [gr-qc].
  





oikonomouimperfect

 V. K. Oikonomou,
  
  arXiv:1703.09009 [gr-qc].
  


oikonomoufluid

 S. Nojiri, S. D. Odintsov and V. K. Oikonomou,
  
  Phys. Lett. B 747 (2015) 310
  
  [arXiv:1506.03307 [gr-qc]].
  
  



Brevik:2016kuy
  I. Brevik and A. V. Timoshkin,
  
  Int. J. Geom. Meth. Mod. Phys.  14 (2017) no.04,  1750061
  
  [arXiv:1612.06689 [gr-qc]].
  
  


ALUL
A. Avelino, Y. Leyva and L. A. Urena-Lopez, Phys. Rev. D 88, 123004
(2013).

CLLPS
N. Cruz, S. Lepe, Y. Leyva, F. Pena and J. Saavedra, Phys. Rev. D
90, 083524 (2014).

LS
 Y. Leyva and M. Sepulveda, arXiv:1611.03833 [gr-qc].

VWM
 H. Velten, J. Wang and X. Meng, Phys. Rev. D 88,
123504 (2013).

VS
 H. Velten and D. Schwarz, Phys. Rev. D 86, 083501
(2012).

VCFCB
 H. Velten, T. R. P. Carames, J. C. Fabris, L. Casarini
and R. C. Batista, Phys. Rev. D 90, 123526 (2014).

VSFZ
 H. Velten, D. J. Schwarz, J. C. Fabris and W. Zimdahl,
Phys. Rev. D 88, 103522 (2013).

V
 H. Velten, Int. J. Geom. Meth. Mod. Phys. 11, 1460013
(2014).

IS
 W. Israel and J. M. Stewart, Annals Phys. 118, 341
(1979).

M
 R. Maartens, astro-ph/9609119.

CLP
 N. Cruz, S. Lepe and F. Pena, Phys. Lett. B 767, 103
(2017).

CCL
 M. Cruz, N. Cruz and S. Lepe, Phys. Lett. B 769, 159
(2017).

LOS
S. Lepe, G. Otalora, J. Saavedra ,
arXiv:1704.05625 [gr-qc].



ZP
 W. Zimdahl and D. Pavon, Phys. Lett. B 521, 133
(2001).



CCMUL
 G. Caldera-Cabral, R. Maartens and L. A. Urena-Lopez,
Phys. Rev. D 79, 063518 (2009); AIP Conf. Proc. 1241, 741 (2010).








Coley:1999uh
  A. A. Coley,
  
  gr-qc/9910074.; 

  
  
G. Leon and C. R. Fadragas,
  
  arXiv:1412.5701 [gr-qc].
  
  




GarciaSalcedo:2012dn
  R. Garcia-Salcedo, T. Gonzalez and I. Quiros,
  
  arXiv:1211.2738 [gr-qc].
  
  



  
Boehmer:2008av
  C. G. Boehmer, G. Caldera-Cabral, R. Lazkoz and R. Maartens,
  
  Phys. Rev. D 78 (2008) 023505
  doi:10.1103/PhysRevD.78.023505
  [arXiv:0801.1565 [gr-qc]].
  
  


  
Boehmer:2014vea
  C. G. Boehmer and N. Chan,
  
  
  arXiv:1409.5585 [gr-qc].
  
  


  
Haba:2016swv
  Z. Haba, A. Stachowski and M. Szydlowski,
  
  JCAP 1607 (2016) no.07,  024
  doi:10.1088/1475-7516/2016/07/024
  [arXiv:1603.07620 [gr-qc]].
  
  


  
Bolotin:2013jpa
  Y. L. Bolotin, A. Kostenko, O. A. Lemets and D. A. Yerokhin,
  
  Int. J. Mod. Phys. D 24 (2014) no.03,  1530007
  doi:10.1142/S0218271815300074
  [arXiv:1310.0085 [astro-ph.CO]].
  
  








Abdelwahab:2007jp
  M. Abdelwahab, S. Carloni and P. K. S. Dunsby,
  
  Class. Quant. Grav.  25 (2008) 135002
  doi:10.1088/0264-9381/25/13/135002
  [arXiv:0706.1375 [gr-qc]].
  
  



Carloni:2007eu
  S. Carloni, S. Capozziello, J. A. Leach and P. K. S. Dunsby,
  
  Class. Quant. Grav.  25 (2008) 035008
  doi:10.1088/0264-9381/25/3/035008
  [gr-qc/0701009].
  
  




Carloni:2004kp
  S. Carloni, P. K. S. Dunsby, S. Capozziello and A. Troisi,
  
  Class. Quant. Grav.  22 (2005) 4839
  doi:10.1088/0264-9381/22/22/011
  [gr-qc/0410046].
  
  




Leon:2014yua
  G. Leon and E. N. Saridakis,
  
  JCAP 1504 (2015) no.04,  031
  doi:10.1088/1475-7516/2015/04/031
  [arXiv:1501.00488 [gr-qc]].
  
  





Xu:2012jf
  C. Xu, E. N. Saridakis and G. Leon,
  
  JCAP 1207 (2012) 005
  doi:10.1088/1475-7516/2012/07/005
  [arXiv:1202.3781 [gr-qc]].
  
  



  
Carloni:2017ucm
  S. Carloni and J. P. Mimoso,
  
  arXiv:1701.00231 [gr-qc].
  
  




  
Leon:2010pu
  G. Leon and E. N. Saridakis,
  
  Class. Quant. Grav.  28 (2011) 065008
  doi:10.1088/0264-9381/28/6/065008
  [arXiv:1007.3956 [gr-qc]].
  
  




Carloni:2007br
  S. Carloni, A. Troisi and P. K. S. Dunsby,
  
  Gen. Rel. Grav.  41 (2009) 1757
  doi:10.1007/s10714-008-0747-9
  [arXiv:0706.0452 [gr-qc]].
  
  


  
Khurshudyan:2016qox
  M. Khurshudyan,
  
  Int. J. Geom. Meth. Mod. Phys.  14 (2016) no.03,  1750041.
  doi:10.1142/S0219887817500414
  


Hohmann:2017jao
  M. Hohmann, L. Jarv and U. Ualikhanova,
  
  arXiv:1706.02376 [gr-qc].
  



Odintsov:2015wwp
  S. D. Odintsov and V. K. Oikonomou,
  
  Phys. Rev. D 93 (2016) no.2,  023517
  doi:10.1103/PhysRevD.93.023517
  [arXiv:1511.04559 [gr-qc]].
  
  



STT

  M. Skugoreva, A. Toporensky and P. Tretyakov,
  
  Grav. Cosmol.  17 (2011) 110
  
  [arXiv:1007.3365 [gr-qc]].
  
  

NO1
S. Nojiri and S. D. Odintsov,
  
  Phys. Rev. D 72 (2005) 023003
  doi:10.1103/PhysRevD.72.023003
  [hep-th/0505215].
  
  



Bamba:2012cp
  K. Bamba, S. Capozziello, S. Nojiri and S. D. Odintsov,
  
  Astrophys. Space Sci.  342 (2012) 155
  doi:10.1007/s10509-012-1181-8
  [arXiv:1205.3421 [gr-qc]].
  
  



Misner:1967uu
  C. W. Misner,
  
  Astrophys. J.  151 (1968) 431.
  doi:10.1086/149448
  
  


Bamba:2016aoo
  K. Bamba,
  
  Int. J. Geom. Meth. Mod. Phys.  13 (2016) no.06,  1630007
 
  [arXiv:1604.02632 [gr-qc]].
  
  




Brevik:2005bj
  I. H. Brevik and O. Gorbunova,
  
  Gen. Rel. Grav.  37 (2005) 2039
  doi:10.1007/s10714-005-0178-9
  [gr-qc/0504001].
  
  



  
Brevik:2007ig
  I. H. Brevik, E. Elizalde, O. Gorbunova and A. V. Timoshkin,
  
  Eur. Phys. J. C 52 (2007) 223
  doi:10.1140/epjc/s10052-007-0357-9
  [arXiv:0706.2072 [gr-qc]].
  
  




Normann:2016jns
  B. D. Normann and I. Brevik,
  
  Entropy, Vol. 18, article 215 (2016)
  doi:10.3390/e18060215
  [arXiv:1601.04519 [gr-qc]].
  
  



BH
 J. D. Barrow and S. Hervik,
  
  Phys. Rev. D 74 (2006) 124017
  doi:10.1103/PhysRevD.74.124017
  [gr-qc/0610013].
  
  





Oikonomou:2014lua
  V. K. Oikonomou and N. Karagiannakis,
  
  Class. Quant. Grav.  32 (2015) no.8,  085001
  doi:10.1088/0264-9381/32/8/085001
  [arXiv:1408.5353 [gr-qc]].
  
  








Balakin:2010xb
  A. B. Balakin and V. V. Bochkarev,
  
  Phys. Rev. D 83 (2011) 024036
  doi:10.1103/PhysRevD.83.024036
  [arXiv:1012.2433 [gr-qc]].;
N. Cruz, S. Lepe and F. Pena,
  
  Phys. Lett. B 663 (2008) 338
  doi:10.1016/j.physletb.2008.04.035
  [arXiv:0804.3777 [hep-ph]].
  
  
  
  








