Ann. Acad. Rom. Sci.

Ser. Math. Appl.

ISSN 2066-6594

Vol. 18, No. 2/2026

DYNAMICS OF THE ΛCDM MODEL OF THE

UNIVERSE FROM THE PERSPECTIVE OF

THE DYNAMICAL SYSTEMS THEORY^∗

‡

Danijela Brankovi´c^†

Zarko Mijajlovi´c

Marija Miki´c^§

ˇ

Communicated by B. Sasu

DOI

10.56082/annalsarscimath.2026.2.5

Abstract

We use the theory of the dynamical systems to study the dynamics

of the ΛCDM model of the universe. A three-dimensional autonomous

dynamical system of a class of Lotka-Volterra dynamical systems, with

density parameters of the universe’s constituents as dependent vari-

ables is inferred. It is assumed that the matter content in our universe

consists of barotropic perfect ﬂuids without mutual interaction (ex-

cept gravitational). The appropriate physical interpretation, as well

as parametrization between the scale expansion factor of the universe

and cosmological density parameters is presented.

Keywords: Lotka-Volterra dynamical systems, Friedmann equations,

ΛCDM model, dynamical systems, stability, equilibriums, Lyapunov func-

tion.

MSC: 34A34, 83F05, 37C75, 34C23.

^∗Accepted for publication on October 21, 2025

^†danijela@etf.bg.ac.rs, School of Electrical Engineering, University of Belgrade,

Bulevar kralja Aleksandra 73, 11120 Belgrade, Serbia

^‡zarkom@matf.bg.ac.rs, Faculty of Mathematics, University of Belgrade, Studentski

trg 16, 11158 Belgrade, Serbia

^§marija.mikic@matf.bg.ac.rs, Faculty of Mathematics, University of Belgrade, Stu-

dentski trg 16, 11158 Belgrade, Serbia

5

Dynamics of the ΛCDM model of the universe

6

1 Introduction

Our starting point is the system of the Friedmann equations [1] with the

cosmological constant Λ

ꢀ ꢁ

2

a˙

8πG

kc²

a²

Λc²

3

Λc²

=

ρ −

+

,

a

3

ꢀ

ꢁ

a¨

4πG

3p

c²

(1)

= −

ρ +

,

a

3

ꢂ

ꢃ

a˙

p

ρ˙ + 3

ρ +

= 0,

a

c²

consisting of nonlinear diﬀerential equations of the ﬁrst and second order

that describes the dynamics of our universe. The ﬁrst equation in (1) is

called the Friedmann equation, while the second and the third one are ref-

ered to the acceleration equation and the ﬂuid equation (also known as

conservation equation), respectively. Cosmological parameters a = a (t),

ρ = ρ (t) and p = p (t) are functions of time variable t and they represent

the scale expansion factor of the universe, the average density of matter

content in the universe and the pressure of matter content in the universe,

respectively. The curvature of the universe is a real number denoted with k.

In this paper we discuss only ﬂat universe (k = 0) and open universe (k < 0).

Derivatives with respect to t are denoted with a dot. Considering physical

meaning of the parameters, we assume that all functions that appear in this

paper are continuously diﬀerentiable as many times as needed.

As it is well known, the system (1) consists from two independent equa-

tions that contain three unknown functions a = a (t), ρ = ρ (t) and p = p (t).

In order to have a determined system, standard approach is assuming an

equation of state for the matter content that builds up the universe. Here we

consider that the matter components in our universe are barotropic perfect

ﬂuids (”matter” and ”radiation”) which do not interact with each other,

more precisely that a single ﬂuid component (for example ”matter”) do not

interact with the other ﬂuid component (”radiation”). By ”matter” it is typ-

ically considered baryonic matter plus cold dark matter, while ”radiation”

refers to relativistic particles, such as neutrinos and photons.

The equation of state for barotropic perfect i-ﬂuid is a linear relation

between p_iand ρ_i

p_i= ω_iρ_ic²,

(2)

where ρ_i= ρ_i(t) is the average density, p_i= p_i(t) is the pressure and ω_i

is the equation of state parameter of i-ﬂuid [2], [3]. In this paper we take

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

7

that ω_iis a constant from the interval [0, 1], that is consistent with the

macroscopic physics [2].

The assumption that a single ﬂuid component do not interact with the

other ﬂuid component is reﬂecting in that the ﬂuid equation is conserved

for every i-ﬂuid [2], [4], i.e.

ꢂ

ꢃ

a˙

p_i

c²

ρ˙_i+ 3

ρ_i+

= 0.

(3)

a

Substituting (2) and the Hubble parameter H(t) = a˙(t)/a(t) in (3), we infer

the following form of the ﬂuid equation for i-ﬂuid

ρ˙_i= −3Hρ_i(1 + ω_i) .

(4)

The equation (4) in the terms of ρ_m(the average density of m-ﬂuid (matter))

and ρ_r(the average density of r-ﬂuid (radiation)), is in the form

ρ˙ = −3Hρ_m(1 + ω_m) ,

m

(5)

ρ˙_r= −3Hρ_r(1 + ω_r) .

Clearly, ρ_m= ρ_m(t) and ρ_r= ρ_r(t) contribute to the average density of

the matter content in the universe, as well as p_m= p_m(t) and p_r= p_r(t)

contribute to the pressure of the matter content in the universe, with

ρ = ρ_m+ ρ_rand p = p_m+ p_r.

(6)

Cosmological parameters that play a central role in this paper are density

parameters of the material in the universe, deﬁned as (see [5])

8πG

3H²(t)

8πG

kc²

a²(t)H²(t)

Ω_m(t) =

ρ_m(t), Ω_k(t) = −

,

(7)

Λc²

Ω_r(t) =

ρ_r(t),

Ω_Λ(t) =

,

3H²(t)

where Ω_m(t), Ω_k(t), Ω_r(t) and Ω_Λ(t) denote matter density, spatial curvature

density, radiation density and cosmological constant density. Throughout

the paper, Ω_i= Ω_i(t) will stand for any of the following density parameters

Ω_m(t), Ω_k(t), Ω_r(t), Ω_Λ(t), unless otherwise noted.

The paper is organized in the following way. In Section 2 we prove that

the system of the Friedmann equations (1) is equivalent to the system of the

ﬁrst order nonlinear diﬀerential equations whose dependent variables are

exactly density parameters of the material in the universe. In Section 3 we

Dynamics of the ΛCDM model of the universe

8

obtained analytical solution of the obtained dynamical system from where

we obtain a new parametrization of density parameters with respect to a (t).

Moreover, in Section 3 we considered in details all two-dimensional and

one-dimensional cases inferred from the initial three-dimensional dynamical

system by allowing that some of the density parameters Ω_iare equal to zero,

which, to our knowledge, is not presented anywhere else in the literature. In

Section 4 we analyzed the stability of equilibriums of the obtained dynamical

system and we discussed the obtained results in the spirit of the universe’s

evolution. In Section 5 the stability analysis of equilibriums of the observed

dynamical system is conducted when ω_m= ω_r, which is, to our knowledge,

also a new result.

We note that in the paper [6] is analyzed the dynamics of the n-dimensio-

nal generalization of the ΛCDM model that is analyzed in this paper. Con-

sequently, the solutions of the ΛCDM model and stability of its equilibriums

in this paper indeed can be obtained from the generalized n-dimensional case

from the paper [6].

However, the results obtained in the paper [6] are more signiﬁcant from

the mathematical aspect, than from the aspect of physics, whereas the sit-

uation is a bit diﬀerent when it comes to this paper. This paper deals with

the exact ΛCDM model and therefore we can, at any time, compare our

conclusions with the current observations, whereas that was not possible

with the n-dimensional universe in the paper [6].

Moreover, every case that is analyzed in the Section 3 is supported with

the appropriate physical interpretation, which is not the case in the paper

[6]. Furthermore, in the Section 4 we analyzed some interesting dynamical

systems, whose phase spaces are illustrated and the conclusions obtained an-

alyzing their phase trajectories are consistent with the current observations.

Those analysis were not implemented in the paper [6].

In addition, in the paper [6] are considered only the cases ω_m= ω_r, for

i = j, whereas in the Section 5 is analyzed the case ω_m= ω_r.

2 Friedmann equations as a dynamical system

By dividing the Friedmann equation in (1) with H²= 0 (H = 0 implies

static universe) and substituting ρ = ρ_m+ ρ_r, we obtain

Ω_m+ Ω_r+ Ω_k+ Ω_Λ= 1.

(8)

According to (7) and (8), as well as the assumption that k ≤ 0, it is clear

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

9

that

0 ≤ Ω_i≤ 1.

(9)

From (8) directly follows

Ω_k= 1 − Ω_m− Ω_r− Ω_Λ.

(10)

In order to infer the dynamical system, we present a method which is pic-

tured in various literature (for example see [2], [4], [7], [8], [9]) in a somewhat

diﬀerent way.

We perform change₀of variables ξ(t) = ln(a(t)), denoting the derivatives

with respect to ξ with . Note that ξ better approximates linearity of time

t than a(t). Diﬀerentiating the density parameters Ω_m, Ω_rand Ω_Λwith

respect to ξ, we infer

˙

dΩ_idt

Ω_i

0

Ω_i=

·

=

.

(11)

dt

dξ

H

Substituting Ω_i= 8πGρ_i/3H²and (4) in (11), for i ∈ {m, r} we obtain

˙

−2H

0

Ω_i= −3Ω_i(1 + ω_i) + Ω_i·

.

(12)

H²

If we replace ρ = ρ_m+ ρ_rand H = a˙/a in the Friedmann equation in (1)

and then diﬀerentiate it with respect to the time t, we infer

8πG

a˙

a³

(ρ˙ + ρ˙_r) + 2kc²·

.

(13)

˙

2HH =

m

3

Furthermore, substitution of (5), (7) and (10) in (13) transforms previous

equation into

˙

−2H

= 2 + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ.

(14)

H²

Replacing (14) into (12), for i ∈ {m, r} we get

0

Ω_i= −3Ω_i(1 + ω_i) + Ω_i(2 + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) . (15)

For Ω_Λ= Λc²/3H², we have

˙

−2H

0

Ω_Λ= Ω_Λ·

,

(16)

(17)

H²

wherefrom, using (14), follows

0

Ω_Λ= Ω_Λ(2 + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) .

Dynamics of the ΛCDM model of the universe

10

Gathering (15) and (17) together, we ﬁnally infer the dynamical system

0

Ω_m= Ω_m(− (1 + 3ω_m) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

0

Ω_r= Ω_r(− (1 + 3ω_r) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

(18)

0

Ω_Λ= Ω_Λ(2 + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

which is the system of the ﬁrst order nonlinear diﬀerential equations whose

dependent variables are density parameters of the material in the universe.

In Theorem 1 we prove that the system of the Friedmann equations (1) is

equivalent to the system of the ﬁrst order nonlinear diﬀerential equations

(18) under certain conditions.

Theorem 1. Assuming p_m= ω_mρ_mc², p_r= ω_rρ_rc²and the identities (6),

the system of the Friedmann equations (1) together with the equations (5)

are equivalent to the system

0

Ω_m= Ω_m(− (1 + 3ω_m) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

0

Ω_r= Ω_r(− (1 + 3ω_r) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

0

Ω_Λ= Ω_Λ(2 + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

(19)

Ω_m+ Ω_r+ Ω_Λ+ Ω_k= 1,

˙

Ω_i

0

Ω_i=

,

H

where Ω_iare deﬁned in (7).

We already inferred that from the system (1) and the equations (5) un-

der the assumptions from Theorem 1 follows the system (19). The opposite

direction is left to be proved. However, we note that Theorem 1 is a spe-

cial case of the Theorem 2.1 in [6] that guarantees the equivalence of the

appropriate n-dimensional dynamical system and the system of Friedmann

equations (1) for n barotropic perfect non-interacting ﬂuids. Thus the proof

of Theorem 1 is omitted here, since it can be obtained from Theorem 2.1 in

[6] for n = 3.

The obtained equivalence guarantees that dynamics of the systems (19)

and (1) are the same, under standard assumptions in cosmology. We note

that in literature the same conclusion is used but only with ω_m= 0 and

ω_r= 1/3, however only one implication is inferred, not the equivalence.

3 Solutions of the dynamical system

We see that the system (19) is of a class of Lotka-Volterra dynamical systems,

which were introduced by Lotka and Volterra (see [10, 11]). Their analytical

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

11

solution is not easy to obtain and so they are mostly solved numerically.

Here we infer analytical solution of the system (19) in various cases,

depending on whether some of the density parameters Ω_iare equal to zero.

Every case that we analyze here corresponds to some model of the universe.

More precisely, it represents some stage, or perhaps some medium stage in

the evolution of the universe. To each case we give the appropriate physical

interpretation [12].

3.1

Three dimensional dynamical system

Here we consider that Ω_m> 0, Ω_r> 0, Ω_Λ> 0 and Ω_k> 0. Therefore, this

case represents an open universe in the matter dominated epoch. Moreover,

it represents an open universe in the phase of the transition between radi-

ation and dark energy that occurred approximately around 5 · 10⁸yrs after

the Big Bang [13].

First, we rewrite the system (19) in the following form

0

2

Ω_m= Ω_m(−(1 + 3ω_m) + (1 + 3ω_r)Ω_r− 2Ω_Λ) + (1 + 3ω_m)Ω_m

,

0

2

Ω_r= Ω_r(−(1 + 3ω_r) + (1 + 3ω_m)Ω_m− 2Ω_Λ) + (1 + 3ω_r)Ω_r,

(20)

0

2

Ω_Λ= Ω_Λ(2 + (1 + 3ω_m)Ω_m+ (1 + 3ω_r)Ω_r) − 2Ω_Λ

,

˙

Ω_i

0

with Ω_m+ Ω_r+ Ω_Λ+ Ω_k= 1 and Ω_i=

. Since Ω_m, Ω_r, Ω_Λand Ω_kare

H

strictly positive and their sum is equal to 1, hence all density parameters

are less than 1.

Note that the density parameters can be observed as a functions of ξ =

ln a. To simplify the system (20), we introduce new variables with

z_m(ξ) = 1/Ω_m(ξ),

z_r(ξ) = 1/Ω_r(ξ),

z_Λ(ξ) = 1/Ω_Λ(ξ).

(21)

Replacing the density parameters in (20) with new variables z_m, z_rand

z_Λ, we infer less complicated system of equations

ꢀ

ꢁ

1 + 3ω_r

2

0

z_m= −z_m−(1 + 3ω_m) +

−

− (1 + 3ω_m),

z_r

z_Λ

ꢀ

ꢁ

1 + 3ω_m

2

0

z_r= −z_r−(1 + 3ω_r) +

−

− (1 + 3ω_r),

(22)

z_m

z_Λ

ꢀ

ꢁ

1 + 3ω_m

1 + 3ω_r

0

z_Λ= −z_Λ2 +

+

+ 2.

z_m

z_r

The next step is solving the system (22) by the elimination method.

Dynamics of the ΛCDM model of the universe

12

The second and the third derivative of z_m₀with respect to the variable ξ,

after substituting the expressions for z_m⁰, z_rand z_Λfrom (22), are in the

0

form

z_m(1 + 6ω_m− 3ω_r) (1 + 3ω_r)

00

z_m= −(1 + 3ω_m)²−

z_r

ꢂ

ꢃ

z_m8 + 12ω_m+ (1 + 3ω_m)²z_Λ

+

,

z_Λ

z_m= − (1 + 3ω_m)³

000

(23)

ꢄ

ꢅ

z_m(1 + 3ω_r) 1 + 27ω_m²+ ω_m(9 − 27ω_r) − 3ω_r+ 9ω_r

2

−

z_r

ꢂ

ꢃ

z_m26 + 72ω_m+ 54ω_m²+ (1 + 3ω_m)³z_Λ

+

.

z_Λ

0

Now, we want to express z_rand z_Λin the terms of z_m, z_mand z_m⁰⁰. In

order to do so, we observe the function F = F (z_m, z_m⁰, z_m⁰⁰, z_r, z_Λ), given

with

ꢄ

ꢅ

ꢄ

ꢅ

ꢄ

ꢅꢅ

F z_m, z_m⁰, z_m⁰⁰, z_r, z_Λ= F₁z_m, z_m⁰, z_r, z_Λ, F₂z_m, z_m⁰⁰, z_r, z_Λ

, (24)

where

ꢀ

ꢁ

1 + 3ω_r

2

F₁= z_m⁰+ z_m−(1 + 3ω_m) +

−

+ (1 + 3ω_m) = 0,

z_r

z_Λ

z_m(1 + 6ω_m− 3ω_r) (1 + 3ω_r)

F₂= z_m⁰⁰+ (1 + 3ω_m)²+

(25)

z_r

ꢂ

ꢃ

z_m8 + 12ω_m+ (1 + 3ω_m)²z_Λ

−

= 0.

z_Λ

According to z_m= 1/Ω_m> 1, ω_r∈ [0, 1] and our assumption that all

functions that appear in this paper are continuously diﬀerentiable as many

times as needed, as well as

ꢆ

∂F

∂F₁

1

2

6 (1 + 3ω_r) (1 + ω_r) z_m

∂z_r

∂F₂

∂z_Λ

∂F₂

= −

= 0,

(26)

2

z_rz_Λ

∂z_r

∂z_Λ

the implicit function theorem allows us to z_rand z_Λexpress as functions of

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

13

0

z_m, z_mand z_m⁰⁰, i.e. to infer the following

ꢄ

3z_m1 + 4ω_r+ 3ω_r

z_r=

ꢅ

2

,

3 (1 + 4ω_m+ 3ω_m²) (z_m− 1) − z_m⁰(4 + 6ω_m) + z_m

00

6z_m(1 + ω_r)

z_Λ= −

,

3 (1 + 3ω_m) (ω_m− ω_r) (1 − z_m) + z_m⁰(1 + 6ω_m− 3ω_r) − z_m

00

(27)

if the denominators of z_rand z_Λare not equal to zero, i.e. if

_m)ξ

z_m= 1+c₁e^{3(1+ω )ξ}+c₂e^(1+3ω

m

and z_m= 1+c₃e^{3(ω −ω )ξ}+c₄e^(1+3ω

,

m

r

(28)

where c₁, c₂, c₃and c₄are arbitrary constants. Since z_m= 0 and ω_r∈ [0, 1],

from (27) follows z_r= 0 and z_Λ= 0. If at least one of the conditions from

(28) is not satisﬁed, from (27) we infer Ω_r= 1/z_r= 0 or Ω_Λ= 1/z_Λ= 0,

which is not consistent with our assumptions in this case. The solutions

for the cases when only one of the density parameters is equal to zero are

analyzed in the Subsect. 3.2. Concretely, the solution for Ω_r= 0 and the

solution for Ω_Λ= 0 are obtained in the equations (45) and (48), respectively.

000

Substituting z_rand z_Λfrom (27) into z_mfrom (23) gives us following

nonhomogeneous linear diﬀerential equation of the third order with constant

coeﬃcients

ꢄ

ꢅ

2

0

z_m⁰⁰⁰− (4 + 9ω_m− 3ω_r) z_m⁰⁰+ 3 1 + 8ω_m+ 9ω_m− 4ω_r− 6ω_mω_rz_m

ꢄ

ꢅ

ꢄ

ꢅ

2

−9 1 + 4ω_m+ 3ω_m

(ω_m− ω_r) z_m+ 9 1 + 4ω_m+ 3ω_m

(ω_m− ω_r) = 0,

whose general solution is

_m−ω_r)ξ

_m)ξ

z_m(ξ) = 1 + c₁e^3(ω

+ c₂e^3(1+ω

+ c₃e^(1+3ω

,

(29)

where c₁, c₂and c₃are strictly positive constants, so that z_m= 1/Ω_m> 1.

Changing z_m⁰(ξ) and z_m⁰⁰(ξ) into (27) gives us z_r(ξ) and z_Λ(ξ). Collecting

all together, solution of the system (22) is

_m−ω_r)ξ

_m)ξ

z_m(ξ) = 1 + c₁e^3(ω

_m)ξ

+ c₂e^3(1+ω

+ c₃e^(1+3ω

,

e^{3(ω −ω}

+ c₁+ c₂e^{3(1+ω )ξ}+ c₃e^{(1+3ω )ξ}

r

z_r(ξ) =

,

(30)

c₁

_m)ξ

e^−3(1+ω

+ c₁e^{−3(1+ω )ξ}+ c₂+ c₃e^−2ξ

r

z_Λ(ξ) =

,

c₂

Dynamics of the ΛCDM model of the universe

14

wherefrom, substituting (21), we infer solutions of the system (20)

1

Ω_m(ξ) =

Ω_r(ξ) =

Ω_Λ(ξ) =

,

_m−ω_r)ξ

_m)ξ

1 + c₁e^3(ω

+ c₂e^3(1+ω

+ c₃e^(1+3ω

c₁

,

(31)

e^{3(ω −ω}

+ c₁+ c₂e^{3(1+ω )ξ}+ c₃e^{(1+3ω )ξ}

_m)ξ

r

c₂

.

_m)ξ

e^−3(1+ω

+ c₁e^{−3(1+ω )ξ}+ c₂+ c₃e^−2ξ

r

It is easy to see from (31) that Ω_i∈ (0, 1). From (31) and (10) we infer the

following dependences

Ω_r(ξ) = c₁e^{3(ω −ω )ξ}Ω_m(ξ),

m

r

Ω_Λ(ξ) = c₂e^{3(1+ω )ξ}Ω_m(ξ),

m

(32)

ꢂ

ꢃ

Ω_k(ξ) = 1 − Ω_m(ξ) 1 + c₁e^3(ω

+ c₂e^3(1+ω

.

_m−ω_r)ξ

_m)ξ

Changing ξ = ln a in (31) and (32), we infer the relations

1

Ω_m(a) =

,

1 + c₁a^{3(ω −ω )}+ c₂a^{3(1+ω )}+ c₃a^1+3ω

m

r

m

c₁

Ω_r(a) =

,

(33)

(34)

a^{3(ω −ω )}+ c₁+ c₂a^{3(1+ω )}+ c₃a^1+3ω

r

m

r

c₂

Ω_Λ(a) =

,

a^{−3(1+ω )}+ c₁a^{−3(1+ω )}+ c₂+ c₃a⁻²

m

r

and

Ω_r(a) = c₁a^{3(ω −ω )}Ω_m(a),

m

r

Ω_Λ(a) = c₂a^{3(1+ω )}Ω_m(a),

m

ꢂ

ꢃ

Ω_k(a) = 1 − Ω_m(a) 1 + c₁a^{3(ω −ω )}+ c₂a^3(1+ω

.

)

m

r

Values of the constants c₁, c₂and c₃we derive using the condition (8)

calculated in the present time moment t = t₀, a(t₀) = 1.

For currently measured mean values of the density parameters Ω_i, we

use notation Ω_i0. According to [14], those values are

Ω_m0= 0.315,

Ω_r0= 0.0000538,

Ω_k0= 0.0007,

Ω_Λ0= 0.685,

(35)

where Ω_m0is the sum of mean values of baryon density and cold dark matter

density of the universe, while for Ω_r0we take mean value of CMB radiation

density of the universe. We note that some values in (35) are calculated

independently, so the identity (8) is not conserved. However, that identity

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

15

must hold in ΛCDM model, therefore here we use normalized values of Ω_i0,

denoted with Ω_i0n, as

Ω_i0

Ω_i0n

=

,

(36)

Ω_m0+ Ω_r0+ Ω_k0+ Ω_Λ0

where Ω_i0are values from (35).

Observing the system (33) in the present time moment t = t₀, a(t₀) = 1,

as well as substituting normalized present mean values of density parame-

ters, we infer the following

c₁≈ 0.000170794,

c₂≈ 2.1746,

c₃≈ 0.00222148.

(37)

Replacing values of c₁, c₂and c₃from (37) into (33), we obtain an approxi-

mate solution of the system (20).

Note that every choice of equation of state parameters ω_m, ω_r∈ [0, 1],

ω_m= ω_r, determine new cosmological model with cosmological constant Λ.

Usually, in cosmology is supposed that ω_m= 0 and ω_r= 1/3. In accor-

dance with that, we substitute those values of the equation of state param-

eters for matter and radiation in system (19) in order to obtain

0

Ω_m= Ω_m(−1 + Ω_m+ 2Ω_r− 2Ω_Λ) ,

0

Ω_r= Ω_r(−2 + Ω_m+ 2Ω_r− 2Ω_Λ) ,

0

Ω_Λ= Ω_Λ(2 + Ω_m+ 2Ω_r− 2Ω_Λ) ,

(38)

Ω_m+ Ω_r+ Ω_Λ+ Ω_k= 1,

˙

Ω_i

0

Ω_i=

.

H

Substituting ω_m= 0 and ω_r= 1/3 in the equations (33) and (34), we obtain

solution of the system (38)

1

Ω_m(a) =

Ω_r(a) =

Ω_Λ(a) =

,

1 + c₁a⁻¹+ c₂a³+ c₃a

c₁

,

(39)

(40)

a + c₁+ c₂a⁴+ c₃a²

c₂

,

a

⁻³+ c₁a⁻⁴+ c₂+ c₃a⁻²

as well as

Ω_r(a) = c₁a⁻¹Ω_m(a),

Ω_Λ(a) = c₂a³Ω_m(a),

ꢄ

ꢅ

Ω_k(a) = 1 − Ω_m(a) 1 + c₁a⁻¹+ c₂a³

.

Dynamics of the ΛCDM model of the universe

16

Furthermore, following the same treatment as in the general case, we obtain

an approximate solution of the system (38), with c₁, c₂and c₃taking values

from (37).

We note that a form of algebraic dependencies (39) and (40) is obtained

in [13], with approach that is independent with the one presented here. In

that way, an algebraic veriﬁcation of this result is inferred in [13].

3.2

Two dimensional dynamical system

Here we consider that only one of four Ω_iis equal to zero. Therefore, we

discuss four diﬀerent cases.

3.2.1 Case Ω_m= 0, Ω_r> 0, Ω_Λ> 0 and Ω_k> 0

This case considers open universe, with both radiation and cosmological

constant. Deﬁciency of matter ﬂuid directly implies to the early universe’s

phase, since a term with properties of dark energy may caused rapid ex-

pansion of the early universe [15], [16]. Nevertheless, that theory is not

consentient with the standard cosmological model. In spite of that, we an-

alyze this case in order to obtain all solutions of the system (19).

In this situation, the system (19) has the following form

0

Ω_r= Ω_r(− (1 + 3ω_r) + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

(41)

0

Ω_Λ= Ω_Λ(2 + Ω_r(1 + 3ω_r) − 2Ω_Λ) ,

˙

Ω_i

0

with Ω_r+Ω_Λ+Ω_k= 1 and Ω_i=

. In the same manner as in Subsect. 3.1,

H

we infer the solution of the system (41)

1

Ω_r(a) =

,

1 + c₁a^{3(1+ω )}+ c₂a^1+3ω

r

(42)

c₁

Ω_Λ(a) =

,

r

a

^{−3(1+ω )}+ c₁+ c₂a⁻²

where c₁and c₂are strictly positive constants. From (42) and Ω_k= 1 −

Ω_r− Ω_Λfollows

Ω_Λ(a) = c₁a^{3(1+ω )}Ω_r(a),

r

ꢂ

ꢃ

(43)

Ω_k(a) = 1 − Ω_r(a) 1 + c₁a^{3(1+ω )}

.

r

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

17

3.2.2

Case Ω_m> 0, Ω_r= 0, Ω_Λ> 0 and Ω_k> 0

According to a lack of radiation, this is an example of an open universe,

around time of density equality of matter and dark energy, i.e. around 10⁹

yrs.

Here, from the system (19) we obtain

0

Ω_m= Ω_m(− (1 + 3ω_m) + Ω_m(1 + 3ω_m) − 2Ω_Λ) ,

(44)

0

Ω_Λ= Ω_Λ(2 + Ω_m(1 + 3ω_m) − 2Ω_Λ) ,

˙

Ω_i

0

with Ω_m+ Ω_Λ+ Ω_k= 1 and Ω_i=

. Following the same treatment as in

H

Subsect. 3.1, we infer the solution of the system (44)

1

Ω_m(a) =

,

1 + c₁a^{3(1+ω )}+ c₂a^1+3ω

m

(45)

c₁

Ω_Λ(a) =

,

a^{−3(1+ω )}+ c₁+ c₂a⁻²

m

where c₁and c₂are strictly positive constants. From (45) and Ω_k= 1 −

Ω_m− Ω_Λfollows

Ω_Λ(a) = c₁a^{3(1+ω )}Ω_m(a),

m

ꢂ

ꢃ

(46)

Ω_k(a) = 1 − Ω_m(a) 1 + c₁a^3(1+ω

.

)

m

Note that, if we replace index m with r in the relations (45) and (46), we

obtain relations (42) and (43). That is because the systems (41) and (44) are

symmetric with respect to indices m and r. Moreover, the same symmetry

holds in the system (19).

3.2.3

Case Ω_m> 0, Ω_r> 0, Ω_Λ= 0 and Ω_k> 0

Here we have an open universe around radiation and matter transition phase,

since at that time dark energy could be neglected.

We consider the system

0

Ω_m= Ω_m(− (1 + 3ω_m) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r)) ,

(47)

0

Ω_r= Ω_r(− (1 + 3ω_r) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r)) ,

Dynamics of the ΛCDM model of the universe

18

˙

Ω_i

0

with Ω_m+Ω_r+Ω_k= 1 and Ω_i=

. In the same manner as in Subsect. 3.1,

H

we obtain the solution of the system (47)

1

Ω_m(a) =

,

1 + c₁a^{3(ω −ω )}+ c₂a^1+3ω

m

r

m

(48)

c₁

Ω_r(a) =

,

a^{3(ω −ω}⁾+ c₁+ c₂a^1+3ω

r

m

r

where c₁and c₂are strictly positive constants. From (48) and Ω_k= 1 −

Ω_m− Ω_rfollows

Ω_r(a) = c₁a^{3(ω −ω )}Ω_m(a),

m

r

ꢂ

ꢃ

(49)

Ω_k(a) = 1 − Ω_m(a) 1 + c₁a^3(ω

.

_m−ω_r)

3.2.4

Case Ω_m> 0, Ω_r> 0, Ω_Λ> 0 and Ω_k= 0

We have a ﬂat universe in the matter dominated stage, more precisely around

5 · 10⁸yrs (radiation and dark energy transition phase).

Substituting Ω_Λ= 1 − Ω_m− Ω_rin (19), we infer

0

Ω_m= Ω_m(−3 (1 + ω_m) + 3 (1 + ω_m) Ω_m+ 3 (1 + ω_r) Ω_r) ,

(50)

0

Ω_r= Ω_r(−3 (1 + ω_r) + 3 (1 + ω_m) Ω_m+ 3 (1 + ω_r) Ω_r) ,

˙

Ω_i

0

with Ω_i=

. Following the procedure presented in Subsect. 3.1, we obtain

H

the solution of the system (50)

1

Ω_m(a) =

,

1 + c₁a^{3(ω −ω )}+ c₂a^3(1+3ω

)

m

r

m

(51)

c₁

Ω_r(a) =

,

a^{3(ω −ω )}+ c₁+ c₂a^{3(1+3ω )}

r

m

r

where c₁and c₂are strictly positive constants. From (51) and Ω_Λ= 1 −

Ω_m− Ω_rfollows

Ω_r(a) = c₁a^{3(ω −ω )}Ω_m(a),

m

r

ꢂ

ꢃ

(52)

Ω_Λ(a) = 1 − Ω_m(a) 1 + c₁a^3(ω

.

_m−ω_r)

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

19

3.3

One dimensional dynamical system

Finally, we observe the system (19) when exactly two of four Ω_iare equal

to zero. That lead to six diﬀerent cases. In each case, the dynamical system

consists of only one diﬀerential equation.

3.3.1

Case Ω_m= 0, Ω_r= 0, Ω_Λ> 0 and Ω_k> 0

Taking into account a lack of matter and radiation, as well as dark energy

domination, this is a future stage of an open universe.

Substituting Ω_m= 0 and Ω_r= 0 into the system (19), we obtain

0

Ω_Λ= Ω_Λ(2 − 2Ω_Λ),

(53)

˙

Ω_Λ

0

with Ω_Λ+ Ω_k= 1 and Ω_Λ

(53) is

=

. The solution of the diﬀerential equation

H

a²

Ω_Λ(a) =

,

(54)

a²+ c₁

where c₁is strictly positive constant. From Ω_k= 1 − Ω_Λand (54), follows

c₁

a²+ c₁

Ω_k(a) =

.

3.3.2

Case Ω_m= 0, Ω_r> 0, Ω_Λ= 0 and Ω_k> 0

In this case, we analyze an open universe with radiation as a dominant ﬂuid,

i.e. an early phase of an open universe.

Replacing Ω_m= 0 and Ω_Λ= 0 into the system (19), we infer

0

Ω_r= Ω_r(− (1 + 3ω_r) + Ω_r(1 + 3ω_r)) ,

(55)

˙

Ω_r

0

with Ω_r+ Ω_k= 1 and Ω_r

derive

=

. Solving (55) and using Ω_k= 1 − Ω_r, we

H

1

Ω_r(a) =

,

1 + c₁a^1+3ω

r

(56)

1

1 + c₁a^1+3ω

Ω_k(a) = 1 −

,

r

where c₁is strictly positive constant.

Dynamics of the ΛCDM model of the universe

20

3.3.3

Case Ω_m= 0, Ω_r> 0, Ω_Λ> 0 and Ω_k= 0

Here we have ﬂat universe with radiation and dark energy. For the same

reason as in 3.2.1, we analyze this case only for the sake of completeness of

this manuscript, despite the absence of physical interpretation in the context

of the ΛCDM model.

Here we observe

0

Ω_r= Ω_r(−3 (1 + ω_r) + 3 (1 + ω_r) Ω_r) ,

(57)

˙

Ω_r

0

with Ω_r+ Ω_Λ= 1 and Ω_r

=

. Solving (57) and replacing its solution

H

into Ω_Λ= 1 − Ω_r, we obtain

1

Ω_r(a) =

,

1 + c₁a^{3(1+ω )}

r

(58)

1

1 + c₁a^{3(1+ω )}

Ω_Λ(a) = 1 −

,

r

where c₁is strictly positive constant.

3.3.4

Case Ω_m> 0, Ω_r= 0, Ω_Λ= 0 and Ω_k> 0

The universe in the description is open and in matter dominated phase.

Since the system (19) is symmetric with respect to indices m and r, the

solution of the diﬀerential equation

0

Ω_m= Ω_m(− (1 + 3ω_m) + Ω_m(1 + 3ω_m)) ,

(59)

˙

Ω_m

0

with Ω_m+ Ω_k= 1 and Ω_m

=

, is the same as the solution of the

H

diﬀerential equation (55), after replacing index r with m, i.e.

1

Ω_m(a) =

,

1 + c₁a^1+3ω

m

(60)

1

Ω_k(a) = 1 −

,

1 + c₁a^1+3ω

m

where c₁is strictly positive constant.

3.3.5

Case Ω_m> 0, Ω_r= 0, Ω_Λ> 0 and Ω_k= 0

Since the presence of both matter and cosmological constant in a ﬂat uni-

verse, as well as the absence of radiation, this suits the present stage of our

universe.

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

21

Again, according to the already mentioned symmetry between the indices

m and r, the solution of the diﬀerential equation

0

Ω_m= Ω_m(−3 (1 + ω_m) + 3 (1 + ω_m) Ω_m) ,

(61)

˙

Ω_m

0

with Ω_m+ Ω_Λ= 1 and Ω_m

=

, is the same as the solution of the

H

diﬀerential equation (57), after replacing index r with m, i.e.

1

Ω_m(a) =

,

1 + c₁a^3(1+ω

)

m

(62)

1

Ω_Λ(a) = 1 −

,

1 + c₁a^3(1+ω

)

m

where c₁is strictly positive constant.

3.3.6

Case Ω_m> 0, Ω_r> 0, Ω_Λ= 0 and Ω_k= 0

In this case, we deal with a ﬂat universe around a time period of transition

between radiation and matter.

Here, we solve

0

Ω_m= Ω_m(−3 (ω_m− ω_r) + 3 (ω_m− ω_r) Ω_m) ,

(63)

˙

Ω_m

0

with Ω_m+ Ω_r= 1 and Ω_m

=

. Solving previous diﬀerential equation

H

and substituting its solution into Ω_r= 1 − Ω_m, we infer

a^3ω

r

Ω_m(a) =

,

a^3ω+ c₁a^3ω

r

m

(64)

a^3ω

r

Ω_r(a) = 1 −

,

a^3ω+ c₁a^3ω

r

m

where c₁is strictly positive constant.

In various literature (for example see [2], [4], [7], [8]) the system (38) or

some similar systems are analyzed by applying theory of dynamical systems.

We present that method here in more details in order to analyze the system

(19), which is more general. That analysis enable us to clearly understand

the evolution of the universe from new perspective. In the next section we

will see how density parameters ruled the expansion of the universe from its

beginning and how their behaviour will have direct impact on the universe’s

future.

Dynamics of the ΛCDM model of the universe

22

4 Evolution of the dynamical system

According to Ω_i∈ [0, 1] and (8), the phase space for the system (19) is

tetrahedron whose vertices are points E₀(0, 0, 0), E₁(1, 0, 0), E₂(0, 1, 0) and

E₃(0, 0, 1), which lie on the axes Ω_m, Ω_rand Ω_Λ, (Figure 1, see [17]).

Since ω_m, ω_r∈ [0, 1] and ω_m= ω_r, the only equilibriums of the sys-

tem (19) are exactly the vertices E₀(0, 0, 0), E₁(1, 0, 0), E₂(0, 1, 0) and

E₃(0, 0, 1).

Figure 1: Phase space portrait of the dynamical system (19) with the values

ω_m= 0 and ω_r= 1/3.

As we will see in the following text, the condition ω_m= ω_rensure that

the eigenvalues of the Jacobian matrix in our equilibrium points E_i, where

i ∈ {0, 1, 2, 3}, have non-zero real parts. Thus, our equilibriums E_iare hy-

perbolic so we can analyze their stability in the framework of the correspond-

ing linearization of the dynamical system (19), since the Hartman–Grobman

theorem gives us topological equivalence between a nonlinear dynamical sys-

tem and its linearization in the neighbourhood of the hyperbolic equilibrium.

Now we inspect the stability of all equilibriums of the system (19).

The equilibrium E₀. The eigenvalues of the Jacobian matrix at the

equilibrium E₀are λ₁= − (1 + 3ω_m), λ₁< 0, λ₂= − (1 + 3ω_r), λ₂< 0

and λ₃= 2, wherefrom we conclude that E₀is a saddle point (unstable

equilibrium). Considering that E₀is the origin, i.e. that Ω_m= 0, Ω_r= 0

and Ω_Λ= 0, we infer Ω_k= 1, wherefrom k < 0. Hence, the universe near

the equilibrium E₀is empty and open universe without Λ, known as the

Milne universe. According to the discussion above, we conclude that this is

an unstable phase of the universe’s evolution.

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

23

The equilibrium E₁. When it comes to the equilibrium E₁, the

eigenvalues of the Jacobian matrix at E₁are λ₁= 1 + 3ω_m, λ₁> 0,

λ₂= 3 (ω_m− ω_r), λ₂= 0 and λ₃= 3 (1 + ω_m), λ₃> 0. We have two

cases, depending on the sign of (ω_m− ω_r). If ω_m> ω_r, then λ₂> 0, which

implies that E₁is an unstable node (unstable equilibrium). If ω_m< ω_r, then

E₁is a saddle point (unstable equilibrium). According to the coordinates

of E₁, i.e. that Ω_m= 1, Ω_r= 0 and Ω_Λ= 0, we infer Ω_k= 0, wherefrom

k = 0. Therefore, the universe near the equilibrium E₁is ﬂat universe and

matter dominated, without Λ, known as the Einstein-de Sitter universe. We

note that, regardless of the sign of (ω_m− ω_r), this phase of the universe’s

evolution is unstable.

The equilibrium E₂. The eigenvalues of the Jacobian matrix at the

equilibrium E₂are symmetric to the eigenvalues at E₁, with respect to the

indices m and r, i.e. λ₁= 1 + 3ω_r, λ₁> 0, λ₂= 3 (ω_r− ω_m), λ₂= 0

and λ₃= 3 (1 + ω_r), λ₃> 0. Again, the nature of the equilibrium depends

on the sign of (ω_m− ω_r). Here, if ω_m> ω_r, then λ₂< 0, which indicates

that E₂is a saddle point (unstable equilibrium). If ω_m< ω_r, then E₂is

an unstable node (unstable equilibrium). Since in this case is Ω_m= 0,

Ω_r= 1 and Ω_Λ= 0, we infer Ω_k= 0, as well as k = 0. To conclude, the

universe near the equilibrium E₂is ﬂat universe and radiation dominated,

without cosmological constant. As for the previous case, the universe near

this equilibrium is unstable.

The equilibrium E₃. Finally, the eigenvalues of the Jacobian matrix

at the equilibrium E₃are λ₁= −3 (1 + ω_m), λ₁< 0, λ₂= −3 (1 + ω_r),

λ₂< 0 and λ₃= −2, wherefrom E₃is a stable node (asymptotically stable

equilibrium). Here we have Ω_m= 0, Ω_r= 0, Ω_Λ= 1 and again k = 0.

Hence, the universe near the equilibrium E₃is ﬂat and Λ-dominated, also

known as the de-Sitter universe. Thus, this is a stable phase in the evolution

of the universe.

To summarize, E₀is a saddle point (unstable equilibrium), E₃is a stable

node (asymptotically stable equilibrium), while the nature of E₁and E₂

depends on the sign of (ω_m− ω_r). Without loss of generality, we analyze

the case ω_m< ω_r. Thus, E₁is a saddle point (unstable equilibrium), whereas

E₂is an unstable node (unstable equilibrium).

Solutions of the system (19) are the orbits that start at E₂as ξ → −∞,

i.e. a(t) → 0₊, and end in E₃as ξ → +∞, i.e. a(t) → +∞. It means

that the early universe was ﬂat and radiation dominated, while the future of

the ﬂat universe is characterised by Λ-domination, which is consistent with

the current observations. The exceptions are the orbits on the boundaries

Dynamics of the ΛCDM model of the universe

24

and in the interior of the triangles ∆E₀E₁E₂and ∆E₀E₁E₃, that lie in

the (Ω_m,Ω_r)-plane and (Ω_m,Ω_Λ)-plane, respectively. Those exceptions are

expected, since they correspond to the trajectories in a universe without

cosmological constant, or to the trajectories in a universe without radiation,

respectively.

4.1

A universe without cosmological constant Λ

Replacing Ω_Λ= 0 into (19), we obtain

0

Ω_m= Ω_m(− (1 + 3ω_m) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r)) ,

(65)

0

Ω_r= Ω_r(− (1 + 3ω_r) + Ω_m(1 + 3ω_m) + Ω_r(1 + 3ω_r)) ,

˙

Ω_i

0

with Ω_i∈ [0, 1], Ω_m+ Ω_r+ Ω_k= 1 and Ω_i=

, which is the same system

H

as (47). The solutions of the system (65) are derived in Subsubsects. 3.2.3,

3.3.2, 3.3.4 and 3.3.6, depending on whether all Ω_i> 0, or otherwise.

The phase space for the system (65) is the triangle in the (Ω_m, Ω_r)-

plane whose vertices are points with coordinates (0, 0), (1, 0) and (0, 1),

which are exactly the coordinates of the points E₀, E₁and E₂in the plane

Ω_Λ= 0. Therefore, for the sake of simplicity we denote our phase space

with ∆E₀E₁E₂(Figure 2, see [17]).

1.0

E₂

0.8

0.6

0.4

0.2

E₀

E₁

1.0

0.0

0.2

0.4

0.6

0.8

Figure 2: Phase space portrait of the dynamical system (65) with the values

ω_m= 0 and ω_r= 1/3.

According to ω_m, ω_r∈ [0, 1] and ω_m< ω_r, the equilibriums of the system

(65) are vertices E₀(0, 0), E₁(1, 0), and E₂(0, 1). Utilizing linear stability

analysis near the equilibrium points, we infer the following.

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

25

The equilibrium E₀. The eigenvalues of the Jacobian matrix at E₀

are λ₁= − (1 + 3ω_m), λ₁< 0 and λ₂= − (1 + 3ω_r), λ₂< 0, wherefrom

we conclude that E₀is a stable node (asymptotically stable equilibrium).

Similarly as in the three dimensional case, we conclude that the universe

near the equilibrium E₀is the Milne universe (empty and open universe

without Λ).

The equilibrium E₁. In this case, the eigenvalues of the Jacobian

matrix at E₁are λ₁= 1 + 3ω_m, λ₁> 0 and λ₂= 3 (ω_m− ω_r), λ₂< 0.

Consequently, E₁is a saddle point (unstable equilibrium). Furthermore, a

universe near the equilibrium E₁is Einstein-de Sitter universe (ﬂat universe

and matter dominated, without Λ).

The equilibrium E₂. Finally, E₂is an unstable node (unstable equi-

librium), given that the eigenvalues of the Jacobian matrix at E₂are λ₁=

1 + 3ω_r, λ₁> 0 and λ₂= 3 (ω_r− ω_m), λ₂> 0. The universe near the equi-

librium E₂is ﬂat universe and radiation dominated, without cosmological

constant.

In conclusion, solutions of the system (65) are the orbits that start at E₂

as ξ → −∞ and end in E₀as ξ → +∞. That is quite interesting, because

it states that in the absence of Λ, nearly ﬂat and radiation dominated uni-

verse will tend to evolve into the Milne universe, which is open and empty

universe. The exceptions are the orbits that connect E₂to E₁and E₁to

E₀, i.e. the ones on the line Ω_r= 1 − Ω_mand on the Ω_m-axis. These or-

bits correlate either to a ﬂat universe without Λ, or to a universe without

radiation and Λ. In the ﬁrst case, they describe the tendency of evolving

ﬂat and radiation dominated universe without cosmological constant into

the Einstein-de Sitter universe, or evolving the Einstein-de Sitter universe

into the Milne universe in the latter case.

4.2

A universe without radiation

Substituting Ω_r= 0 into (19), we infer

0

Ω_m= Ω_m(− (1 + 3ω_m) + Ω_m(1 + 3ω_m) − 2Ω_Λ) ,

(66)

0

Ω_Λ= Ω_Λ(2 + Ω_m(1 + 3ω_m) − 2Ω_Λ) ,

˙

Ω_i

0

with Ω_m+ Ω_Λ+ Ω_k= 1 and Ω_i=

, which is the same system as (44).

H

The solutions of the system (66) are obtained in Subsubsects. 3.2.2, 3.3.1,

3.3.4 and 3.3.5, depending on whether all Ω_i> 0, or otherwise.

Dynamics of the ΛCDM model of the universe

26

The phase space for the system (66) is the triangle in the (Ω_m, Ω_Λ)-plane

whose vertices are points with coordinates (0, 0), (1, 0) and (0, 1), which are

exactly the coordinates of the points E₀, E₁and E₃in the plane Ω_r= 0.

Thus, for the sake of simplicity we denote our phase space with ∆E₀E₁E₃

(Figure 3, see [17]).

1.0

E₃

0.8

0.6

0.4

0.2

E₀

E₁

1.0

0.0

0.2

0.4

0.6

0.8

Figure 3: Phase space portrait of the dynamical system (66) with the value

ω_m= 0.

Since ω_m∈ [0, 1], the equilibriums of the system (66) are the vertices

E₀(0, 0), E₁(1, 0), and E₃(0, 1). Using the Hartman-Grobman theorem,

we infer the following.

The equilibrium E₀. In this case, the eigenvalues of the Jacobian

matrix at E₀are λ₁= − (1 + 3ω_m), λ₁< 0 and λ₂= 2. Consequently, E₀is

a saddle point (unstable equilibrium). In the same manner as in the three

dimensional case, we conclude that the universe near the equilibrium E₀is

the Milne universe (empty and open universe without Λ).

The equilibrium E₁. The equilibrium E₁is an unstable node (unstable

equilibrium), given that the eigenvalues of the Jacobian matrix at E₁are

λ₁= 1 + 3ω_m, λ₁> 0 and λ₂= 3 (1 + ω_m), λ₂> 0. The universe near the

equilibrium E₁is the Einstein-de Sitter universe (ﬂat universe and matter

dominated, without Λ).

The equilibrium E₃. The eigenvalues of the Jacobian matrix at E₃are

λ₁= −2 and λ₂= −3 (1 + ω_m), λ₂< 0, wherefrom we conclude that E₃is a

stable node (asymptotically stable equilibrium). Moreover, the universe near

the equilibrium E₃is the de-Sitter universe (ﬂat universe and Λ dominated).

Solutions of the system (66) are the orbits that start at E₁as ξ → −∞

and end in E₃as ξ → +∞. It implies that the Einstein-de Sitter universe

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

27

will tend to evolve to the de-Sitter universe. The exceptions are the orbits

that connect E₁to E₀and E₀to E₃, i.e. the ones on the Ω_m-axis and on

the Ω_Λ-axis. These orbits correlate either to a universe without radiation

and Λ, or to a universe without radiation and matter. In the ﬁrst case,

they describe the tendency of evolving the Einstein-de Sitter universe into

the Milne universe (we noticed this in Subsect. 4.1), or evolving the Milne

universe into the de-Sitter universe in the latter case.

5 Stability analysis of equilibriums when ω_m= ω_r

If we let ω_m= ω_r, then E₁is not hyperbolic equilibrium (λ₂= 3 (ω_m− ω_r) =

0). Therefore we cannot determine its stability using the Hartman-Grobman

theorem. Instead, we will determine its stability using the Lyapunov stabil-

ity theorem. Concretely, we will show that the equilibrium E₁is unstable

on (0, +∞)³. Consequently, the equilibrium E₁would also be unstable on

the tetrahedron on Figure 1, i.e. on the phase space of the system (19).

Let suppose that the Lyapunov function is in the form V₁(Ω_m, Ω_r, Ω_Λ) :=

Ω_m⁻³Ω_rΩ_Λ(see [18]). Then V₁(Ω_m, Ω_r, Ω_Λ) = 3 (1 + ω_m) Ω_m⁻³Ω_rΩ_Λ.

We see that function V₁is diﬀerentiable on (0, +∞)³, V₁(E₁) = 0 and

V₁(Ω_m, Ω_r, Ω_Λ) > 0 for every (Ω_m, Ω_r, Ω_Λ) ∈ (0, +∞)³. Since ω_m∈ [0, 1],

2

˙

3

˙

then V₁(Ω_m, Ω_r, Ω_Λ) > 0 for every (Ω_m, Ω_r, Ω_Λ) ∈ (0, +∞) . Consequently,

the function V₁is a Lyapunov function for the equilibrium E₁and E₁is an

unstable equilibrium on (0, +∞)³.

Moreover, assuming the condition ω_m= ω_r, the same result holds for

the equilibrium E₂, i.e. the non-hyperbolic equilibrium E₂is unstable on

(0, +∞)³and consequently, it is unstable on the appropriate tetrahedron. In

order to see that, we use the Lyapunov function in the form V₂(Ω_m, Ω_r, Ω_Λ) :=

Ω_m²Ω_r⁻³Ω_Λ.

To summarize, if ω_m= ω_rthen both equilibriums E₁and E₂are unstable

on the phase space of the system (19).

6 Conclusion

In this paper we discussed the dynamics of the universe in the framework

of the ΛCDM model. We established the equivalence between the system

of the Friedmann equations (1) and the dynamical system (19), under some

assumptions that are standard in cosmology. Furthermore, all analytical

solutions of the system (19) are derived and represent new parameterizations

Dynamics of the ΛCDM model of the universe

28

of the density parameters with respect to the scale factor. All solutions of the

system (19) are connected with appropriate stage of the universe’s evolution.

We demonstrated that linear stability analysis near the equilibrium points

indicates some interesting behavior with respect to the dynamics of the

universe. Moreover, in order to obtain stability analysis of equilibriums

when ω_m= ω_r, we found the appropriate Lyapunov functions.

Note. All computations are checked using the Wolfram Mathematica

package.

Acknowledgements. The authors would like to thank Jelena Kati´c

and Jovana Nikoli´c for careful reading of the paper and the constructive

and helpful remarks.

References

¨

[1] A. Friedmann, Uber die M¨oglichkeit einer Welt mit konstanter negativer

Kru¨mmung des Raumes, Z. Phys. 21 (1924), 326-332.

[2] S. Bahamonde, C.G. B¨ohmer, S. Carloni, E.J. Copeland, W. Fang and

N. Tamanini, Dynamical systems applied to cosmology: Dark energy

and modiﬁed gravity, Phys. Rep. 775–777 (2018), 1-122.

[3] M. Hobson, G. Efstathiou and A. Lasenby, General Relativity: An

Introduction for Physicists, Cambridge University Press, Cambridge,

2006.

[4] J. Perez, A. Fu¨zfa, T. Carletti, L. M´elot and S. Guedezounme, The Jun-

gle Universe: Coupled cosmological models in a Lotka-Volterra frame-

work, Gen. Relativ. Gravitation 46 (2014), 1753.

[5] A.R. Liddle and D.H. Lyth, Cosmological Inﬂation and Large-Scale

Structure, Cambridge University Press, 2000.

[6] D. Brankovi´c, Friedmann equations as n-dimensional dynamical sys-

tem, An. St. Univ. Ovidius Constanta 31 (2023), 23-37.

[7] R. Garc´ıa-Salcedo, T. Gonzales, F.A. Horta-Rangel, I. Quiros and D.

Sanchez-Guzm´an, Introduction to the application of the dynamical sys-

tems theory in the study of the dynamics of cosmological models of dark

energy, Eur. J. Phys. 36 (2015), 025008.

ˇ

D. Brankovi´c, Z. Mijajlovi´c, M. Miki´c

29

[8] J.-P. Uzan and R. Lehoucq, A dynamical study of the Friedmann equa-

tions, Eur. J. Phys. 22 (2001), 371.

[9] M. Goliath and G.F.R. Ellis, Homogeneous cosmologies with a cosmo-

logical constant, Phys. Rev. D 60 (1999), 023502.

[10] A.J. Lotka, Elements of Physical Biology, Williams-Wilkins, Baltimore,

1925.

[11] V. Volterra, Le¸con sur la Theorie Math´ematique de la Lutte pour la

Vie, Gauthier-Villars, Paris, 1931.

[12] D. Brankovi´c, Cosmic time for multi-component universe, Serb. Astron.

J. 201 (2020), 15-23.

ˇ

[13] Z. Mijajlovi´c and D. Brankovi´c, On Algebraic Dependence of Cosmo-

logical Parameters, Gravit. Cosmol. 29 (2023), 456–467.

[14] S. Navas et al. (Particle Data Group), Phys. Rev. D 110 (2024), 030001.

[15] F. Niedermann and M.S. Sloth, New early dark energy, Phys. Rev. D

103 (2021), L041303.

[16] V. Poulin, T.L. Smith, T. Karwal and M. Kamionkowski, Early Dark

Energy can Resolve the Hubble Tension, Phys. Rev. Lett. 122 (2019),

221301.

[17] Mathematica ﬁles 1712.03107 [gr-qc],

https://www.christianboehmer.co.uk/review.html (accessed 29

April 2025).

[18] M.L. Zeeman, Extinction in competitive Lotka-Volterra systems, Proc.

Amer. Math. Soc. 123 (1995), 87-96.