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 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(\hat{R})=\hat{R}+\gamma \hat{R}^2$ in the Palatini formalism, when dynamics is determined by the corresponding piece-wise smooth dynamical system. As an effect we obtain a degenerated 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 $\Omega_\gamma=3\gamma H_0^2$ is equal $9.70\times 10^{-11}$. We consider model in both Jordan and Einstein frames. We show that after transition to the Einstein frame we obtain both 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 a aleksander.stachowski@uj.edu.pl b marek.szydlowski@uj.edu.pl c andrzej.borowiec@ift.uni.wroc.pl Starobinsky [2]. In the original Starobinsky model the term R 2 /6M 2 was motivated by the conformal anomaly in the quantum gravity. Problem of inflation in f (R) cosmological model are strictly related with the choice of frames. Authors [1] show that CMB spectra in the 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 a 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 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]. 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][13][14][15][16][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 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 the modern cosmology, the most popular trend is to explain the dark energy and the dark matter in terms of some substances, which nature is unknown up to now. Albert Einstein was representing the opposite relational point of view on description of gravity, in which all substantial forms should be eliminated. Such a point of view is called anti-substantialism. 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 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.
In the Palatini gravity action for f (R) gravity is introduced to be whereR = g µνR µν (Γ ) is the generalized Ricci scalar andR µν (Γ ) 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 δLm δ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 [21]. In eq. (3)∇ α means the covariant derivative calculated with respect to Γ . In order to solve equation (3) it is convenient to introduce a new metric for which the connection Γ = Γ L−C (h) is a Levi-Civita connection. As a consequence in dim M = 4 one gets i.e. that 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 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 be also included to the model by introducing a curvature fluid . From the conservation condition T µ ν;µ = 0 we obtain that ρ = ρ 0 a −3(1+w) . Therefore trace T reads as 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 first three terms in the power series decomposition of an arbitrary function f (R). In fact, since the termsR 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, whereR 0 is a constant and [R 0 ] = [R]. In this case the power series expansion reads: 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 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 section 2, we consider the Palatini approach in the Jordan and Einstein frame. In section 3, we present some generalities concerning dynamical systems of Newtonian type, and their relations with the Palatini-Starobinsky model. Section 4, is devoted to 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 [22] to these systems using the notion of elasticity of the potential function with respect the scale factor. In section 5, we will analyze the singularities in the Starobinsky model in the Palatini formalism. This system requires the form of piece-wise smooth dynamical system. Statistical analysis of the model is presented in section 6. In section 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 [23] 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 =ḡ µνR µν = Φ −1R andḡ µνR = g µνR . The gtrace of (25a) produces a new structural equation Now equations (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,24]).
In FRW case, one gets the metricḡ µν in the following form where dt = Φ(t)  Ensuing cosmological equations (in the case of the barotropic matter) are given by wherē and w =p m /ρ m = p m /ρ m . In this case, the conservation equations has the following forṁ Let us consider the Starobinsky-Palatini model in the above framework. The potentialŪ is described by the following formulā 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)) 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 The evolution of Φ(ā), at the beginning of the inflation epoch, is presented in figure 2.
For γ ≈ 0, the potentialŪ can be approximated as U = −ρ 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 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 [25]. The relation between γ and the approximate number of e-foldings N is presented in figure 3.
The condition for appearing of the inflation requires the value of the parameter γ to be close to zero, hence the influence of the parameter λ for the evolution of the universe is negligible.
In figure 4 it is demonstrated the typical evolution ofρ m (ā) at the beginning of the inflation epoch. The typical evolution ofρ Φ , at the beginning of the inflation epoch, is presented in figure 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 function H during the inflation epoch.
The conservation equation forρ Φ can be written aṡ wherep Φ is an effective pressure. In this case the equation of state for the dark energy is expressed by the following formulā where the function w(a) is given by the expression The diagram of coefficient of equation of state w(a), at the beginning the inflation epoch, is presented in figure  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, We obtain Brans-Dicke action with the coupling parameter ω = − 3 2 in the Jordan frame. 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 the cosmological ground it means that the scale factor a(t) evolves according to the Friedmann equation (13). It has been recently shown that cosmological data favor the value ω ≈ −1 on the 3σ level [26].

Singularities in cosmological dynamical systems of Newtonian type
There is a class of cosmological models, which dynamics can be reduced to the 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 assumes the form of 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 that ρ(t) = ρ(a(t)) as well as p(t) = p(a(t)), i.e. that both energy density as well as pressure depends 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 an equivalent forṁ where In (53) ρ(a) plays the model role of 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 acceleration equation. It is easily to check thaẗ where V (a) is given by (53) provided that conservation equation (51) is fulfilled. Due to equation (56) the evolution of a universe can be interpreted as a motion of a fictitious particle of unit mass in the potential V (a). Here a(t) plays the role of a positional variable. Equation of motion (56) assumes the form analogous to the Newtonian equation of motion.
If we know the form of effective energy density then we can construct the form of potential V (a), which determine 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 potential V (a) contains all information needed to construction of a phase space portrait. In this case the phase space is two-dimensional (a,ȧ) :ȧ In a general case of arbitrary potential, a dynamical system which describes the evolution of a universe takes the forṁ We shall study the system above using 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. The equations (58) and (59) can be rewritten in dimensionless variables if we replace effective energy density ρ eff by 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=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 where det A is determinant of linearization matrix calculated at the critical points, i.e.
From equation (64) and (65) one can conclude that only admissible critical points are the saddle type if ∂a 2 | a=a0 < 0 or centers type if ∂ 2 V (a) ∂a 2 | a=a0 > 0. If a shape of the potential function is known (from the knowledge of effective energy density), then it is possible to calculate cosmological functions in exact form a deceleration parameter, an effective barotropic factor a parameter of deviation from de Sitter universe [22] h(t) ≡ −(q(t) + 1) = 1 2 (note that if V (a) = − Λa 2 6 , h(t) = 0), effective matter density and pressure and, finally, a 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 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 [22] the crucial role is played by the parameter h(t). Note that a cosmological sense of 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 discovery of an accelerated phase in the expansion of our Universe, many theoretical possibilities for future singularity are seriously considered. If we assume that the Universe expands following the Friedmann equation, then this phase of expansion is driven by dark energy -hypothetical fluid, which violates the strong energy condition. Many of new types of singularities were classified by Nojiri et al. [27]. 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 ρ:  [28] or divergent [29]. -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 [30], 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 or/and its derivative can diverge at isolated points (value of the scale factor). Therefore mentioned before classification has application only for a single component of piece-wise smooth potential. In our model the dynamical system describing evolution of a universe belongs to the class of a piecewise smooth dynamical sys-tems. 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 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 degenerated. Formally it is composed of two types III singularities: one in the future and another one in the past. If we considered the evolution of the universe before this singularity we detect isolated singularity of type III in the future. Conversely if we consider the evolution after the singularity then going back in time we meet type III singularity in the past. Finally, at the finite scale factor both singularities will meet together. For description of behavior near the singularity one considers t = t(a) relation. This relation has a horizontal inflection point and it is natural to expand this relation in the 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 a is undefined. It this case, the degenerated singularity of type III is called sewn (non-isolated) singularity. We, therefore, obtain [31] t − t s ± 1 2 Above formula combine two types of behavior near the freeze singularities in the future a − a sing ∝ −(t sing − t) 1/2 for t → t sing − and in the past a − a sing ∝ +(t − t sing ) 1/2 for t → t sing + .
(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 < 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 function while its first derivative has a discontinuity. Therefore, the corresponding dynamical system represents a piece-wise smooth 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 Jan Sbierski [32]. In the context of FLRW cosmological models, Sbierski's methodology was considered in [33].

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 the 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 (81) From equation (81), we can find an expression for a value of the scale factor for the freeze singularity The relation between a freeze and positive Ω γ is presented in Figure 11.
The sudden singularity appears when b = 0. This provides the following algebraic equation The above equation can be rewritten as 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 Figure 12.
+ Ω ch + Ω k . We can rewrite dynamical system (58)-(59) as where We can treated the dynamical system (86)-(87) as a sewn dynamical system [34,35]. 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 , dynamical system (86)-(87) can be rewritten to the corresponding form where V 1 = V (−η(a − a s ) + 1) and η(a) notes the Heaviside function. For a > a sing , in the analogous way, we get the following equations where V 2 = V η(a − a s ). The phase portraits, for dynamical system (86)-(87), are presented in Figure 13 and 15. Figure 13 provides the phase portrait for positive Ω γ while Figure 15 provides the phase portrait for negative Ω γ . In Figure 13 there are two critical points labelled '1' and '2' at the finite domain. They are both the saddles. These critical points correspond to a maximum of the potential function. The saddle '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 to 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. Diagram of the evolution of scale factor for closed orbit is presented by Figure 14. It is also interesting that trajectories in 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 formula  Fig. 11 Diagram of the relation between a sing and positive Ωγ . Note that in the limit Ωγ → 0 the singularity overlaps with a big-bang singularity.  Fig. 12 Diagram of the relation between a sing and negative Ωγ . Note that in the limit Ωγ → 0 the singularity overlaps with a big-bang singularity. (3.13) in [1]) with respect to Ω γ is shown in Figure 16. This figure illustrates the number of e-foldings is too small for to obtain the inflation effect.

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-Paczyński test, measurements of CMB and lensing by Planck and low by WMAP. The likelihood function for observations of supernovae of type Ia [37] is given by the following expression 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 log 10 D L + 25 (where the luminosity distance is D L = c(1 + z) 6dF Galaxy Redshift Survey measurements at redshift z = 0.1 [39], and WiggleZ measurements at redshift z = 0.44, 0.60, 0.73 [40] have the following likelihood function where r s (z d ) is the sound horizon at the drag epoch [41,42]. For the Alcock-Paczynski test [43,44] 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 [45][46][47][48][49][50][51][52][53].
The likelihood function for measurements of the Hubble parameter H(z) of galaxies from [54][55][56] is given by the expression In this paper, we use the likelihood function for observations of CMB [9] and lensing by Planck, and lowpolarization from the WMAP (WP) in the following 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 [41]. The total likelihood function is expressed in the following form To estimate model parameters, we use our own code CosmoDarkBox. The Metropolis-Hastings algorithm [57,58] is used in this code. Table 1 completes the values of parameters for the best fit with errors. Figure 17 and Figure 18 show the intersection of a likelihood function with the 68% and 95% confidence level projections on the planes (Ω γ , Ω m,0 ) and (Ω γ , H 0 ).
In this paper, we use Bayesian information criterion (BIC) [59,60], for comparison 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 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=BIC-BIC ΛCDM is equal 6.407. The evidence for the model is strong [60] if ∆BIC is more 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 piecewisesmooth dynamical systems (PWS), i.e., systems whose the phase space is partitioned into different regions, each of them associated to a different smooth functional 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 H 0 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. 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 the inflation as an endogenous effect of the dynamics in the Palatini formalism. While the effect of inflation appears in the model under consideration the number of sufficient e-folds is not achieved and the additional effect of amplification is required. Note that this type of inflation is a realisation of the idea of singular inflation [61][62][63][64]. In our model inflation is driven by the freeze degenerated singularity (the extension of type III isolated singularity). We show that 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 mimicking, the evolution of the universe contain all information which are needed for 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 called this type of sewn singularities in analogy to sewn dynamical systems [65,66].
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. Figure 13). In the first domain the evolution starts from big-bang followed by the deceleration phase then 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 correspond to the present day epoch called 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 the 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 know 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 section 7.4 in [1]). In our case, statistical analysis more favors model with f R > 0 (Ω γ > 0) than with f R < 0 (Ω γ < 0). In other words, statistical analysis more 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 obtain 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 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 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 similarly as the original (metric) Starobinsky model is very important for the explanation of the 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.
-The 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 like the Starobinsky potential and has the minimum for Φ = 1 + 8γλ. -While the freeze double singularities appear in our model in the Jordan frame there is no such singularities in the dynamics of the model in the Einstein frame.
-If Ω γ is small, then a sing = − for positive Ω γ . These values defines 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 cases of 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 contracting model before the bounce with the initial singularity in the expanding model. The sewn evolutionary scenarios reveal the presence of 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 the dynamical system theory [67] equivalence assumes the form of topological equivalence establish by homeomorphism). There is only a quantitative difference related with the presence of the non-isolated freeze singularity. The scale of appearance of this type singularity can be also estimated and in terms of 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 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 [68].