EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 3, Article Number 6581
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Asymptotic Analysis of Varying Light Speed
Cosmologies

Dimitrios Trachilis

College of Engineering and Technology, American University of the Middle East,
Egaila 54200, Kuwait

Abstract. We investigate potential singularities in isotropic cosmological models featuring a time-
varying speed of light and the gravitational constant. The governing field equations generally
simplify into two-dimensional systems, subject to analysis through dynamical systems techniques
within phase space and the application of asymptotic splitting methods. In the broader context,
our findings reveal the existence of initially expanding closed models that undergo subsequent
recollapse towards a future singularity, and open universes exhibiting perpetual expansion into the
future. The specific characteristics of these singularities are also expounded upon.

2020 Mathematics Subject Classifications: 83F05, 34C20, 83C05, 34A34

Key Words and Phrases: Method of asymptotic splittings, Varying light speed cosmologies,
Asymptotic solutions

1. Introduction

Varying Speed of Light (VSL) cosmologies have garnered significant attention as al-
ternatives to cosmological inflation, offering an alternative framework for addressing the
challenges posed by the standard model (refer to [1] and associated references for further
details). Rather than embracing the inflationary paradigm positing a phase of superlu-
minal expansion in the early universe, these cosmologies postulate an accelerated speed
of light during this epoch. In such cosmological frameworks, the conventional enigmas
encountered in standard early universe cosmology are circumvented, albeit at the expense
of violating the general covariance of the underlying gravitational theory. This response
refrains from delving into a discourse on the foundational aspects of VSL theories and
the conceptual intricacies stemming from the variable nature of light speed (see [1] for an
in-depth discussion).

In this paper, we present initial findings from [2] involving the application of the asymp-
totic splittings method, as introduced in [3], to analyze singularities that may manifest
in Variable Speed of Light (VSL) cosmologies. Specifically, our focus is directed towards

DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6581

Email address: Dimitrios.Trachilis@aum.edu.kw (D. Trachilis)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 2 of 12

the VSL model proposed in [4] and further elucidated in [5]. A noteworthy feature of
this model lies in the assumption of minimal coupling within the framework of Einstein’s
equations. This assumption posits that a time-variable speed of light (c) should not in-
duce alterations in the curvature terms within Einstein’s equations in the cosmological
context; thus, the foundational equations of Einstein persist. As noted by Barrow in [6],
changes in the speed of light occur within the local Lorentzian frames associated with cos-
mological expansion, representing a special-relativistic effect. Consequently, the resulting
theory lacks covariance, necessitating the adoption of a specific time coordinate. Opting
for comoving proper time as this specific choice, the Friedman equations maintain their
structure while accommodating variations in both c(t) and G(t).

The plan of this paper is as follows. In the next two sections, we use the Friedman
equations with a power representation for the speed of light to derive an autonomous dy-
namical system in suitable variables. Then, we apply the method of asymptotic splittings
to study the evolution of this model. We show that there exist asymptotic solutions of
the unknown variables near the initial singularity. In sections 4-5, apart from the speed of
light, we consider an additional power representation of G and apply the same asymptotic
method to find the admissible asymptotic solutions of this generalized model. In section
6, we connect the general solutions found from the method of asymptotic splittings to the
flatness problem described in [7] and [8].

We start with the Friedman equations,

(
ȧ

a
)2 = −kc2

a2
+

8πG

3
ρ, (1)

ä

a
= −8πG

6
(3γ − 2)ρ, (2)

where the equation of state has the form,

p = (γ − 1)ρc2, (3)

with 0 < γ ≤ 2. Here a = a(t) is the scale factor, p is the fluid pressure, ρ is the fluid
density, and k is the constant curvature normalized to k = 0,+1 or −1, for the flat, closed,
or open space, respectively.

We consider the basic Friedman cosmological equations (1) and (2) for varying speed
of light c and for both constant and varying Newtonian gravitational G to derive the asso-
ciated conservation law of the theory. We show that these equations form an autonomous
dynamical system in suitable variables. Further, we apply the method of asymptotic split-
tings given in [3] to examine the early evolution of the varying speed of light cosmology.
In particular, we study separately 2+1 cases with respect to the varying c and G, namely
that the subcases of varying c with constant G, varying c and G, and constant c with
varying G. The last subcase seems not to belong to a varying speed of light cosmology,
but in fact it provides an equivalent problem. For every previous representation of c and
G we find valid asymptotic solutions near the initial singularity t0.



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 3 of 12

2. Reduction to two dimensions with constant G

We assume that the varying speed of light c has a representation of the form,

c = c0a
n, n ∈ R, (4)

the system of units of which satisfies 8πG = 1 = c0
2. Differentiating eq. (1) and using eq.

(2), we find the conservation law,

ρ̇+ 3γρ
ȧ

a
=

3k

4πG

cċ

a2
. (5)

Obviously, if ċ = 0, that is n = 0, we obtain the same equations as in general relativity
because the right-hand side of (5) vanishes. Setting

x =
1

a
, (6)

H =
ȧ

a
, (7)

the system (1), (2), and (5) implies,

ẋ = −xH, (8)

Ḣ = −H2 − 3γ − 2

6
ρ, (9)

ρ̇ = −3γρH + 6knx2−2nH, (10)

(11)

together with the constraint,

H2 + kx2−2n =
1

3
ρ. (12)

Using the last equation, we can eliminate x or ρ. If we eliminate x, the remaining 2-
dimensional system takes the form,

Ḣ = −H2 − 3γ − 2

6
ρ,

ρ̇ = −(3γ − 2n)ρH − 6nH3. (13)

3. Asymptotic forms of H and ρ with constant G

In this section we will apply the method of asymptotic splittings [3] (see also [9],[10]) to
the system (13) in order to find the dominant features on approach to the initial singularity,
taken at t = 0. In fact, the singularity can be replaced by any t = t0 using the variable
τ = t − t0. System (13) can be written as an autonomous system of the form, Below
we are interested in finding asymptotic solutions that correspond to varying light speed
cosmologies, thus we consider only the n ̸= 0 cases. We start with,

ẋ = f(x), x = (H, ρ), (14)



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 4 of 12

and can be decomposed asymptotically into a dominant part and another, subdominant
part

f(x) = f
(0)
i (x) + f

(sub)
i (x), (15)

in 9 different ways (i = 1, 2, 3, · · · , 9). A solution to the system (13), which is called
dominant near the singularity, is of the form,

x(t) = atp = (αtp, βtq), (16)

where a = (α, β) ∈ C2 and p = (p, q) ∈ Q2. The pair (a,p) is called a dominant balance
of the vector field f(x) if it corresponds to a dominant solution of the system (13) near
the singularity.

The first possible asymptotic decomposition is

f(x) = f
(0)
1 (x) + f

(sub)
1 (x), (17)

with dominant part,

f
(0)
1 (x) =

(
−H2,−(3γ − 2n)ρH

)
, (18)

and subdominant part,

f
(sub)
1 (x) =

(
− 3γ − 2

6
ρ,−6nH3

)
, (19)

Substituting
H = αtp and ρ = βtq (20)

into the dominant system ẋ = f
(0)
1 (x) we evaluate the dominant balance (a,p). In partic-

ular, we find that,

(a,p) =
(
(1, 0), (−1, q)

)
, q ≤ −1. (21)

Further, we need to determine whether the term (19) is indeed subdominant by evaluating
the limit of this term divided by tp−1 when t → 0. For values of q ≤ 2, we find,

f
(sub)
1 (x)

tp−1
→ (0, 0), when t → 0, (22)

and therefore the decomposition (17) is acceptable asymptotically.
Next, we evaluate the Kovalevskaya matrix, given by,

K = Df
(0 )
1 (a)− diag(p), (23)

and we find that there is always one eigenvalue λ = −1 that corresponds to an arbitrary
constant of the dominant behaviour of the solution given by the Fuchsian series expansions,

H =

∞∑
j=0

cj1t
p−j , ρ =

∞∑
j=0

cj2t
q−j , (24)



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 5 of 12

where c01 = α = 1 and c02 = β = 0. Substituting the series (24) into the original system
(13), we find that there is no solution for H and ρ unless n = 0. This is because of our
final test for admission of the solutions, namely that, the Fredholm’s alternative,

vT (K − jI)cj = 0, (25)

where v denote the eigenvectors of K and I is the identity matrix. In particular, in
our first decomposition, the expression (25) is true only for n = 0. Therefore, the first
decomposition does not give any valid asymptotic solution of the initial system around
the singularity.

The second asymptotic decomposition is

f(x) = f
(0)
2 (x) + f

(sub)
2 (x), (26)

with dominant part,

f
(0)
2 (x) =

(
−H2,−6nH3

)
, (27)

and subdominant part,

f
(sub)
2 (x) =

(
− 3γ − 2

6
ρ,−(3γ − 2n)ρH

)
, (28)

Substituting H = αtp and ρ = βtq into the second dominant system, we find the dominant
balance,

(a,p) =
(
(1, 3n), (−1,−2)

)
. (29)

The corresponding subdominant term (28) satisfies

f
(sub)
2 (x)

tp−1
=

(2− 3γ

2
n, 3n(2n− 3γ)

)
, (30)

and therefore approaches (0, 0) if and only if (γ = 2
3 , n = 1) or (γ ∈ R, n = 0). In principle,

we can accept only the first pair. However, the expression (25) does not give any acceptable
solution for the combination γ = 2

3 and n = 1. Thus the second decomposition does not
lead to a fully acceptable solution as well.

The third possible asymptotic decomposition is

f(x) = f
(0)
3 (x) + f

(sub)
3 (x). (31)

The dominant part is,

f
(0)
3 (x) =

(
−H2,−(3γ − 2n)ρH − 6nH3

)
, (32)

and the subdominant part is,

f
(sub)
3 (x) =

(
− 3γ − 2

6
ρ, 0

)
. (33)



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 6 of 12

Substituting (20) into the third dominant system, we find,

(a,p) =
(
(1, 0), (−1,−2)

)
and (a,p) =

(
(1, 3), (−1,−2)

)
. (34)

In case that (a,p) =
(
(1, 0), (−1,−2)

)
, we find no asymptotic solutions with n ̸= 0, but

in case that (a,p) =
(
(1, 3), (−1,−2)

)
, we obtain the special solution,

H = t−1, ρ = 3t−2, where k = 0, γ =
2

3
, n = −1

2
. (35)

In addition, the corresponding subdominant term (33) satisfies

f
(sub)
3 (x)

tp−1
→ (0, 0), (36)

when t → 0, in this case.
The next three possible asymptotic decompositions,

f(x) = f
(0)
j (x) + f

(sub)
j (x), j = 4, 5, 6, (37)

for which,

f
(0)
4 (x) =

(
− 3γ − 2

6
ρ,−(3γ − 2n)ρH

)
, (38)

f
(0)
5 (x) =

(
− 3γ − 2

6
ρ,−6nH3

)
, (39)

f
(0)
6 (x) =

(
− 3γ − 2

6
ρ,−(3γ − 2n)ρH − 6nH3

)
, (40)

with corresponding subdominant parts,

f
(sub)
4 (x) =

(
−H2,−6nH3

)
, (41)

f
(sub)
5 (x) =

(
−H2,−(3γ − 2n)ρH

)
, (42)

f
(sub)
6 (x) =

(
−H2, 0

)
, (43)

give no asymptotic solution due to the fact that (41)-(43) do not satisfy the conditions

f
(sub)
j (x)

tp−1
→ (0, 0), j = 4, 5, 6, when t → 0. (44)

We apply the method of asymptotic splittings for the remaining three asymptotic
decompositions, the last of which is the all-terms dominant decomposition. We have,

f(x) = f
(0)
j (x) + f

(sub)
j (x), j = 7, 8, 9, (45)



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 7 of 12

where,

f
(0)
7 (x) =

(
−H2 − 3γ − 2

6
ρ,−(3γ − 2n)ρH

)
, (46)

f
(0)
8 (x) =

(
−H2 − 3γ − 2

6
ρ,−6nH3

)
, (47)

f
(0)
9 (x) =

(
−H2 − 3γ − 2

6
ρ,−(3γ − 2n)ρH − 6nH3

)
, (48)

and,

f
(sub)
7 (x) =

(
0,−6nH3

)
, (49)

f
(sub)
8 (x) =

(
0,−(3γ − 2n)ρH

)
, (50)

f
(sub)
9 (x) =

(
0, 0

)
, (51)

In every previous decomposition (46)-(48), we find that p = −1 and q = −2. However,
the possible values of α and β that are accepted from the decomposed systems vary, and
give the following six special asymptotic solutions for H and ρ in (13),

H =
1

1− c
t−1, ρ =

3c

(1− c)2
t−2, where k ̸= 0, γ = 0, and n = c satisfies

−18c3 + 89c2 − 106c− 1 = 0, (52)

H =
c− 1

108c
t−1, ρ =

(c− 1)3

419904c2
t−2, where (k = 0, n = c = 109, γ =

218

3
), or

(k < 0, n = c = −53± 6
√
78, γ =

2c

3
), (53)

H =
1

1− c
t−1, ρ =

6c

(1− c)2(2− 3γ)
t−2, where k = 0,±1, γ ̸= 0,

2

3
, n = c ̸= 1,

and (3γ − 2)2(c2 + 34c+ 1) + 36c(c− 1)(9γ + 2c− 8) = 0, (54)

H =
1

1− c
t−1, ρ =

3

(1− c)2
t−2, where k = 0, γ =

2− 2n

3
, n = c ̸= 1,

and c3 − 2c2 + 19c+ 18 = 0, (55)

H =
2

3γ
t−1, ρ =

4

3γ2
t−2, where k = 0, γ ̸= 0,

2

3
, n ̸= 1,

and γ2(3γ − 2)(3γ − 7) + 8(3γ + 4n)(3γ − 2n) = 0, (56)

and the general asymptotic solution,

H = t−1 + c11t
−2 + · · · , ρ = 3t−2 + c12t

−3 + · · · , where c12 = 3c11,

and k = 0,±1, γ =
2

3
, n = −1

2
. (57)

Solution (57) contains two arbitrary constants, c11 and one more that corresponds to the
arbitrary position of the singularity (here without loss of generality, we consider t0 = 0).



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 8 of 12

Notice that the special solution (35) found in the third decomposition is just a special case
of (57). Our various series solutions found in all nine asymptotic decompositions satisfy the
phase space of the initial system (13) that follows from Eq. (12). In particular, for closed
models, the phase portrait of (13) is given by the set

{
(H, ρ) ∈ R2 : ρ ≥ 0, ρ− 3H2 > 0

}
,

and for open models it is given by
{
(H, ρ) ∈ R2 : ρ ≥ 0, ρ− 3H2 < 0

}
.

4. Reduction to two dimensions with varying G

In this section we assume that both functions c = c(t) and G = G(t) have a power-law
dependence of the form,

c = c0a
n, G = G0a

m, n,m ∈ R, (58)

where the system of units requires 8πG0 = c0
2 = 1. In this general case, if we differentiate

(1) and use (2), we obtain the generalized conservation law,

ρ̇+ 3γρ
ȧ

a
= −ρ

Ġ

G
+

3k

4πG

cċ

a2
. (59)

Setting again x = 1
a and H = ȧ

a , our initial system (1), (2), together with (59), becomes,

ẋ = −xH, (60)

Ḣ = −H2 − 3γ − 2

6
ρx−m, (61)

ρ̇ = −(m+ 3γ)ρH + 6knxm+2−2nH, , (62)

subject to the constraint,

H2 + kx2−2n =
1

3
ρx−m. (63)

Using (63) we can only eliminate ρ and get a two-dimensional system with respect to
(x,H). The system that we find is,

ẋ = −xH,

Ḣ = −3γ

2
H2 − 3γ − 2

2
kx2−2n. (64)

Notice that the value of m is not present in (64). However, it must be present in the
solution of ρ due to the constraint (63).

5. Asymptotic forms with varying G

In this section we will apply the method of asymptotic splittings to the system (64)
generated by the Friedman equations with a power-law representation of c and G. There
are three different decomposed systems to apply the method of asymptotic splittings, that
is,

fgen(x) = f
(0)
gen,j(x) + f

(sub)
gen,j (x), j = 1, 2, 3, (65)



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 9 of 12

where,

f
(0)
gen,1(x) =

(
− xH,−3γ

2
H2

)
, (66)

f
(0)
gen,2(x) =

(
− xH,−3γ − 2

2
kx2−2n

)
, (67)

f
(0)
gen,3(x) =

(
− xH,−3γ

2
H2 − 3γ − 2

2
kx2−2n

)
, (68)

and,

f
(sub)
gen,1(x) =

(
0,−3γ − 2

2
kx2−2n

)
, (69)

f
(sub)
gen,2(x) =

(
0,−3γ

2
H2

)
, (70)

f
(sub)
gen,3(x) =

(
0, 0

)
, (71)

The second decomposition (67) does not give any asymptotic solution to the problem,
because the compatibility condition (25) is not valid in this case. However, the first de-
composition and the all-terms dominant case provide the following admissible asymptotic
solution of the system (64) near the initial singularity,

H =
2

3γ
t−1, x = c11t

− 2
3γ , where c11 ̸= 0,m ∈ R, n < 1,

with (k = 0, γ ̸= 0) or (k = 0,±1, γ =
2

3
), (72)

In fact, solution (72) appears on both first and all-terms dominant decompositions. This
solution for H and x belongs to the phase space of the system given in(64), which is the
set

{
(x,H) ∈ R2 : x ≥ 0, H2 + kx2−2n > 0

}
.

Taking further into account the constraint (63), which represents the first Friedman
equation (1), we can have the asymptotic solutions for H and ρ. In particular, the solution
found in (72) imply,

H =
2

3γ
t−1, ρ =

4

3γ2
cm11t

−2− 2m
3γ + 3kc2+m−2n

11 t
4−4n
3γ

− 2m
3γ , c11 ̸= 0,m ∈ R, n < 1,

with (k = 0, γ ̸= 0) or (k = 0,±1, γ =
2

3
), (73)

It is interesting to mention that the analysis we presented for the varying G case can
also be used for the special case (n = 0,m ∈ R). Taking into account Einstein’s field
equation,

Rij −
1

2
gijR =

8πG

c4
Tij (74)

and (58), the choice to keep G constant as c varies or to keep c constant as G varies can be
considered equivalent. In both cases the ratio 8πG

c4
takes the form aν = am−4n. Of course,



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 10 of 12

the two Friedman equations (1) and (2) lead to a different physical meaning, because c is
multiplying the scale factor term a−2 while G is multiplying the density term ρ. Returning
back to our solutions, if n = 0 and m ̸= 0, the right-hand sides of the Friedman equations,
that we started with, do not depend only on the variation of a and ρ, but they depend
on the variation of G as well, and it is clear that the solution (72) will remain invariant
under n = 0 (n can be any number less than 1), while (73) will give,

H =
2

3γ
t−1, ρ =

4

3γ2
cm11t

−2− 2m
3γ + 3kc2+m

11 t
4
3γ

− 2m
3γ , c11 ̸= 0,m ̸= 0, n = 0,

with (k = 0, γ ̸= 0) or (k = 0,±1, γ =
2

3
). (75)

6. Asymptotic solutions to the flatness problem

Starting with the generalized conservation law (59) with c = c0a
n and G = G0a

m we
obtain,

(am+3γρ)· = 6knc20a
2n+3γ−3ȧ. (76)

Integrating (76), we get,

ρ =
6kna2n−2−m

2n+ 3γ − 2
+Ma−3γ−m, if 2n ̸= 2− 3γ, (77)

or
ρ = 3k(2− 3γ)a−3γ−m ln a+Ma−3γ−m, if 2n = 2− 3γ. (78)

In particular, taking into account the general asymptotic solutions (57) and (73), we are
interested in choosing γ = 2

3 , because this is the only choice that allows k to be 0 or ±1.
Then, (77) and (78) read,

ρ = 3kna2n−2−m +Ma−2−m, if n ̸= 0, (79)

and
ρ = Ma−2−m, if n = 0, (80)

respectively. Next, returning to the first Friedman equation (1) we have,

H2 =
M

3
a−2 + (n− 1)ka2(n−1), if n ̸= 0, (81)

and

H2 = (
M

3
− k)a−2, if n = 0, (82)

Considering (81) at large a, it is clear that the first term dominates the constant curvature
term whenever n ≤ 0, while in (82) k does not play any important role as it can be neglected
by specific values of the integrating constant M . Therefore, in both cases, the flatness
problem can be solved if n ≤ 0. Notice that the value of m does not have any contribution



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 11 of 12

to the results. This is because G only multiplies ρ in (1) and their product eliminates the
constant m.

Our expressions (57) and (73) for constant and varying G respectively, accept solutions
for which n ≤ 0, so they can also be considered as asymptotic solutions to the flatness
problem. In particular, (57) remains the same, but (73) and its special case (75) become,

H = t−1, ρ =
1

3
cm11t

−2−m + 3kc2+m−2n
11 t2−2n−2m, c11 ̸= 0,m ∈ R, n ≤ 0, (83)

and

H = t−1, ρ =
1

3
cm11t

−2−m + 3kc2+m
11 t2−m, c11 ̸= 0,m ̸= 0, n = 0, (84)

respectively.

7. Conclusion

Comparing the results for both constant and varying G by taking only into account
the asymptotic solutions that express a varying speed of light cosmological theory (n ̸= 0),
we find that, in both cases, the associated general solutions (57) and (73) are not invariant
with respect to the time-symmetry t → −t, H → −H, ρ → ρ. This is due to the form of
the energy density ρ in the solutions. In particular, in (57) the expression for ρ contains
odd powers of t. Also, in (73) there are infinitely many combinations of the parameters
m, n, and γ that do not respect the symmetry ρ → ρ (this is also true for (75) when n = 0
and m ̸= 0). This means that for early times the existence of singularity related to the
dynamical systems (13) and (64) creates a time-asymmetry.

In this paper, we have considered the behavior of the VSL cosmological models near the
initial singularity for both G =constant and varying G = G(t). Taking into consideration
various asymptotic conditions that each decomposed system must satisfy in order to have
admissible solutions near the initial singularity, we have presented the form of the series
solutions of the variables of the systems (13) and (64) that are accepted.

It turns out that in both cases of G, the Hubble parameter behaves as t−1 near the
initial singularity. By applying various asymptotic arguments, we were able to built solu-
tions to the field equations in both systems (13) and (64) in the form of a Fuchsian formal
series expansion which are consistent with all other constraints, and are asymptotically
dominant to leading order. Thus, we conclude that the given solutions correspond to
early time attractors of all such homogeneous and isotropic VSL cosmologies. The same
solutions for n ≤ 0 can also be used for the solution to the flatness problem.

References

[1] J Magueijo. New varying speed of light theories. Reports on Progress in Physics,
66:2025–2068, 2003.

[2] J Miritzis and S Cotsakis. Singularities of varying light speed cosmologies. Journal
of Physics: Conference Series, 68:012019, 2007.



D. Trachilis / Eur. J. Pure Appl. Math, 18 (3) (2025), 6581 12 of 12

[3] S Cotsakis and J D Barrow. The dominant balance at cosmological singularities.
Journal of Physics: Conference Series, 68:012004, 2007.

[4] A Albrecht and J Magueijo. Time varying speed of light as a solution to cosmological
puzzles. Physics Letters D, 59:043516, 1999.

[5] J D Barrow and J Magueijo. Varying-α theories and solutions to the cosmological
problems. Physics Letters B, 443:104–110, 1998.

[6] J D Barrow. Cosmologies with varying light speed. Physical Review D, 59:043515,
1999.

[7] J D Barrow. Unusual features of varying speed of light cosmologies. Physics Letters
B, 564:1–7, 2003.

[8] A Goriely. Integrability and Nonintegrability of Dynamical Systems. World Scientific,
University of Arizona, USA, 2001.

[9] S Cotsakis S Kadry G Kolionis and A Tsokaros. Asymptotic vacua with higher
derivatives. Physics Letters B, 755:387–392, 2016.

[10] S Cotsakis G Kolionis and A Tsokaros. The initial state of generalized radiation
universes. Physics Letters B, 721:1–6, 2013.