Starobinsky cosmological model in Palatini formalism

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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f({\hat{R}})={\hat{R}}+\gamma {\hat{R}}^2$$\end{document} 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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _\gamma =3\gamma H_0^2$$\end{document} is equal 9.70×10-11\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$9.70\times 10^{-11}$$\end{document}. 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.


Introduction
The main aim of the paper is the construction of the Starobinsky model with a squared termR 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 R 2 /6M 2 was motivated by the conformal anomaly in the quantum gravity. The problem of a e-mail: aleksander.stachowski@uj.edu.pl b e-mail: marek.szydlowski@uj.edu.pl c e-mail: andrzej.borowiec@ift.uni.wroc.pl 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ä(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 n s = 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 + R 2 6M 2 + λ n 2n R n (3M 2 ) n−1 . It was demonstrated that the slow-roll inflation can be achieved as long as the dimensionless coupling λ n is sufficiently small [11].
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][17][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][27][28][29][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 whereR = 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 where T μν = − 2 √ −g δL m δg μν is matter energy-momentum tensor, 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, (4) for which the connection = LC (h) is a Levi-Civita connection. As a consequence in dim M = 4 one gets 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 where T = g μν T μν . Because of cosmological applications we assume that the metric g is FRW metric 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 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 model by introducing a curvature fluid ρ k = − k 2 a −2 , with the barotropic factor w = − 1 3 ( p k = − 1 3 ρ k ). From the conservation condition T μ ν;μ = 0 we obtain ρ = ρ 0 a −3(1+w) . Therefore the trace T reads 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, induced by the first three terms in the power series decomposition of an arbitrary function f (R 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) 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: where r,0 = ρ r,0 ,0 = ρ ,0 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.

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] S(g μν , λ ρσ , χ) = Introducing a scalar field = f (χ ) and taking into account the constraint χ =R, one gets the action (22) in the following form: where the potential U ( ) is defined by The Palatini variation of the action (23) gives rise to the following equations of motion: Equation (25b) implies that the connectionˆ is a metric connection for a new metricḡ μν = g μν ; thusR μν =R μν ,R = g μνR μν = −1R andḡ μνR = g μνR . The g-trace of (25a) produces a new structural equation Now Eqs. (25a) and (25c) take the following form: where we introduceŪ (φ) = U (φ)/ 2 ,T μν = −1 T μν and the structural equation can be replaced by In this case, the action for the metricḡ μν and scalar field is given by where we have to take into account a non-minimal coupling between andḡ μν [17,37]). In FRW case, one gets the metricḡ μν in the following form: where dt = (t) 1 2 d t and new scale factorā(t) = (t) 1 2 a(t). Ensuing cosmological equations (in the case of the barotropic matter) are given by and w =p m /ρ m = p m /ρ m . In this case, the conservation equations has the following form: Let us consider the Starobinsky-Palatini model in the above framework. The potentialŪ is described by the following formula: Figure 1 presents the relationρ ( ). Note that the functionρ has the same shape like the Starobinsky potential. The functionρ ( ) has the minimum for In general, the scalar field (ā) is given by (cf. (11)) Becauseρ m = −2 ρ m ,p m = −2 p m , and taking into account (34) one gets the algebraic equation determining the function (ā) for a given barotropic factor w. This provides an implicit dependence (ā). 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 where y = − 1 2 . The evolution of (ā), at the beginning of the inflation epoch, is presented in Fig. 2.
For γ ≈ 0, the potentialŪ can be approximated asŪ = −ρ m + 1 4γ . In this case the Friedmann equation can be written as In the case ofρ m = 0,ρ is constant and the Friedmann equation has the following form: 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: 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 s 2 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. The condition for appearing of the inflation is for the value of the parameter γ to be close to zero, hence the influence of the parameter λ on the evolution of the universe is negligible.
In Fig. 4 the typical evolution is demonstrated ofρ m (ā) at the beginning of the inflation epoch. The typical evolution ofρ , at the beginning of the inflation epoch, is presented in Fig. 5. Note that, for the late time universe,ρ can be approximated as a constant. Figure 6 presents the evolution of the scale factorā(t) during the inflation. Figure 7 shows the Hubble functionH during the inflation epoch.
The conservation equation forρ can be writteṅ wherep is an effective pressure. In this case the equation of state for the dark energy is expressed by the following formula: where the function w(a) is given by the expression 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 where R is the metric Ricci scalar, = f (R),R = χ ( ).
We obtain the Brans-Dicke action with the coupling parameter ω = − 3 2 in the Jordan frame. The equations of motion take the form In this case the dynamics of the metric g is exactly the same as described by the original Palatini equations (2)-(6). On cosmological grounds it means that the scale factor a(t) evolves according to the Friedmann equation (13). It has recently been shown that cosmological data favor the value ω ≈ −1 on the 3σ level [39].

Singularities in cosmological dynamical systems of Newtonian type
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 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, where dot denotes differentiation with respect to the cosmic time t, H ≡ȧ a 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 tȯ It would be convenient to rewrite (49) in the equivalent forṁ where In (53) ρ(a) plays the model role of an effective energy density. For example for the standard cosmological model (12) where ρ eff = ρ m + and ρ m = ρ m,0 a −3 . Equation (50) is equivalent tö which is called the acceleration equation. It is easily to check thaẗ 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,ȧ). 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, In the general case of an arbitrary potential, the dynamical system which describes the evolution of a universe takes the forṁ 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 x 2 2 + V (a) = − k 2 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 t → τ = |H 0 |t and 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 x 0 = 0, ∂ V (a) ∂a | a=a 0 . 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 where det A is the determinant of the linearization matrix calculated at the critical points, i.e.
From Eqs. (64) and (65) ∂a 2 | a=a 0 > 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, the deceleration parameter, the effective barotropic factor the parameter of deviation from de Sitter universe [35] 6 , h(t) = 0), the effective matter density and pressure and, finally, the Ricci scalar curvature for the FRW metric (48), From the formulas above one can observe that the most of them depend on the quantity 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) ∝ a 2 , 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.

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Ḣ 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Ḣ 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 a s can appear for which ∂ V ∂a (a s ) 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 dt da = 1 Ha is zero. For the freeze singularity, the scale factor remains constant a s , ρ and H blow up andä is undefined. It this case, the degenerate singularity of type III is called sewn (non-isolated) singularity. We, therefore, obtain [44] t − t s ± 1 2 The above formula combines two types of behavior near the freeze singularities in the future, a − a sing ∝ −(t sing − t) 1/2 for t → t sing − (77) and in the past 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 < a s } and {a : a > a s }. 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 = a min . In this singularity the potential itself is a continuous func-tion 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].

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 multiplicative expression b 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 model when the multiplicative expression b 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 where K ∈ [0, 3). The solution of the above equation is From Eq. (81), we can find an expression for a value of the scale factor for the freeze singularity (82) The relation between a freeze and positive γ is presented in Fig. 11.
The sudden singularity appears when b = 0. This leaves us with the following algebraic equation: which, in fact, becomes a (degenerate) critical point and a bounce at the same time. The relation between a sing and negative γ is presented in Fig. 12.
+ ch + k . We can rewrite dynamical system (58)-(59) as  (1) and (2) present two static Einstein universes. The phase portrait belongs to the class of sewn dynamical systems [49] where ≡ d dσ = b+ d 2 b d dτ 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 < a sing and the second part is for a > a sing . Both parts are glued along the singularity.
For a > a sing , in an analogous way, we get the following equations: where  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 t init to t fin of this inflation period N = H init (t fin − t init ) (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. 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-Paczyń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: 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: where r s (z d ) 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: where AP(z) th ≡ H (z) z z 0 dz H (z ) and AP(z i ) obs are observational data [58][59][60][61][62][63][64][65][66].
The likelihood function for measurements of the Hubble parameter H (z) of galaxies from [67][68][69] is given by the expression 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: where C is the covariance matrix with the errors, x is a vector of the acoustic scale l A , the shift parameter R and b h 2 where where z * is the redshift of the epoch of the recombination [54]. The total likelihood function is expressed in the following form: To estimate model parameters, we use our own code Cos-moDarkBox. The Metropolis-Hastings algorithm [70,71] is used in this code.   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 ( γ , H 0 ) 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 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.  [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.

Conclusions
In this paper, we demonstrated that evolution of the Starobinsky model with a quadratic term R 2 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][75][76][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 = a sing . 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 < a sing }. 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 RR > 0 gives rise to the negative values of M 2 . We do 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Ū ( ) has the same shape as the Starobinsky potential and has the minimum for = 1 + 8γ λ. 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 z freeze = −1/3 γ . -We estimated the model parameters using astronomical data and conclude that positive γ is favored by the best fit value; still the model withoutR 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].