#### Starobinsky cosmological model in Palatini formalism

Eur. Phys. J. C
Starobinsky cosmological model in Palatini formalism
Aleksander Stachowski 2
Marek Szydłowski 1 2
Andrzej Borowiec 0
0 Institute for Theoretical Physics, Wroclaw University , pl. Maxa Borna 9, 50-204 Wrocław , Poland
1 Mark Kac Complex Systems Research Centre, Jagiellonian University , Ł ojasiewicza 11, 30-348 Kraków , Poland
2 Astronomical Observatory, Jagiellonian University , Orla 171, 30-244 Kraków , Poland
We classify singularities in FRW cosmologies, which dynamics can be reduced to the dynamical system of the Newtonian type. This classification is performed in terms of the geometry of a potential function if it has poles. At the sewn singularity, which is of a finite scale factor type, the singularity in the past meets the singularity in the future. We show that such singularities appear in the Starobinsky model in f ( Rˆ ) = Rˆ + γ Rˆ 2 in the Palatini formalism, when dynamics is determined by the corresponding piecewise-smooth dynamical system. As an effect we obtain a degenerate singularity. Analytical calculations are given for the cosmological model with matter and the cosmological constant. The dynamics of model is also studied using dynamical system methods. From the phase portraits we find generic evolutionary scenarios of the evolution of the universe. For this model, the best fit value of γ = 3γ H02 is equal 9.70 × 10−11. We consider a model in both Jordan and Einstein frames. We show that after transition to the Einstein frame we obtain both the form of the potential of the scalar field and the decaying Lambda term.
1 Introduction
The main aim of the paper is the construction of the
Starobinsky model with a squared term Rˆ 2 in the Palatini
formalism and the investigation of cosmological implications of
this model. In this model the inflation phase of evolution
of the universe can be obtained by the modification of
general relativity in the framework of f ( Rˆ ) modified gravity
theories [1]. In this context, historically the first theory of
inflation was proposed by Starobinsky [2]. In the original
Starobinsky model the term R2/6M 2 was motivated by the
conformal anomaly in the quantum gravity. The problem of
inflation in an f ( R) cosmological model is strictly related
with the choice of frames. The authors of [1] show that
CMB spectra in both Einstein and Jordan frames are
different functions of the number of e-foldings until the end of
inflation.
Inflation is a hypothesis about the existence of a short
but very fast (of exponential type) accelerated growth of the
scale factor a(t ) during the early evolution of the universe,
after the Big-Bang but before the radiation-dominated epoch
[3, 4]. It implies a¨ (t ) > 0. Irregularities in the early epoch
may lead to the formation of structures in the universe due
to the appearance of inflation.
Starobinsky [2] was the first who proposed a very simple
theoretical model with one parameter M (energy scale M )
of such inflation and which is in good agreement with
astronomical data and CMB observation. The Starobinsky model
is representing the simplest version of f ( R) gravity theories
which have been developed considerably in the last decade
[1, 5, 6], whose extra term in the Lagrangian is quadratic in
the scalar curvature. This model predicts the value of spectral
index ns = 0.9603 ± 0.0073, at the 68% CL, with
deviation from scale-invariance of the primordial power spectrum
[7, 8].
The Starobinsky model is also compatible with Planck
2015 data [9] and nicely predicts the number N = 50–60
e-folds between the start and the end of inflation [10].
It has been recently investigated some generalization of
the Starobinsky inflationary model with a polynomial form
of f ( R) = R + 6RM22 + 2λnn (3MR2n)n−1 . It was demonstrated
that the slow-roll inflation can be achieved as long as the
dimensionless coupling λn is sufficiently small [11].
The Starobinsky model becomes generic because the
smallness of the dimensionless coupling constant λn does
not imply that fine-tuning is necessary [11]. The Starobinsky
model was developed in many papers [8, 12–17].
In this paper we develop the idea of endogenous
inflation as an effect of modification of the FRW equation after
the formulation of f (R) cosmological model in the Palatini
formalism.
We are looking for an inflation mechanism as a pure
dynamical mechanism driven by the presence of the
additional term (square of the Ricci scalar) in the Lagrangian,
without necessity of the choice of a frame (Einstein vs.
Jordan frame) [16–18].
In modern cosmology, a most popular trend is to explain
the dark energy and the dark matter in terms of some
substances, of which the nature is unknown up to now.
Einstein was representing the opposite relational point of view
on the description of gravity, in which all substantial forms
should be eliminated. Such a point of view is called
antisubstantialism. Extended f (R) gravity models [6,19] offer
intrinsic or geometric models of both dark matter and dark
energy—the key elements of Standard Cosmological Model.
Therefore, the Einstein idea of relational gravity, in which
dark matter and dark energy can be interpreted as geometric
objects, is naturally realized in f (Rˆ ) extended gravity. The
methods of dynamical system in the context of investigation
dynamics of f (R) models are used since Carroll [19,20].
Unfortunately, the metric formulation of extended
gravity gives rise to fourth order field equations. To avoid this
difficulty, the Palatini formalism can be apply where both
the metric g and the symmetric connection are assumed
to be independent dynamical variables. In consequence, one
gets a system of second order partial differential equations.
The Palatini approach reveals that the early universe inherits
properties of the global CDM evolution.
The Palatini approach has become of some interest lately.
An excellent review of the Palatini f (R) theories can be
found in Olmo’s paper [21]. He has published many other
papers on this topic, namely, about the scalar–tensor
representation of the Palatini theories [22,23]. The other important
papers were on the existence of non-singular solutions in the
Palatini gravity [24,25]. Some more recent papers
concentrate on studying black holes and their singularities in the
Palatini approach [26–30]. Other work which is important
to mention is Flanagan’s papers on the choice of a
conformal frame [31,32]. Pannia et al. considered the impact of the
Starobinsky model in compact stars [33].
In the Palatini gravity action for f (Rˆ ) gravity is introduced
to be
1
S = Sg + Sm = 2
√−g f (Rˆ )d4x + Sm,
(
1
)
where Rˆ = gμν Rˆμν ( ) is the generalized Ricci scalar and
Rˆμν ( ) is the Ricci tensor of a torsionless connection . In
this paper, we assume that 8π G = c = 1. The equation
of motion obtained from the first order Palatini formalism
reduces to
(
2
)
(
3
)
(
4
)
(
5
)
(
6
)
(
7
)
(
8
)
(
9
)
2 δLm is matter energy-momentum tensor,
where Tμν = − √−g δgμν
i.e. one assumes that the matter minimally couples to the
metric. As a consequence the energy-momentum tensor is
conserved, i.e.: ∇μTμν = 0 [34]. In Eq. (
3
) ∇ˆ α means the
covariant derivative calculated with respect to . In order to
solve Eq. (
3
) it is convenient to introduce a new metric,
√
hhμν =
√−g f (Rˆ )gμν
for which the connection = LC(h) is a Levi-Civita
connection. As a consequence in dim M = 4 one gets
hμν = f (Rˆ )gμν ,
i.e. both metrics are related by the conformal factor. For this
reason one should assume that the conformal factor f (Rˆ ) =
0, so it has strictly positive or negative values.
Taking the trace of (
2
), we obtain additional so called
structural equation
f (Rˆ )Rˆ − 2 f (Rˆ ) = T .
where T = gμν Tμν . Because of cosmological applications
we assume that the metric g is FRW metric
ds2 = −dt 2 +a2(t )
1
1 − kr 2 dr 2 + r 2(dθ 2 + sin2 θ dφ2) ,
1
f (Rˆ )Rˆμν − 2 f (Rˆ )gμν = Tμν ,
∇ˆ α(
√−g f (Rˆ )gμν ) = 0,
where a(t ) is the scale factor, k is a constant of spatial
curvature (k = 0, ±1), t is the cosmological time. For simplicity
of presentation we consider the flat model (k = 0).
As a source of gravity we assume a perfect fluid, with the
energy-momentum tensor
T μ
ν = diag(−ρ , p, p, p),
where p = wρ, w = const is a form of the equation of
state (w = 0 for dust and w = 1/3 for radiation). Formally,
effects of the spatial curvature can also be included into the
k a−2, with
model by introducing a curvature fluid ρk = − 2
the barotropic factor w = − 31 ( pk = − 3 ρk). From the
con1
μ
servation condition Tν;μ = 0 we obtain ρ = ρ0a−3(1+w).
Therefore the trace T reads
T =
ρi,0(3wi − 1)a(t )−3(1+wi ).
i
f (Rˆ ) = Rˆ + γ Rˆ 2,
induced by the first three terms in the power series
decomposition of an arbitrary function f (R). In fact, since the terms
Rˆ n have different physical dimensions, i.e. [Rˆ n] = [Rˆ m ] for
n = m, one should take instead the function Rˆ0 f (Rˆ /Rˆ0)
for constructing our Lagrangian, where Rˆ0 is a constant and
[Rˆ0] = [Rˆ ]. In this case the power series expansion reads
Rˆ0 f (Rˆ /Rˆ0) = Rˆ0 n=0 αn(Rˆ /Rˆ0)n = n=0 α˜ n Rˆ n, where
the coefficients αn are dimensionless, while [α˜ n] = [Rˆ ]1−n
are dimension full.
From the other hand the Lagrangian (
10
) can be viewed as
a simplest deviation, by the quadratic Starobinsky term, from
the Lagrangian Rˆ which provides the standard cosmological
model a.k.a. CDM model. A corresponding solution of the
structural equation (
6
)
Rˆ = −T ≡ 4ρ ,0 + ρm,0a−3.
is, in fact, exactly the same as for the CDM model, i.e.
with γ = 0. However, the Friedmann equation of the CDM
model (with dust matter, dark energy and radiation)
b2
b
×
+
d 2
b + 2
(K − 3)(K + 1)
2b
r,0a−4
+
k ,
(
10
)
(
11
)
(
12
)
(
13
)
(
14
)
(
15
)
(
16
)
ρr,0a−4 + ρm,0a−3 + ρ ,0
is now hardly affected by the presence of quadratic term.
More exactly a counterpart of the above formula in the model
under consideration looks as follows:
γ ( m,0a−3 +
,0)2
+ ( m,0a−3 +
,0)
In what follows we consider visible and dark matter ρm in
the form of dust w = 0, dark energy ρ with w = −1 and
radiation ρr with w = 1/3.
Because a form of the function f (Rˆ ) is unknown, one
needs to probe it via ensuing cosmological models. Here
we choose the simplest modification of the general relativity
Lagrangian,
ρ ,0 ,
,0 = 3H02
,
m,0 = 3ρHm,002 ,
(
17
)
(
18
)
(
19
)
(
20
)
(
21
)
b = f (Rˆ ) = 1 + 2 γ ( m,0a−3 + 4
1 db
d = H dt = −2 γ ( m,0a−3 +
,0),
,0)(3 − K )
From the above one can check that the standard model (
12
)
can be reconstructed in the limit γ → 0. The study of this
generalized Friedmann equation is a main subject of our
research.
The paper is organized as follows. In Sect. 2, we consider
the Palatini approach in the Jordan and Einstein frame. In
Sect. 3, we present some generalities concerning
dynamical systems of Newtonian type, and their relations with the
Palatini–Starobinsky model. Section 4, is devoted to the
classification of cosmological singularities with special attention
on Newtonian type systems represented by potential function
V (a). We adopt the Fernandes-Jambrina and Lazkoz
classification of singularities [35] to these systems using the notion
of elasticity of the potential function with respect the scale
factor. In Sect. 5, we will analyze the singularities in the
Starobinsky model in the Palatini formalism. This system
requires the form of piecewise-smooth dynamical system.
Statistical analysis of the model is presented in Sect. 6. In
Sect. 7, we shall summarize obtained results and draw some
conclusions.
2 The Palatini approach in different frames (Jordan vs.
Einstein frame)
Because the effect of acceleration can depend on a choice of a
frame [36] this section is devoted to showing the existence of
the inflation effect if the model is considered in the Einstein
frame.
The action (
1
) is dynamically equivalent to the first order
Palatini gravitational action, provided that f (Rˆ ) = 0 [1,6,
17]
1
S(gμν , ρλσ , χ ) = 2
d4x √−g f (χ )(Rˆ − χ ) + f (χ )
+ Sm (gμν , ψ ),
Introducing a scalar field = f (χ ) and taking into account
the constraint χ = Rˆ , one gets the action (
22
) in the following
form:
1
S(gμν , ρλσ , ) = 2
d4x √−g
Rˆ − U ( )
+ Sm (gμν , ψ ),
(
22
)
(
23
)
where the potential U ( ) is defined by
U f ( ) ≡ U ( ) = χ ( )
− f (χ ( ))
with = d df(χχ) and Rˆ ≡ χ = dUd( ) .
The Palatini variation of the action (
23
) gives rise to the
following equations of motion:
1
Rˆμν − 2 gμν Rˆ
1
+ 2 gμνU ( ) − Tμν = 0,
∇ˆ λ(
√−g gμν ) = 0,
Rˆ − U ( ) = 0.
Equation (25b) implies that the connection ˆ is a metric
connection for a new metric g¯μν = gμν ; thus Rˆμν = R¯μν , R¯ =
g¯μν R¯μν = −1 Rˆ and g¯μν R¯ = gμν Rˆ . The g-trace of (25a)
produces a new structural equation
2U ( ) − U ( )
= T .
Now Eqs. (25a) and (25c) take the following form:
1 1
R¯μν − 2 g¯μν R¯ = T¯μν − 2 g¯μνU¯ ( ),
R¯ − ( 2 U¯ ( )) = 0,
where we introduce U¯ (φ) = U (φ)/ 2, T¯μν =
the structural equation can be replaced by
−1Tμν and
U¯ ( ) + T¯ = 0 .
In this case, the action for the metric g¯μν and scalar field
is given by
1
S(g¯μν , ) = 2
d4x
−g¯ R¯ − U¯ ( ) + Sm ( −1g¯μν , ψ ),
where we have to take into account a non-minimal coupling
between and g¯μν
T¯ μν = − √−g¯ δg¯μν
2
δ
Sm = (ρ¯+ p¯)u¯μu¯ν + p¯ g¯μν =
−3T μν ,
(
31
)
u¯μ = − 21 uμ, ρ¯ = −2ρ , p¯ = −2 p, T¯μν =
−1Tμν , T¯ = −2T (see e.g. [17,37]).
In FRW case, one gets the metric g¯μν in the following
form:
ds¯2 = −dt¯2 + a¯ 2(t ) dr 2 + r 2(dθ 2 + sin2 θ dφ2) ,
(
32
)
1
where1 dt¯ = (t ) 2 d t and new scale factor a¯ (t¯) =
(t¯) 2 a(t¯). Ensuing cosmological equations (in the case of
the barotropic matter) are given by
a
¨
3H¯ 2 = ρ¯ + ρ¯m , 6 a¯ = 2ρ¯ − ρ¯m (1 + 3w)
¯
where
ρ¯ = 21 U¯ ( ), ρ¯m = ρ0a¯ −3(1+w) 21 (3w−1)
and w = p¯m/ρ¯m = pm/ρm. In this case, the conservation
equations has the following form:
ρ˙¯m + 3H¯ ρ¯m(1 + w) = −ρ¯˙ .
Let us consider the Starobinsky–Palatini model in the above
framework. The potential U¯ is described by the following
formula:
U¯ ( ) = 2ρ¯ ( ) =
1
4γ + 2λ
1 1 1
2 − 2γ
1
+ 4γ . (
36
)
Figure 1 presents the relation ρ¯ ( ). Note that the
function ρ¯ has the same shape like the Starobinsky potential.
The function ρ¯ ( ) has the minimum for
min = 1 + 8γ λ.
In general, the scalar field
(a¯ ) is given by (cf. (
11
))
= 1 + 2γ Rˆ = 1 + 8γ λ + 2γρm − 6γ pm .
Because ρ¯m = −2ρm , p¯m =
account (
34
) one gets
−2 pm , and taking into
2γ (1 − 3w)ρ0a¯ −3(1+w) 23 (w+1) −
+ 1 + 8γ λ = 0. (
39
)
the algebraic equation determining the function (a¯ ) for a
given barotropic factor w. This provides an implicit
dependence (a¯ ). In order to get it more explicit one needs to solve
(
39
) for some interesting values w. For example in the case
of dust we obtain the third order polynomial equation
1
2γ + 4
y3 − 21γ y + ρ0wa¯ −3 = 0
1
where y = − 2 .
The evolution of (a¯ ), at the beginning of the inflation
epoch, is presented in Fig. 2.
For γ ≈ 0, the potential U¯ can be approximated as U¯ =
−ρ¯m + 41γ . In this case the Friedmann equation can be written
as
(
33
)
(
34
)
(
35
)
(
37
)
(
38
)
(
40
)
3H¯ 2 = ρ¯2m + 81γ .
1.4 1068
1.2 1068
1.0 1068
8.0 1067
6.0 1067
4.0 1067
2.0 1067
5
10
15
20
Fig. 2 Illustration of the typical evolution of with respect to ln(a¯ )
at the beginning of the inflation epoch. We assume that γ = 1.16 ×
10−69 s2 and a¯0 = 1 at the beginning of the inflation epoch
ln a
In the case of ρ¯m = 0, ρ¯ is constant and the Friedmann
equation has the following form:
1
3 H¯ 2 = 8γ .
In this model the inflation phenomenon appears when the
value of the parameter γ is close to zero and the matter ρ¯m is
negligible with comparison to ρ¯ . In this case the
approximate number of e-foldings is given by the following formula:
N = Hinit(t¯fin − t¯init) =
t¯fin − t¯init .
√24γ
The number of e-folds N should be equal 50 ∼ 60 in the
inflation epoch [10]. In this model we obtain N = 60, when
γ = 1.16×10−69 s2 and the timescale of the inflation is equal
10−32 s [38]. The relation between γ and the approximate
number of e-foldings N is presented in Fig. 3.
(
41
)
(
42
)
ρ m ln a
4 10105
3 10105
2 10105
ln a
ρ ln a
3.090 10107
3.085 10107
3.080 10107
1.0
0.5
0.5
1.0
1.5
2.0
ln a
Fig. 5 Illustration of the typical evolution of ρ¯φ with respect to ln(a¯ )
at the beginning of the inflation epoch. We assume that γ = 1.16 ×
10−69 s2 and a¯0 = 1 at the beginning of the inflation epoch. ρ¯ is
expressed in units of s2kMmp2c2 . Note that during the inflation ρ¯φ ≈ const
a t
7
6
5
4
3
2
1
Fig. 6 Illustration of the typical evolution of a¯ with respect to t¯ at
the beginning of the inflation epoch. We assume that γ = 1.16 ×
10−69 s2 and a¯0 = 1 at the beginning of the inflation epoch. The time
t¯ is expressed in seconds
p
¯
= w(a)ρ¯ ,
where the function w(a) is given by the expression
w(a) = −1 − √3√
ρ
˙
¯
ρ¯m + ρ¯ ρ
1 d ln ρ
= −1 − 3 H¯ dt¯
The diagram of the coefficient of equation of state w(a),
at the beginning the inflation epoch, is presented in Fig. 8.
Note that the function w(a), for the late time, is a constant
and equal −1.
The action (
23
) can be rewritten in the Jordan frame
(gμν , ) as
4 √
d x
−g
R + 2
3
∂μ
∂ μ
− U ( ) ,
t
(
44
)
H ln a
Fig. 7 Illustration of the typical evolution of H¯ with respect to ln(a¯ )
at the beginning of the inflation epoch. We assume that γ = 1.16 ×
10−69 s2 and a¯0 = 1 at the beginning of the inflation epoch. H¯ is
expressed in units of s kMmpc . Note that, for the late time, H¯ can be treated
as a constant
ln a
ln a
Fig. 8 Illustration of the typical evolution of wφ with respect to ln(a¯ ).
We assume that γ = 1.16 × 10−69 s2 and a¯0 = 1 at the beginning of
the inflation epoch. Note that during the inflation wφ ≈ −1
where R is the metric Ricci scalar, = f ( Rˆ ), Rˆ = χ ( ).
We obtain the Brans–Dicke action with the coupling
parameter ω = − 23 in the Jordan frame. The equations of
motion take the form
1
Rμν − 2 gμν R
3
− 4
gμν ∇σ
∇
σ
3
+ 2
∇μ
∇ν
+gμν
3
R −
− ∇μ∇ν
1
+ 2 gμν U (φ) − κ Tμν = 0 , (47a)
3 Singularities in cosmological dynamical systems of
Newtonian type
which is called the acceleration equation. It is easily to check
that
There is a class of cosmological models, of which the
dynamics can be reduced to a dynamical system of the Newtonian
type. Let consider a homogeneous and isotropic universe with
a spatially flat space-time metric of the form
ds2 = dt 2 − a2(t ) dr 2 + r 2(dθ 2 + sin2 θ dφ2) ,
(
48
)
where a(t ) is the scale factor and t is the cosmological time.
Let us consider the energy-momentum tensor Tνμ for the
perfect fluid with energy density ρ(t ) and pressure p(t ) as a
source of gravity. In this case the Einstein equations assume
the form of the Friedmann equations,
(
56
)
(
57
)
(
58
)
(
59
)
(
60
)
(
61
)
(
62
)
where V (a) is given by (
53
) provided that the conservation
equation (
51
) is fulfilled.
Due to Eq. (
56
) the evolution of the universe can be
interpreted as the motion of a fictitious particle of unit mass in the
potential V (a). Here a(t ) plays the role of a position
variable. The equation of motion (
56
) assumes a form analogous
to the Newtonian equation of motion.
If we know the form of the effective energy density then
we can construct the form of the potential V (a), which
determines the whole dynamics in the phase space (a, a˙ ). In this
space the Friedmann equation (
52
) plays the role of a first
integral and determines the phase space curves representing
the evolutionary paths of the cosmological models. The
diagram of the potential V (a) contains all information needed
to construction a phase space portrait. In this case the phase
space is two-dimensional,
(a, a˙ ) : a˙22 + V (a) = − k2 .
In the general case of an arbitrary potential, the dynamical
system which describes the evolution of a universe takes the
form
a˙ = x ,
x˙ = −
∂ V (a)
∂a
We shall study the system above using the theory of
piecewise-smooth dynamical systems. Therefore it is
assumed that the potential function, except some isolated
(singular) points, belongs to the class C 2(R+).
The lines x22 + V (a) = − k2 represent possible evolutions
of the universe for different initial conditions. Equations (
58
)
and (
59
) can be rewritten in terms of dimensionless variables
if we replace the effective energy density ρeff by the density
parameter:
where dot denotes differentiation with respect to the cosmic
time t , H ≡ aa˙ is the Hubble function. In our notation we use
the natural system of units in which 8π G = c = 1.
We assume ρ(t ) = ρ(a(t )) as well as p(t ) = p(a(t )),
i.e. both energy density and pressure depend on the cosmic
time through the scale factor a(t ). The conservation condition
T;μμν = 0 reduces to
It would be convenient to rewrite (
49
) in the equivalent form
In (
53
) ρ(a) plays the model role of an effective energy
density. For example for the standard cosmological model (
12
)
ρm,0a−3 + ρ ,0 ,
and ρm = ρm,0a−3. Equation (
50
) is
3a2
ρ = 3H 2 = a˙2 ,
2a¨ a˙ 2 ,
p = − a − a2
ρ˙ = −3H (ρ + p).
a˙ 2 = −2V (a),
where
V (a) = −
ρ(a)a2
6
Any cosmological model can be identified by its form of
the potential function V (a) depending on the scale factor a.
From the Newtonian form of the dynamical system (
58
)–(
59
)
one can see that all critical points correspond to vanishing of
r.h.s. of the dynamical system x0 = 0, ∂V (a)
∂a |a=a0 .
Therefore all critical points are localized on the x -axis, i.e. they
represent a static universe.
Because of the Newtonian form of the dynamical system
the character of critical points is determined from the
characteristic equation of the form
2
a + det A|x0=0, ∂V∂a(a) |a0 =0 = 0,
where det A is the determinant of the linearization matrix
calculated at the critical points, i.e.
det A =
∂2V (a)
∂a2 |a0, ∂V∂a(a) |a0 =0.
From Eqs. (
64
) and (
65
) one can conclude that only
admissible critical points are of saddle type if ∂2V (a)
∂a2 |a=a0 < 0 or
of center type if ∂2∂Va(2a) |a=a0 > 0.
If the shape of the potential function is known (from
knowledge of the effective energy density), then it is
possible to calculate the cosmological functions in exact form,
t =
a
da
√−2V (a)
,
H (a) = ±
−
2V (a)
a2
,
aa¨ 1 d ln(−V )
q = − a˙ 2 = 2 d ln a
the deceleration parameter, the effective barotropic factor
where t → τ = |H0|t and
V˜ (a) = −
.
peff 1
weff(a(t )) = ρeff = − 3
d ln(−V )
d ln a
+ 1 ,
the parameter of deviation from de Sitter universe [35]
1 d ln(−V )
h(t ) ≡ −(q(t ) + 1) = 2 d ln a
− 1
(note that if V (a) = − 6a2 , h(t ) = 0), the effective matter
density and pressure
ρeff = −
,
peff =
2V (a)
a2
d ln(−V )
d ln a
+ 1 ,
and, finally, the Ricci scalar curvature for the FRW metric
(
48
),
R =
6V (a)
a2
d ln(−V )
d ln a
+ 2 .
From the formulas above one can observe that the most of
them depend on the quantity
(
72
)
(
73
)
(
74
)
(
75
)
(
63
)
(
64
)
(
65
)
(
66
)
(
67
)
(
68
)
(
69
)
(
70
)
(
71
)
Iν (a) =
d ln(−V )
d ln a
This quantity measures the elasticity of the potential
function, i.e. indicates how the potential V (a) changes if the scale
factor a changes. For example, for the de Sitter universe
−V (a) ∝ a2, the rate of growth of the potential is 2% as
the rate of growth of the scale factor is 1%.
In the classification of the cosmological singularities by
Fernandez-Jambrina and Lazkoz [35] a crucial role is played
by the parameter h(t ). Note that in a cosmological sense this
parameter is
In this approach the classification of singularities is based
on generalized power and asymptotic expansion of the
barotropic index w in the equation of state (or
equivalently of the deceleration parameter q) in terms of the time
coordinate.
4 Degenerated singularities—new type (VI) of singularity—sewn singularities
Recently, due to the discovery of an accelerated phase in the
expansion of our universe, many theoretical possibilities for
future singularities are seriously considered. If we assume
that the universe expands following the Friedmann equation,
then this phase of expansion is driven by dark energy—a
hypothetical fluid, which violates the strong energy
condition. Many of the new types of singularities were classified
by Nojiri et al. [40]. Following their classification the type of
singularity depends on the singular behavior of the
cosmological quantities like the scale factor a, the Hubble parameter
H , the pressure p and the energy density ρ:
– Type 0: ‘Big crunch’. In this type, the scale factor a is
vanishing and there is blow-up of the Hubble parameter
H , energy density ρ and pressure p.
– Type I: ‘Big rip’. In this type, the scale factor a, energy
density ρ and pressure p are blown up.
– Type II: ‘Sudden’. The scale factor a, energy density ρ
and Hubble parameter H are finite and H˙ and the pressure
p are divergent.
– Type III: ‘Big freeze’. The scale factor a is finite and the
Hubble parameter H , energy density ρ and pressure p
are blown up [41] or divergent [42].
– Type IV. The scale factor a, Hubble parameter H , energy
density ρ, pressure p and H˙ are finite but higher
derivatives of the scale factor a diverge.
– Type V. The scale factor a is finite but the energy density
ρ and pressure p vanish.
Following Królak [43], big rip and big crunch singularities
are strong whereas sudden, big freeze and type IV are weak
singularities.
In the model under consideration the potential function
and/or its derivative can diverge at isolated points (value
of the scale factor). Therefore the classification mentioned
before has application only for a single component of
piecewise-smooth potential. In our model the dynamical
system describing the evolution of a universe belongs to the
class of a piecewise-smooth dynamical systems. As a
consequence new types of singularities at finite scale factor as
can appear for which ∂∂Va (as ) does not exist (is not
determined). This implies that the classification of singularities
should be extended to the case of non-isolated
singularities.
Let us illustrate this idea on the example of a freeze
singularity in the Starobinsky model with the Palatini formalism
(previous section). Such a singularity has a complex
character and in analogy to the critical point we called it degenerate.
Formally it is composed of two types III singularities: one
in the future and another one in the past. If we consider the
evolution of the universe before this singularity we detect an
isolated singularity of type III in the future. Conversely if we
consider the evolution after the singularity, then going back
in time we meet a type III singularity in the past. Finally,
at the finite scale factor the two singularities will meet. For
a description of behavior near the singularity one considers
the t = t (a) relation. This relation has a horizontal inflection
point and it is natural to expand this relation in a Taylor series
near this point at which ddat = H1a is zero. For the freeze
singularity, the scale factor remains constant as, ρ and H blow
up and a¨ is undefined. It this case, the degenerate
singularity of type III is called sewn (non-isolated) singularity. We,
therefore, obtain [44]
t − ts
1 d2t 2
± 2 da2 |a=asing (a − asing) .
(
76
)
0.0002
0.0004
t
Fig. 10 Illustration of a sewn sudden singularity. The model with
negative γ has a mirror symmetry with respect to the cosmological time.
Note that the spike on the diagram shows a discontinuity of the function
∂∂Va . Note the existence of a bounce at t = 0
The above formula combines two types of behavior near the
freeze singularities in the future,
a − asing ∝ −(tsing − t )1/2 for t → tsing−
and in the past
a − asing ∝ +(t − tsing)1/2 for t → tsing+ .
(
77
)
(
78
)
Figure 9 illustrates the behavior of the scale factor in
cosmological time in neighborhood of a pole of the potential
function. Diagram of a(t ) is constructed from the dynamics
in two disjoint region {a : a < as} and {a : a > as}. Figure 10
presents the behavior of the scale factor in the cosmological
time in a neighborhood of the sudden singularity.
In the model under consideration another type of sewn
singularity also appears. It is a composite singularity with
two sudden singularities glued at the bounce when a = amin.
In this singularity the potential itself is a continuous
function while its first derivative has a discontinuity. Therefore,
the corresponding dynamical system represents a
piecewisesmooth dynamical system.
The problem of C 0 metric extension beyond the future
Cauchy horizon, when the second derivative of the metric is
inextendible, was discussed in work of Sbierski [45]. In the
context of FLRW cosmological models, Sbierski’s
methodology was considered in [46].
afsing
5 Singularities in the Starobinsky model in the Palatini formalism
In our model, one finds two types of singularities, which
are a consequence of the Palatini formalism: the freeze and
sudden singularity. The freeze singularity appears when the
b
multiplicative expression b+d/2 , in the Friedmann equation
(
13
), is equal to infinity. So we get a condition for the freeze
singularity: 2b+d = 0, which produces a pole in the potential
function. It appears that the sudden singularity appears in our
b
model when the multiplicative expression b+d/2 vanishes.
This condition is equivalent to the case b = 0.
The freeze singularity in our model is a solution of the
algebraic equation
2b + d = 0
⇒ f (K ,
,0,
γ ) = 0
or
where K ∈ [0, 3).
The solution of the above equation is
Kfreeze = 3 + 3 γ ( m+1
,0)
,0
.
1
,0)
,0
+ 1 = 0,
From Eq. (
81
), we can find an expression for a value of the
scale factor for the freeze singularity
afreeze =
8
1 −
,0 +
,0
1
γ ( m+
,0)
.
The relation between afreeze and positive γ is presented in
Fig. 11.
The sudden singularity appears when b = 0. This leaves
us with the following algebraic equation:
1 + 2 γ ( m,0a−3 +
,0)(K + 1) = 0.
(83)
x = −
∂ V (a)
∂ a
,
which, in fact, becomes a (degenerate) critical point and a
bounce at the same time. The relation between asing and
negative γ is presented in Fig. 12.
Let V = − a22 γ c2h (K −32)b(K +1) + ch + k . We can
rewrite dynamical system (
58
)–(
59
) as
asuddsing
where ≡ ddσ = b+b d2 ddτ is a new parametrization of time.
We can treat the dynamical system (86)–(87) as a sewn
dynamical system [47, 48]. In this case, we divide the phase
portrait into two parts: the first part is for a < asing and the
second part is for a > asing. Both parts are glued along the
singularity.
For a < asing, dynamical system (86)–(87) can be
rewritten in the corresponding form,
where V2 = V η(a − as ). The phase portraits, for dynamical
system (86)–(87), are presented in Figs. 13 and 14. Figure
x = −
∂ V1(a)
∂ a
,
x = −
∂ V2(a)
∂ a
,
where V1 = V (−η(a −as ) + 1) and η(a) notes the Heaviside
function.
For a > asing, in an analogous way, we get the following
equations:
20
15
10
5
5
10
15
k
1
(88)
(89)
(90)
(91)
Fig. 14 The phase portrait of the system (86)–(87) for negative γ .
The scale factor a is in logarithmic scale. The red trajectories
represent a spatially flat universe. Trajectories under the top red trajectory
and below the bottom red trajectory represent models with a negative
spatial curvature. Trajectories between the top and bottom red
trajectory represent models with the positive spatial curvature. The dashed
line b = 0 corresponds to the sudden singularity. The shaded region
represents trajectories with b < 0. If we assume that f (R) > 0 then
this region can be removed. The phase portrait possesses the symmetry
a˙ → −a˙ and in consequence this singularity presents a bounce. This
symmetry can be used to identify the corresponding points on the b-line.
The critical point (
1
) represents the static Einstein universe. The phase
portrait belongs to the class of sewn dynamical systems [49]
13 shows the phase portrait for positive γ , while Fig. 14
shows the phase portrait for negative γ .
In Fig. 13 there are two critical points labeled ‘1’ and ‘2’
at the finite domain. They are both saddle points. These
critical points correspond to a maximum of the potential
function. The saddle point ‘2’ possesses the homoclinic closed
orbit starting from it and returning to it. This orbit
represents an emerging universe from the static Einstein
universe and approaching it again. During the evolution this
universe (orbit) goes two times through the freeze
singularity. The region bounded by the homoclinic orbit contains
closed orbits representing the oscillating universes. A
diagram of the evolution of scale factor for closed orbit is
presented by Fig. 15. It is also interesting that trajectories in
the neighborhood of straight vertical line of freeze
singularities undergo short time inflation x = const. The
characteristic number of e-foldings from tinit to tfin of this inflation
period N = Hinit(tfin − tinit) (see Eq. (3.13) in [1]) with
respect to γ is shown in Fig. 16. This figure illustrates
the number of e-foldings is too small to obtain the inflation
effect.
a σ
H (zi )obs − H (zi )th 2
In this paper, we use the likelihood function for observations
of CMB [9] and lensing by Planck, and low- polarization
from the WMAP (WP) in the following form:
ln LCMB+lensing = − 21 (xth − xobs)C−1(xth − xobs),
where C is the covariance matrix with the errors, x is a vector
of the acoustic scale l A, the shift parameter R and bh2 where
BAO observations such as Sloan Digital Sky Survey
Release 7 (SDSS DR7) dataset at z = 0.275 [51], 6dF Galaxy
Redshift Survey measurements at redshift z = 0.1 [52], and
WiggleZ measurements at redshift z = 0.44, 0.60, 0.73 [53]
have the following likelihood function:
1
ln LBAO = − 2
dobs
rs (zd )
− DV (z)
C−1 dobs
rs (zd )
,
− DV (z)
Fig. 15 Illustration of the evolution of a(σ ) for closed orbit which is
contained by the homoclinic orbit, where σ = b+b d2 t is a
reparametrization of time. We choose s×Mpc/(100×km) as a unit of σ
0
Fig. 16 Diagram of the relation between positive γ and the
approximate number of e-foldings N = Hinit(tfin − tinit) from tinit to tfin
6 Observations
In this paper we perform statistical analysis using the
following astronomical observations: observations of 580
supernovae of type Ia, BAO, measurements of H (z) for galaxies,
Alcock–Paczyn´ski test, measurements of CMB and lensing
by Planck and low by WMAP.
The likelihood function for observations of supernovae of
type Ia [50] is given by the following expression:
1
ln LSNIa = − 2 [ A − B2/C + ln(C/(2π ))],
(92)
where A = (μobs − μth)C−1(μobs − μth), B = C−1(μobs −
μth), C = TrC−1 and C is a covariance matrix for
observations of supernovae of type Ia. The distance modulus is
defined by the formula μobs = m − M (where m is the
apparent magnitude and M is the absolute magnitude of
observations of supernovae of type Ia) and μth = 5 log10 DL + 25
(where the luminosity distance is DL = c(1 + z) 0z Hd(zz) ).
where rs (zd ) is the sound horizon at the drag epoch [54,55].
For the Alcock–Paczynski test [56,57] we used the
following expression for the likelihood function:
1
ln L AP = − 2
i
A Pth(zi ) − A Pobs(zi ) 2
.
π
l A = rs (z∗) c
0
R =
m,0 H02
z∗ dz
H (z )
z∗ dz
0
H (z )
,
where z∗ is the redshift of the epoch of the recombination
[54].
The total likelihood function is expressed in the following
form:
Ltot = LSNIa LBAO LAP L H(z) LCMB+lensing.
To estimate model parameters, we use our own code
CosmoDarkBox. The Metropolis–Hastings algorithm [70,71] is
used in this code.
Table 1 The best fit and errors for the estimated model for the
positive γ with m,0 from the interval (0.27, 0.33), γ from the interval
(0.0, 2.6 × 10−9) and H0 from the interval [66.0 (km/(s Mpc)), 70.0
(km/(s Mpc))]. b,0 is assumed as 0.048468. The redshift of matter–
radiation equality is assumed as 3395. H0, in the table, is expressed in
km/(s Mpc). The value of reduced χ 2 of the best fit of our model is
equal 0.187066 (for the CDM model 0.186814)
95% CL
Fig. 17 The intersection of the likelihood functions of two model
parameters ( γ , m,0) with the marked 68 and 95% confidence levels
Fig. 18 The intersection of the likelihood functions of two model
parameters ( γ , H0) with the marked 68 and 95% confidence levels
Table 1 shows the values of parameters for the best fit with
errors. Figures 17 and 18 show the intersection of a likelihood
function with the 68 and 95% confidence level projections
on the ( γ , m,0) and ( γ , H0) planes.
In this paper, we use the Bayesian information criterion
(BIC) [72, 73], for comparison of our model with the CDM
model. The expression for BIC is defined as
BIC = χ 2 + j ln n,
where χ 2 is the value of χ 2 in the best fit, j is the number of
model parameters (our model has three parameters, CDM
model has two parameters) and n is the number of data points
(n = 625) which are used in the estimation.
For our model, the value of BIC is equal 135.668 and for
the CDM model BIC CDM = 129.261. So
BIC=BICBIC CDM is equal 6.407. The evidence for the model is
(100)
strong [73] if BIC is higher than 6. So, in comparison to our
model, the evidence in favor of the CDM model is strong,
but we cannot absolutely reject our model.
7 Conclusions
In this paper, we demonstrated that evolution of the
Starobinsky model with a quadratic term R2 gives rise to the
description of dynamics in terms of piecewise-smooth dynamical
systems, i.e., systems whose the phase space is partitioned
into different regions, each of them associated to a different
smooth functional form of the system of a Newtonian type.
Different regions of the phase space correspond to different
forms of the potential separated by singularities of the type
of poles.
Our idea was to obtain inflation as an endogenous effect
of the dynamics in the Palatini formalism. While the effect of
inflation appears in the model under consideration a sufficient
number of e-folds are not achieved and the additional effect
of amplification is required. Note that this type of inflation is
a realization of the idea of singular inflation [74–77]. In our
model inflation is driven by the freeze degenerate singularity
(the extension of a type III isolated singularity).
We show that the dynamics of the model can be analyzed
in terms of two-dimensional dynamical systems of the
Newtonian type. In this approach, in the diagram of the
potential of a fictitious particle, the evolution of the universe
contains all information which is needed for an investigation of
singularities in the model. Note that they are not isolated
singularities which were classified into five types but rather
double singularities glued in one point of the evolution at
a = asing. Appearance of such types of singularities is
typical for piecewise-smooth dynamics describing the model
evolution. We call this type sewn singularities in analogy to
sewn dynamical systems [78, 79].
We investigated the model with f ( Rˆ ) = Rˆ + γ Rˆ 2, where
γ assumes the positive or negative values. While the
dynamics of this class of models depend crucially on the sign of the
parameter γ in the early universe for the late time we obtain
the behavior consistent with the CDM model.
Note that in the model with positive γ , the phase space is a
sum of two disjoint domains which boundary represents the
double freeze singularity (cf. Fig. 13). In the first domain the
evolution starts from a big bang followed by the deceleration
phase; then it changes to acceleration (early acceleration ≡
inflation) after reaching a maximum of the potential function.
In the second domain, on the right from the vertical line of the
freeze singularity, the universe decelerates and after reaching
another maximum starts to accelerate again. This last eternal
acceleration corresponds to the present day epoch called the
dark energy domination epoch. Two phases of deceleration
and two phases of acceleration are key ingredients of our
model. While the first phase models a transition from the
matter domination epoch to inflation the second phase models
a transition from the second matter dominated epoch toward
the present day acceleration.
As De Felice and Tsujikawa have noted [1, p. 24] the
applications of f ( R) theories should be focused on construct of
viable cosmological models, for which a sequence of
radiation, matter and accelerating epochs is realized. All these
epochs are also presented in the model under consideration
but, for negative γ (negative squared M 2 for the scalar field),
some difficulties appear in the interpretation of the phase
space domain {a : a < asing}. The size of this domain will
depend on the value of the parameter γ and this domain
vanishes as we are going toward γ equal zero.
On the other hand it is well known that violation of
condition f > 0 gives rise to the negative values of M 2. We do
Rˆ Rˆ
not assume this condition but we require that f R > 0 to avoid
ˆ
the appearance of ghosts (see Sect. 7.4 in [1]). In our case,
statistical analysis favors a model with f R > 0 ( γ > 0)
ˆ
rather than a model with f R < 0 ( γ < 0). In other words,
ˆ
statistical analysis favors the case without ghosts.
In order to obtain deeper insight into the model we have
also performed complementary investigations in the Einstein
frame. In this case we find that the model is reduced to the
FRW cosmological model with the selfinteracting scalar field
and the vanishing part of the kinetic energy. Therefore from
the Palatini formulation we obtain directly the form of the
potential and the (implicit) functional dependence between
the scalar field and the scale factor. Moreover, we obtain the
parametrization of the decaying cosmological constant.
Due to a time-dependent cosmological constant the model
evolution can be described in terms of an interaction between
the matter and the decaying lambda terms. We study how
the energy is transferred between the sectors and how the
standard scaling relation for matter is modified.
We pointed out that the consideration of the Starobinsky
model in the Einstein frame gives rise to new interesting
properties from the cosmological point of view; similar to
the original (metric) the Starobinsky model is very
important for the explanation of inflation. The model under the
consideration gives rise analogously to the running
cosmological term. This fact seems to be interesting in the context
of an explanation of the cosmological constant problem.
Detailed conclusions coming from our analysis are the
following:
– We show that the interaction between two sectors: the
matter and the decaying vacuum, appears naturally in the
Einstein frame. For the model formulated in the Jordan
frame this interaction is absent.
– Inflation appears in our model formulated in the Einstein
frame, when the parameter γ is close to zero and the
density of matter is negligible in comparison to ρ¯ .
– In our model in the Einstein frame, the potential U¯ ( )
has the same shape as the Starobinsky potential and has
the minimum for = 1 + 8γ λ.
– While the freeze double singularities appear in our model
in the Jordan frame there are no such singularities in the
dynamics of the model in the Einstein frame.
1/3
– If
γ is small, then asing =
for
negative γ and asing = 8 ,0+1−γ ( ,m01+ ,0) for positive
γ . These values define the natural scale at which
singularities appear in the model under consideration with the
negative or positive value of γ parameter. It seems to be
natural to identify this scale with a cut off at which the
model can be treated as some kind of effective theory.
– In both the cases of a negative and positive γ one deals
with a finite scale factor singularity. For negative γ it
is a double sudden singularity which meets the future
singularity of a contracting model before the bounce with
the initial singularity in the expanding model. The sewn
evolutionary scenarios reveal the presence of a bounce
during the cosmic evolution.
– In the context of the Starobinsky model in the Palatini
formalism we found a new type of double singularity beyond
the well-known classification of isolated singularities.
– The phase portrait for the model with a positive value
of γ is equivalent to the phase portrait of the CDM
model (following dynamical system theory [80]
equivalence assumes the form of topological equivalence
established by a homeomorphism). There is only a quantitative
difference related with the presence of the non-isolated
freeze singularity. The scale of the appearance of this type
singularity can also be estimated and be cast in terms of
the redshift zfreeze = γ−1/3.
– We estimated the model parameters using astronomical
data and conclude that positive γ is favored by the best
fit value; still the model without Rˆ 2 term is statistically
admitted.
In our model, the best fit value of γ is equal 9.70 ×
10−11 and positive γ parameter belongs to the interval
(0, 2.2143 × 10−9) at 2-σ level. This mean that the
positive value of γ is more favored by astronomical data than
the negative value of γ . The difference between values of
BIC for our model and the CDM model is equal 6.407.
So, in comparison to our model, the evidence in favor of the
CDM model is strong. But one cannot absolutely reject the
model.
Note added in proof After completing the paper we found
a paper by Faraoni and Cardini where freeze singularities
have been analyzed in a different context, both from point
particle and cosmological perspectives [81].
Acknowledgements AB and MS would like to thank Salvatore
Capozziello for encouraging us to consider this model in different
frames. The work has been supported by the Polish National Science
Centre (NCN), project DEC-2013/09/B/ ST2/03455.
Open Access This article is distributed under the terms of the Creative
Commons Attribution 4.0 International License (http://creativecomm
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit
to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made.
Funded by SCOAP3.
1. A. De Felice , S. Tsujikawa, f(R) Theories . Living Rev. Rel . 13 , 3 ( 2010 ). arXiv: 1002 . 4928
2. A.A. Starobinsky , A new type of isotropic cosmological models without singularity . Phys. Lett. B 91 , 99 - 102 ( 1980 )
3. A.H. Guth , The inflationary universe: a possible solution to the horizon and flatness problems . Phys. Rev. D 23 , 347 - 356 ( 1981 )
4. A.D. Linde , Eternally existing selfreproducing chaotic inflationary universe . Phys. Lett. B 175 , 395 - 400 ( 1986 )
5. S. Nojiri , S.D. Odintsov , Introduction to modified gravity and gravitational alternative for dark energy . Int. J. Geom. Method Mod. Phys. 4 , 115 ( 2007 ). arXiv: hep-th/0601213
6. T.P. Sotiriou , V. Faraoni , f(R) Theories of gravity . Rev. Mod. Phys . 82 , 451 - 497 ( 2010 ). arXiv: 0805 . 1726
7. V.F. Mukhanov , G.V. Chibisov , Quantum fluctuations and a nonsingular universe . JETP Lett . 33 , 532 - 535 ( 1981 ). Pisma Zh. Eksp. Teor. Fiz . 33 , 549 ( 1981 )
8. A.A. Starobinsky , The perturbation spectrum evolving from a nonsingular initially de-Sitter cosmology and the microwave background anisotropy . Sov. Astron. Lett. 9 , 302 ( 1983 )
9. Planck Collaboration , P.A.R. Ade , Planck 2015 results . XIV. Dark energy and modified gravity . Astron. Astrophys . 594 , A14 ( 2016 ). arXiv: 1502 . 01590
10. C. Cheng, Q.-G. Huang, Y.-Z. Ma , Constraints on single-field inflation with WMAP. SPT and ACT data-a last-minute stand before Planck . JCAP 1307 , 018 ( 2013 ). arXiv: 1303 . 4497
11. Q.-G. Huang , A polynomial f(R) inflation model . JCAP 1402 , 035 ( 2014 ). arXiv: 1309 . 3514
12. L.A. Kofman , A.D. Linde , A.A. Starobinsky , Inflationary universe generated by the combined action of a scalar field and gravitational vacuum polarization . Phys. Lett. B 157 , 361 - 367 ( 1985 )
13. S.V. Ketov , A.A. Starobinsky , Embedding (R + R2 ) -inflation into supergravity . Phys. Rev. D 83 , 063512 ( 2011 ). arXiv: 1011 . 0240
14. S.A. Appleby , R.A. Battye , A.A. Starobinsky , Curing singularities in cosmological evolution of F(R) gravity . JCAP 1006 , 005 ( 2010 ). arXiv: 0909 . 1737
15. S. Capozziello , M. De Laurentis , S. Nojiri , S.D. Odintsov , Classifying and avoiding singularities in the alternative gravity dark energy models . Phys. Rev. D 79 , 124007 ( 2009 ). arXiv: 0903 . 2753
16. A. Alho , S. Carloni , C. Uggla , On dynamical systems approaches and methods in f (R) cosmology . JCAP 1608 ( 08 ), 064 ( 2016 ). arXiv: 1607 . 05715
17. S. Capozziello , M.F. De Laurentis , L. Fatibene , M. Ferraris , S. Garruto , Extended cosmologies . SIGMA 12 , 006 ( 2016 ). arXiv: 1509 . 08008
18. S. Capozziello, P. Martin-Moruno , C. Rubano , Physical nonequivalence of the Jordan and Einstein frames . Phys. Lett. B 689 , 117 - 121 ( 2010 ). arXiv: 1003 . 5394
19. S.M. Carroll , A. De Felice , V. Duvvuri , D.A. Easson , M. Trodden , M.S. Turner , The cosmology of generalized modified gravity models . Phys. Rev. D 71 , 063513 ( 2005 ). arXiv: astro-ph/0410031
20. A. Borowiec , M. Kamionka , A. Kurek , M. Szydlowski , Cosmic acceleration from modified gravity with Palatini formalism . JCAP 1202 , 027 ( 2012 ). arXiv: 1109 . 3420
21. G.J. Olmo , Palatini approach to modified gravity: f(R) theories and beyond . Int. J. Mod. Phys. D 20 , 413 - 462 ( 2011 ). arXiv: 1101 . 3864
22. G.J. Olmo , Post-Newtonian constraints on f(R) cosmologies in metric and Palatini formalism . Phys. Rev. D 72 , 083505 ( 2005 ). arXiv: gr-qc/0505135
23. G.J. Olmo , The gravity Lagrangian according to solar system experiments . Phys. Rev. Lett . 95 , 261102 ( 2005 ). arXiv: gr-qc/0505101
24. C. Barragan , G.J. Olmo , Isotropic and anisotropic bouncing cosmologies in Palatini gravity . Phys. Rev. D 82 , 084015 ( 2010 ). arXiv: 1005 . 4136
25. C. Barragan , G.J. Olmo , H. Sanchis-Alepuz , Bouncing cosmologies in Palatini f(R) gravity . Phys. Rev. D 80 , 024016 ( 2009 ). arXiv: 0907 . 0318
26. C. Bejarano , G.J. Olmo , D. Rubiera-Garcia, What is a singular black hole beyond general relativity? Phys. Rev. D 95 , 064043 ( 2017 ). arXiv: 1702 . 01292
27. C. Bambi , A. Cardenas-Avendano , G.J. Olmo , D. Rubiera-Garcia, Wormholes and nonsingular spacetimes in Palatini f (R) gravity . Phys. Rev. D 93 , 064016 ( 2016 ). arXiv: 1511 . 03755
28. G.J. Olmo , D. Rubiera-Garcia, Nonsingular black holes in f (R) theories . Universe 1 , 173 - 185 ( 2015 ). arXiv: 1509 . 02430
29. G.J. Olmo , D. Rubiera-Garcia, Nonsingular black holes in quadratic Palatini gravity . Eur. Phys. J. C 72 , 2098 ( 2012 ). arXiv: 1112 . 0475
30. G.J. Olmo , D. Rubiera-Garcia, Palatini f (R) black holes in nonlinear electrodynamics . Phys. Rev. D 84 , 124059 ( 2011 ). arXiv: 1110 . 0850
31. E.E. Flanagan , Palatini form of 1/R gravity . Phys. Rev. Lett . 92 , 071101 ( 2004 ). arXiv: astro-ph/0308111
32. E.E. Flanagan, The conformal frame freedom in theories of gravitation . Class. Quantum Gravity 21 , 3817 ( 2004 ). arXiv: gr-qc/0403063
33. F.A. Teppa Pannia , F. Garcia , S.E. Perez Berliaffa , M. Orellana , G.E. Romero , Structure of compact stars in R-squared Palatini gravity . Gen. Relativ. Gravity 49 , 25 ( 2017 ). arXiv: 1607 . 03508
34. T. Koivisto, Covariant conservation of energy momentum in modified gravities . Class. Quantum Gravity 23 , 4289 - 4296 ( 2006 ). arXiv: gr-qc/0505128
35. L. Fernandez-Jambrina , R. Lazkoz , Classification of cosmological milestones . Phys. Rev. D 74 , 064030 ( 2006 ). arXiv: gr-qc/0607073
36. S. Bahamonde , S.D. Odintsov , V.K. Oikonomou , P.V. Tretyakov , Deceleration versus acceleration universe in different frames of F (R) gravity . Phys. Lett. B 766 , 225 - 230 ( 2017 ). arXiv: 1701 . 02381
37. M.P. Dabrowski , J. Garecki , D.B. Blaschke , Conformal transformations and conformal invariance in gravitation . Ann. Phys. 18 , 13 - 32 ( 2009 ). arXiv: 0806 . 2683
38. BICEP2 Collaboration, Ade, P.A.R. , et al.: Detection of B-mode polarization at Degree Angular Scales by BICEP2 . Phys. Rev. Lett . 112 ( 24 ), 241101 ( 2014 ). arXiv: 1403 . 3985
39. O. Hrycyna , M. Szydlowski , M. Kamionka , Dynamics and cosmological constraints on Brans-Dicke cosmology . Phys. Rev. D 90 ( 12 ), 124040 ( 2014 ). arXiv: 1404 . 7112
40. S. Nojiri , S.D. Odintsov , S. Tsujikawa , Properties of singularities in (phantom) dark energy universe . Phys. Rev. D 71 , 063004 ( 2005 ). arXiv: hep-th/0501025
41. J.D. Barrow , Sudden future singularities . Class. Quantum Gravity 21 , L79 - L82 ( 2004 ). arXiv: gr-qc/0403084
42. M. Bouhmadi-Lopez , P.F. Gonzalez-Diaz , P. Martin-Moruno , Worse than a big rip? Phys . Lett. B 659 , 1 - 5 ( 2008 ). arXiv:gr-qc/0612135
43. A. Krolak , Towards the proof of the cosmic censorship hypothesis . Class. Quantum Gravity 3 , 267 - 280 ( 1986 )
44. A. Borowiec , A. Stachowski , M. Szydlowski , Inflationary cosmology with Chaplygin gas in Palatini formalism . JCAP 1601 ( 01 ), 040 ( 2016 ). arXiv: 1512 . 01199
45. Sbierski , J.: The C0-inextendibility of the Schwarzschild spacetime and the spacelike diameter in Lorentzian geometry ( 2015 ). arXiv: 1507 . 00601
46. Galloway , G. , Ling , E. : Some remarks on the C0-(in)extendibility of spacetimes ( 2016 ). arXiv: 1610 . 03008
47. O. Hrycyna , M. Szydlowski , Non-minimally coupled scalar field cosmology on the phase plane . JCAP 0904 , 026 ( 2009 ). arXiv: 0812 . 5096
48. G .F .R. Ellis , E. Platts , D. Sloan , Current observations with a decaying cosmological constant allow for chaotic cyclic cosmology . JCAP 1604 ( 04 ), 026 ( 2016 ). arXiv: 1511 . 03076
49. Bautin , N.N. , Leontovich , I.A . (eds.): Methods and Techniques for Qualitative Analysis of Dynamical Systems on the Plane . Nauka, Moscow ( 1976 ). (In Russian)
50. N. Suzuki et al., The Hubble space telescope cluster supernova survey: V. Improving the dark energy constraints above z>1 and building an early-type-hosted supernova sample . Astrophys. J . 746 , 85 ( 2012 ). arXiv: 1105 . 3470
51. S.D.S.S. Collaboration , W.J. Percival et al., Baryon acoustic oscillations in the sloan digital sky survey data release 7 galaxy sample . Mon. Not. R. Astron. Soc . 401 , 2148 - 2168 ( 2010 ). arXiv: 0907 . 1660
52. F. Beutler , C. Blake , M. Colless , D.H. Jones , L. Staveley-Smith , L. Campbell , Q. Parker , W. Saunders , F. Watson , The 6dF galaxy survey: baryon acoustic oscillations and the local Hubble constant . Mon. Not. R. Astron. Soc . 416 , 3017 - 3032 ( 2011 ). arXiv: 1106 . 3366
53. C. Blake et al., The WiggleZ dark energy survey: Joint measurements of the expansion and growth history at z<1 . Mon. Not. R. Astron . Soc. 425 , 405 - 414 ( 2012 ). arXiv: 1204 . 3674
54. W. Hu, N. Sugiyama , Small scale cosmological perturbations: an analytic approach . Astrophys. J. 471 , 542 - 570 ( 1996 ). arXiv: astro-ph/9510117
55. D.J. Eisenstein , W. Hu, Baryonic features in the matter transfer function . Astrophys. J . 496 , 605 ( 1998 ). arXiv: astro-ph/9709112
56. C. Alcock , B. Paczynski , An evolution free test for non-zero cosmological constant . Nature 281 , 358 - 359 ( 1979 )
57. M. Lopez-Corredoira , Alcock-Paczynski cosmological test . Astrophys. J . 781 ( 2 ), 96 ( 2014 ). arXiv: 1312 . 0003
58. P.M. Sutter , G. Lavaux , B.D. Wandelt , D.H. Weinberg , A first application of the Alcock-Paczynski test to stacked cosmic voids . Astrophys. J . 761 , 187 ( 2012 ). arXiv: 1208 . 1058
59. C. Blake et al., The WiggleZ dark energy survey: measuring the cosmic expansion history using the Alcock-Paczynski test and distant supernovae . Mon. Not. R. Astron. Soc . 418 , 1725 - 1735 ( 2011 ). arXiv: 1108 . 2637
60. N.P. Ross et al., The 2dF-SDSS LRG and QSO survey: the 2-point correlation function and redshift-space distortions . Mon. Not. R. Astron. Soc . 381 , 573 - 588 ( 2007 ). arXiv: astro-ph/0612400
61. C. Marinoni , A. Buzzi , A geometric measure of dark energy with pairs of galaxies . Nature 468 ( 7323 ), 539 - 541 ( 2010 )
62. J. da Angela, P.J. Outram , T. Shanks, Constraining beta(z) and Omega 0(m) from redshift-space distortions in z 3 galaxy surveys . Mon. Not. R. Astron. Soc . 361 , 879 - 886 ( 2005 ). arXiv: astro-ph/0505469
63. P.J. Outram , T. Shanks , B.J. Boyle , S.M. Croom , F. Hoyle , N.S. Loaring , L. Miller , R.J. Smith, The 2df qso redshift survey. 13. A measurement of lambda from the qso power spectrum . Mon. Not. R. Astron. Soc . 348 , 745 ( 2004 ). arXiv: astro-ph/0310873
64. L. Anderson et al., The clustering of galaxies in the SDSS-III baryon oscillation spectroscopic survey: baryon acoustic oscillations in the data release 9 spectroscopic galaxy sample . Mon. Not. R. Astron. Soc . 427 ( 4 ), 3435 - 3467 ( 2013 ). arXiv: 1203 . 6594
65. I. Paris et al., The Sloan Digital Sky Survey quasar catalog: ninth data release . Astron. Astrophys . 548 , A66 ( 2012 ). arXiv: 1210 . 5166
66. S.D.S.S. Collaboration , D.P. Schneider et al., The Sloan Digital Sky Survey quasar catalog V. Seventh data release . Astron. J . 139 , 2360 - 2373 ( 2010 ). arXiv: 1004 . 1167
67. J. Simon , L. Verde , R. Jimenez , Constraints on the redshift dependence of the dark energy potential . Phys. Rev. D 71 , 123001 ( 2005 ). arXiv: astro-ph/0412269
68. D. Stern , R. Jimenez , L. Verde , M. Kamionkowski , S.A. Stanford , Cosmic chronometers: constraining the equation of state of dark energy. I: H(z) measurements . JCAP 1002 , 008 ( 2010 ). arXiv: 0907 . 3149
69. M. Moresco et al., Improved constraints on the expansion rate of the Universe up to z 1.1 from the spectroscopic evolution of cosmic chronometers . JCAP 1208 , 006 ( 2012 ). arXiv: 1201 . 3609
70. N. Metropolis , A.W. Rosenbluth , M.N. Rosenbluth , A.H. Teller , E. Teller, Equation of state calculations by fast computing machines . J. Chem. Phys . 21 , 1087 - 1092 ( 1953 )
71. W.K. Hastings , Monte carlo sampling methods using markov chains and their applications . Biometrika 57 , 97 - 109 ( 1970 )
72. G. Schwarz, Estimating the dimension of a model . Ann. Stat . 6 ( 2 ), 461 - 464 ( 1978 )
73. R.E. Kass , A.E. Raftery , Bayes factors . J. Am. Stat. Assoc . 90 , 773 - 795 ( 1995 )
74. J.D. Barrow , A.A.H. Graham , Singular inflation . Phys. Rev. D 91 ( 8 ), 083513 ( 2015 ). arXiv: 1501 . 04090
75. Bamba , K. : Inflationary universe in fluid description ( 2016 ). arXiv: 1601 . 04773
76. S. Nojiri , S.D. Odintsov , V.K. Oikonomou , Singular inflation from generalized equation of state fluids . Phys. Lett. B 747 , 310 - 320 ( 2015 ). arXiv: 1506 . 03307
77. S.D. Odintsov , V.K. Oikonomou , Singular inflationary universe from F (R) gravity . Phys. Rev. D 92 ( 12 ), 124024 ( 2015 ). arXiv: 1510 . 04333
78. Z.T. Zhusubaliyev , E. Mosekilde, Bifurcations and Chaos in Piecewise-Smooth Dynamical Systems (World Scientific, Singapore , 2003 )
79. R. Leine , H. Nijmeijer , Dynamics and Bifurcations of Non-Smooth Mechanical Systems (Springer, Berlin, 2004 )
80. L. Perko et al., Differential Equations and Dynamical Systems , vol. 7 , 3rd edn., Texts in Applied Mathematics (Springer, New York, 2001 )
81. Faraoni , V. , Cardini , A.M.: Analogues of glacial valley profiles in particle mechanics and in cosmology ( 2016 ). arXiv: 1608 . 02542