Do sewn singularities falsify the Palatini cosmology?

We investigate further (cf. arXiv:1512.01199, JCAP01 (2016) 040) Starobinsky cosmological model $R+\gamma R^2$ in the Palatini formalism with Chaplygin gas and baryonic matter as a source. For this aim we use dynamical system theory. The dynamics is reduced to the 2D sewn dynamical system of a Newtonian type (a piecewise-smooth dynamical system). We classify all evolutional paths in the model as well as trajectories in the phase space. We demonstrate the presence of a degenerate freeze singularity (glued freeze type singularities) for the positive $\gamma$. In this case it is a generic feature of early evolution of the universe. We point out that a degenerate type III of singularity can be considered as an endogenous model of inflation between the matter dominating epoch and the dark energy phase. We also investigate cosmological models with negative $\gamma$. It is demonstrated that $\gamma$ equal zero is a bifurcation parameter and dynamics qualitatively changes in comparison to positive $\gamma$. Instead of the big bang the sudden bounce singularity of a finite scale factor appears and there is a generic class of bouncing solutions sewn along the line $a=a_{\text{sing}}$. And we argue that the presence of sudden singularities in an evolutional scenario of the Universe falsifies the negative $\gamma$ in the Palatini cosmology. Only very small values of $\Omega_{\gamma}$ parameter are admissible if we requires that agreements physics with the $\Lambda$CDM model. From the statistical analysis of astronomical observations, we deduce that the case of negative values of $\Omega_\gamma$ can be rejected even if it may fit better to the data.

R + γR 2 in the Palatini formalism with Chaplygin gas and baryonic matter as a source. For this aim we use dynamical system theory. The dynamics is reduced to the 2D sewn dynamical system of a Newtonian type (a piecewise-smooth dynamical system). We classify all evolutional paths in the model as well as trajectories in the phase space. We demonstrate the presence of a degenerate freeze singularity (glued freeze type singularities) for the positive γ. In this case it is a generic feature of early evolution of the universe. We point out that a degenerate type III of singularity can be considered as an endogenous model of inflation between the matter dominating epoch and the dark energy phase. We also investigate cosmological models with negative γ. It is demonstrated that γ equal zero is a bifurcation parameter and dynamics qualitatively changes in comparison to positive γ. Instead of the big bang the sudden bounce singularity of a finite scale factor appears and there is a generic class of bouncing solutions sewn along the line a = a sing . And we argue that the presence of sudden singularities in an evolutional scenario of the Universe falsifies the negative γ in the Palatini cosmology. Only very small values of Ω γ parameter are admissible if we requires that agreements physics with the ΛCDM model. From the statistical analysis of astronomical observations, we deduce that the case of negative values of Ω γ can be rejected even if it may fit better to the data.

Introduction: Cosmology with Chaplygin gas in Palatini formalism
Today's modern cosmology suffers problems which the standard theory, that is the ΛCDM model derived from Einstein's general relativity, is not able to explain. There are such issues like dark matter and dark energy, origin and source of inflation or large scale structure which are widely investigated from many different points of view. Satisfactory explanations have not been set so far which makes searching and examination of new models still desirable. Due to that fact we have proposed [1] and examined a model which modifies the standard one in two different ways. The first modification refers to the gravitational Lagrangian which we enriched with the so-called Starobinsky term γR 2 keeping in mind that the parameter γ is small and being about to estimate by observational data and analytical analysis. We have used the Palatini formalism [2][3][4], which treats a metric and a connection as independent dynamical variables. The connection is used to construct Riemann and Ricci tensors, while contraction with the metric provides (generalized) Ricci scalar (see e.g. [3,5]).
Palatini approach to the description of gravitational field was originally introduced by Einstein himself [6, p. 415] but historical misunderstanding decides its name in this context [6, p. 191, 485]. The main idea of this formalism is to treated the connection Γ appeared in the definition of the Ricci tensor as a variable independent of the spacetime metric g. Therefore there is no special reasons to apply the Palatini variational principle in GR if we have metric formalism. The situation changes if we considered Extended Theory of Gravity (ETG) because in these cases both metric and Palatini variational principle satisfy different field equations which in order can give rise to different physics [6, p. 486, 769]. Application of the Palatini approach in the context of cosmological investigations has been the subject of many papers [6, p. 211, 721, 722, 843, 1129]. The Newtonian potential can be obtained if we considered the weak-field limit of ETG and its relations with a conformal factor [6, p. 794]. In particular, it has been shown [7,8] that vacuum solutions of Palatini gravity differ from GR by the presence of cosmological constant. Due to this fact the values of cosmological constant which are admitted by solar system tests are many order of magnitude bigger than the values obtained from cosmological estimations [9].
Palatini approach is very important in cosmology as the equations of motion are the second order differential equations in comparison to the standard metric approach in which the field equations give rise to the fourth order ones (see e.g. [1] and references therein). The second modification was taken to the matter part of the modified Einstein equation. Instead of considering perfect fluid with barotropic equation of state p = ωρ we have studied the so-called generalized Chaplygin gas which has also, similarly to cosmological constant, a negative pressure. It has gained a lot of attention in cosmology recently [10][11][12][13][14][15][16][17][18][19][20] as it combines dark energy (cosmological constant) and dark matter into one component. Moreover, this is the only fluid known up to now which has a supersymmetric generalization [21,22]. It seems to be a very important and interesting object to study. We would like to mention that in [1] we have obtained a very good agreement of our model with observational data. We were also able to find an upper bound on the value of the parameter γ in order to locate a singularity, which has appeared in the model, before the recombination epoch. It turn out that the singularity is of the type III and provides an intermediate inflation phase to the evolution of our Universe. That is, the model provides four cosmic evolution phases while the Big Bang singularity is preserved: the decelerating phase dominated by matter, an intermediate inflation phase corresponding to the type III singularity, a phase of matter domination (decelerating phase) and finally, the phase of acceleration of the current universe. Now on, for the reader convenience, let us shortly remind some properties of the approach that we are using.
The general action of the Palatini f (R)-gravity is written in the standard way where f (R) is function of the Ricci scalarR = g µνR µν (Γ). One should notice that the Palatini scalarR is constructed with both objects, that is, the metric and connection. The action S m is a matter action which is independent of the connection but includes other scalar fields and depends on g µν . Because the Lagrangian for matter does not depend on the connection [1] in the our model with barotropic matter satisfying equation of state p = p(ρ) can be postulated like in GR in the following form [23] Variation of the total action with respect to the metric gives rise to the following field equations where prime denotes differentiation with respect toR while the energy-momentum tensor T µν is obtained by the variation of the matter action with respect to g µν . We have also used the geometric units 8πG = c = 1. After taking g -trace of (1.3) one obtains the structural equation in the form where T is trace of the energy-momentum tensor. We will discuss that equations a bit later. The variation with respect to the independent connectionΓ provideŝ which is the indication that the connection is the Levi-Civita connection of the metric h µν = f (R)g µν . That is, the metric h µν is a conformally related metric to the physical metric g µν . It gives rise to the conclusion that the conformal factor f (R), later labeled by b, must be a non-negative function. One should also notice that in the scalar-tensor representation of the f (R) gravity in the Palatini approach the scalar field φ is represented by f (R). Moreover, converting the action to the Einstein frame one gets the negative coupling to the matter part for the case of f (R) < 0. Equation (1.4) for different choices of function f (R) leads to the algebraic equation depending onR. Such an algebraic equation may happen to be difficult to solve. Postulating the dependenceR(T ), the structural equation becomes a linear differential equation. If we putR(T ) = −T in an analogous way to general relativity then we obtain differential equation in the form One gets immediately a simple solution and therefore this class of functions f (R) defines a full range of choices basing on an analogy to general relativity. The energy-momentum tensor T µν satisfies the metric covariant conservation law ∇ µ T µν = 0 since we are considering the Palatini f (R) gravity as a metric theory [5]. Hence, the continuity equation is givenρ + 3H(ρ + p) = 0, (1.8) where˙≡ d dt denotes the differentiation with respect the cosmological time, ρ and p are energy density and pressure, respectively, of the energy momentum for the perfect fluid. The variable H = d(log a) dt denotes as usually the Hubble parameter. For completeness, the form of the equation of state p = p(ρ) should be postulated in order to obtain the scale factor a(t) dependence for the pressure and energy density. Due to that fact, we would like to consider a very special case which is perfect fluid for the generalized Chaplygin gas as a source of gravity [11]. The equations of state is taken as where A is positive constant and 0 ≤ α ≤ 1. Note that negative pressure of generalized Chaplygin gas, for A > 0, suggests naively that as ρ goes to zero, p diverges to minus infinity. But this reasoning is not true. The pressure as well as energy density satisfy the continuity condition which gives rise to the relation (1.10) We parametrize ρ(a) dependence through the physical parameters ρ ch,0 and assume 0 ≤ , today c 2 s = αA s < 1 as the consequence of 0 < α ≤ 1. Note that ρ(a) has a lower limit, namely ρ(a) ≥ ρ ch,0 (A s ) 1 1+α . Note that the Chaplygin gas does not violate null energy condition ρ + p ≥ 0 because ch,0 (1 − A s ) > 0 in the case considered. Therefore, for small values of the scale factor, the density ρ(a) behaves like a dust matter ρ(a) ≈ a −3 . Instead, for the large scale factor one gets ρ(a) = A 1/(1+α) , i.e. the effect of the cosmological constant. We use the idea of the Chaplygin gas [25] because such an equation of state interpolates between the matter dominating phase and the quintessence epoch at which the Λ is dominating. There is also an intermediate phase mimicking the Zeldovich stiff matter domination. Moreover, for α = 0 Chaplygin gas corresponds exactly to the presence of the cosmological constant (dark energy) and dust (dark and baryonic matter) 2 . We recall that in the modified gravity framework 'fluid dark energy' can be replaced by the cosmological constant ensuing from the modification of the gravitational action. Therefore, in the model under consideration, the null energy condition is not violated and the bounce is a consequence of the modification of the Einstein equations, i.e. the presence of the additional term γR 2 in the Lagrangian, when γ is negative.
Since the exact form of a function f (R) is not known, one needs to consider some effective theories probing theoretical possibilities of this approach to gravity. Usually some a priori truncated polynomial form with respect to the scalar field and their inverse are proposed. It enables us to study how different problems of contemporary cosmology like dark energy and dark matter issues can be solved [1,5,11,[26][27][28][29]. Now let us choose the simplest modification of the general relativity Lagrangian already mentioned as a simple solution of the structural equation, that is, let us consider f (R) =R + γR 2 (1.11) for which one deals with the relation Finally, after substitution formulas (1.9) and (1.10) we obtain the following a(t)-dependence for the Palatini scalarR where and ρ ch,0 is the present value of ρ ch , k = −1, 0, +1 is the space curvature and A s = A ρ 1+α ch,0 [1,26,27]. This model was examined by us [1] where we demonstrated how the Palatini formulation modifies evolutionary scenarios with respect to the ΛCDM model so it can be considered as a natural extension of the standard cosmological model. In this paper we focus our attention on a study of the dynamics of the considered model methods provided by dynamical systems. We show that the dynamics of the model can be considered as a two-dimensional dynamical system of a Newtonian type ( [27,28]). It turns out that the phase space structure is more complicated that for the standard dynamical system because of the presence of the degenerate singularity which belongs to type III [30,31]. This singularity has an intermediate character [32] and shared evolutionary paths on two ΛCDM types of evolution (two-phases model with matter and dark energy domination epochs). The mathematical model of such dynamics is formulated with the help of notion of 'sewn dynamical systems' [33]. Following this approach the full trajectories are sewn trajectories of two cuts of half-trajectories along the singularity. We have found that a weak degenerate singularity of type III appears. This is a generic feature of the dynamics in the early universe.
Barrow and Graham have introduced the concept of singular inflation [34]. In this context the presence of a freeze singularity in the early evolution of the universe opens a discussion of modeling the inflation through a singularity of type III.
Let us summarize the work that we are going to present. The main goal of the paper is the investigation of the dynamics of homogeneous and isotropic cosmological model with the LagrangianR + γR 2 in the Palatini formalism. It will be demonstrated that the dynamics can be in general reduced to the form of a two-dimensional dynamical system of a Newtonian type. This enables us to investigate the dynamics in details in the configuration space as well as the phase space where the phase portrait revealing the global dynamics can be constructed. Due to this representation of dynamics it is possible to study all evolutionary paths for all admissible initial conditions.

Dynamical system approach in study of evolution of the Universe
Our case belongs to a class of cosmological models of modified gravity whose dynamics can be reduced to the form of a two-dimensional dynamical system in a Newtonian form [27,28]. Lagrangians have the form of the one for natural mechanical systems. Therefore, the Hamiltonian H has a kinetic term quadratic in momenta and the potential as a function of state variables. The motion of the system is along the energy levels H = E = const. Due to this reduction the Universe dynamics can be treated as a particle of unit mass moving in the one-dimensional potential. This enables us to classify all evolutionary paths in the configuration space and in the phase one as well. Dynamical systems of a Newtonian type are special because they describe the evolution of conservative systems like in classical mechanics. For such systems the construction of a phase portrait can be obtained directly from a functional form of the potential.
In cosmological implications a positional variable constitutes the scale factor a(t) while the localization of the critical points as well as their type are determined by a shape of the potential V (x). Let us remind some commonly used terminology and properties: 1. A static universe is represented by a critical point of the systemẋ = y,ẏ = − ∂V ∂x and it always lies on the x-axis,, that is, y = y 0 = 0, x = x 0 .
2. We will say that the point (x 0 , 0) is a critical point of a Newtonian system if that is a critical point of the function of the potential V (x), it means: is total energy of the system. Spatially flat models will admit the case y =ẋ; E = 0 while the ones with the spatial curvature k = 0 (constant) have E = − k 2 .
3. A critical point (x 0 , 0) belongs to a saddle type if it is a strict local maximum of the potential V (x).
4. If (x 0 , 0) is a strict local minimum of the analytic function V (x) then one deals with a center.
From the above it should be clear that the shape of potential function determines the critical points and its stability. The integral of energy levels defines the algebraic curves in the phase space (x, y) which are mimicking the evolution of the system with time. The eigenvalues of the linearized matrix satisfy the characteristic equation of the form λ 2 + ∂ 2 ∂x 2 V | x=x 0 = 0. They may be real or imaginary with vanishing real parts (non-hyperbolic critical points) and hence there are another critical points beyond the centers and saddles. The centers are structurally unstable [35] in opposite to saddles which represent structurally stable critical points.
3 Classification of the trajectories representing evolution of the model

Classification of the possible evolutional scenarios in the configurational space
In order to classify all evolutionary paths in the configuration space we treat the cosmic evolution as a simple mechanical system with the natural form of Lagrangian L =ȧ where from now on we will use dot as˙≡ d dτ and τ is rescaled cosmological time, that is, As it was already mentioned, the equation (3.1) has a simple mechanical interpretation of the evolutionary of the universe in term of positional variable a(t). Hence the model dynamics has dynamics of particle motion of unit mass in the potentialṼ over the energy level.
Such an interpretation enable us to examine admissible trajectories and their classification in configuration and phase space. All information about dynamics are coded in the geometry of the potential function.
The diagram of the potential function (3.2) for typical values of model parameters is presented in fig. 1. The shaded region represents a non-physical domain forbidden for motion of a classical system for whichȧ 2 ≥ 0. We also define function V (a) in the form The motion of system takes place in the configuration space {a : a ≥ 0} over the 'energy' level E = 0, i.e., the Hamiltonian is of the form H(p, a) = 1 2 p 2 a + V (a) = 0 ≡ E = 0. Different energy levels determine corresponding types of evolution ( fig. 2).

The boundary of the domain admissible for motion is
Of course it is set with the boundary Note that domain E − V < 0 beyond this boundary is forbidden for classical motion. Let us classify all evolutionary scenarios in the configuration space: 1. O 1 -oscillating universes with initial singularities; 2. O 2 -'oscillatory solutions' without the initial and final singularity but with the freeze singularity; 3. B -bouncing solutions; 4. E 1 , E 2 -solutions representing the static Einstein universe; 5. A 1 -the Einstein-de Sitter universe starting from the initial singularity and approaching asymptotically static Einstein universe; 6. A 2 -a universe starting asymptotically from the Einstein universe, next it undergoes the freeze singularity and approaches to a maximum size. After approaching this state it collapses to the Einstein solution E 1 through the freeze singularity; 7. A 1 -expanding universe from the initial singularity toward to the Einstein universe E 2 with an intermediate state of the freeze singularity; 8. EM -an expanding and emerging universe from a static E 2 solution (Lemaitre-Eddington type of solution); 9. I -an inflectional model (the relation a(t) possesses an inflection point), an expanding universe from the initial singularity undergoing the freeze type of singularity.
The last two solutions EM and I lie above the maximum of the potentialṼ . Moreover, one deals with a singularity at which the accelerationä is undefined because left-hand side limit of the derivative of the potential is positive while the right-hand side limit is negative. This kind of a singularity should be treated as sewing two singularities: one type III singularity and one type III singularity with reverse time. In other words, this special type of singularities goes beyond the classification of four types of finite-time singularities [36,37].

Phase portrait from the potential
Due to Hamiltonian formulation of dynamics it is possible also to classify all evolution paths in the phase space by constructing the phase portrait of the system It would be useful to look on the dynamical problem from the point of view of sewn dynamical systems [38,39]. Our strategy is following. For construction of a global phase portrait one divides dynamics in two parts, that is, one for the scale factor a < a fsing and other for a > a fsing . Such a construction divides the configuration space which is glued along the singularity.
The re-parametrization τ → σ can be performed and corresponding dynamical systems assumes the following form for a < a fsinġ a = x, (3.8) where V 1 = V (−η(a − a s ) + 1) with respect to the new time σ and η(a) notes the Heaviside function.
In an analogous way for the domain configuration space {a : a > a fsing } we havė where η is the Heaviside function. The phase portrait is presented in fig. 3 and belongs to the class of sewn dynamical systems. Recently such systems have been applied to the modeling of cosmological evolution and inflation [30,31,40]. The classification of all possible evolutionary paths in the phase space completed our previously classification in the configuration space.
Note that trajectory represented the evolution of our Universe is located in a close neighborhood of a trajectory of the spatially flat model, i.e. a universe is expanding and starting from the initial singularity, going through a freeze singularity and after accelerating phase is going toward the de Sitter attractor.

Classification of the trajectories for Ω γ < 0
For completeness let us consider the case of the spatially flat model with Ω γ < 0. Because the value of γ = 0 is a bifurcation parameter (the dynamics qualitatively changes under the change of sign γ) we consider separately both these cases.
In this case one can apply a simple method of the classification of the evolutionary paths based on the consideration of a boundary of the domain admissible for a motion in the configuration space H 2 ≡ 0.
From (1.14) the equation of the boundary curve assumes the following form where Ω ch = Ω ch,0 3As K  Finally or 14) The diagram of Ω γ (a) function derived from (3.12) is shown in fig. 4. The functions (3.13) and (3.14) are always negative and their approximations for large and small value of the scale factor are, respectively: From asymptotes for the small scale factor the effects of curvature and non-zero Ω γ are negligible. The dynamics of the model is equivalent to the ΛCDM dynamics. This means that matter effect dominates. From asymptotes for the large scale factor we obtain the effects of both matter and curvature are negligible. The dynamics corresponds to the accelerating phase caused by dark energy. Now on, let us consider levels of Ω γ = const ≤ 0. In consequence we obtain two types of possible evolutionary paths O -oscillating models without the initial and final singularities; B -models with a bounce instead of an initial singularity like in the case of Ω γ > 0. Note that class of bouncing models as well as oscillating ones is generic. In fig. 5 we plot the diagram of the H 2 (a) relation. It demonstrates that the admissible domain for the motion of classical trajectories is the union of two separated and disjoint domains occupied by the trajectories. In the domain situated on the left in which H 2 > 0 the motion is bounded by a Big-Bang singularity: a = 0, H = ∞, and a maximum value of the scale factor on the Ω k,0 .
In the domain {a : a > a min } of the configuration space there is the bouncing type of the Universe evolution. Of course the value a min depends on the initial conditions, i.e., Ω k,0 . This relation is illustrated in the diagram of function V (a) shown in fig. 6. Taking different energy levels depending on the value of Ω k,0 one gets different evolutionary scenarios. In the domain of the configuration space {a : a < a max } we have an oscillating solution with the big-bang which after reaching the maximum size will re-collapse to the second singularity.
In turn in the domain {a : a > a sing : b(a sing ) = 0} there are two types of trajectories: one representing evolution of the oscillating closed models without initial and final singularities while the second one for flat (Ω k,0 = 0), closed (Ω k,0 > 0) and open (Ω k,0 < 0) models with the bounce.
There is a vertical line, which separates two disjoint domains of the configuration space. The equation a = a sing can be simply obtained from the definition of b function, namely In the special case of the Chaplygin gas (α = 1) we obtain −2Ω γ Ω ch,0 (3A s ) and hence one deals with the algebraic equation of second order where β 2 = 12Ω 2 γ Ω 2 ch,0 A s Because K ∈ [0, 3), the above equation has one solution which is or, in the terms of the scale factor a sing .
(3.23) Equation (3.21) has a real solution for K ∈ [0, 3) if the parameter β 2 > 4. On the line {b = 0} there is a singularity point (see fig. 7), that is, a = a sing , a sing = ∞. If we go back to the original time thenȧ = 0. It is a singularity of type II called a sudden singularity at which a, ρ eff , and Hubble parameter remain finite (Ḣ diverges). They are past singularities like a big-demarrage arising in models with the generalized Chaplygin gas. Figure 7 is the phase portrait of the original system under re-parametrization of time t → σ : dσ = |b| |b+d/2| dt. The function σ = σ(t) is drawn in the Fig. 8. For this case the dynamical system is expressed by Of course this re-parametrization is singular on the line {b = 0}. One should notice that there is no inverse transformation t = t(σ). Figure 9 illustrates that the function b(a) changes the sign if it passes through the zero on the a-axis. On the other hand the function (b + d/2) is not singular (see fig. 10). Since the re-parametrization is a non-smooth function when b near the singularity of the finite scale factor. The typical trajectory starts this critical point (representing a sudden singularity) and evolves toward the Sitter universe where effects of the curvature are also negligible. On the phase portrait there is also the critical point of the saddle type. This critical point corresponds to the maximum of the potential V = V (a). That is, it is a decreasing function of the argument so the universe is decelerating one. Therefore the trajectories on the right-hand side of the saddle point represents accelerating models. Finally, from the comparison of ΛCDM dynamics with our model it is seen that the initial singularity is replaced by the sudden singularity in the past like the demarrage singu- larities appeared in the cosmological models with the Chaplygin gas.
In fig. 11 it is shown the diagram of the scale factor for the flat models which the generic description of the evolution near sudden singularity. On the phase portrait this singularity lies on the circle at infinity a 2 + a 2 = ∞, where a = a sud.sing. . If we return to the original cosmological time t thenȧ at this critical point at the infinity is correspondingȧ = 0. In the consequence the Hubble parameter is finite. For the FLRW model with initial singularity this follows divergence of conformal time

Degenerate singularity of type III as a model of endogenous intermediate inflation
From the formula (1.14) it is possible to detected singularity of type III, called also a freeze singularity. Our analysis shows that this type singularity is a generic property of the early evolution of the universe. If we consider singularities in FLRW models, which is filled of perfect fluid with effective energy density ρ eff and pressure p eff then all singularities can be classified on the four groups [37]. The first class is a Big Rip (type I) singularity, where the energy density, pressure, and the scale factor diverge. Sudden singularity (type II) happens when the scale factor and effective energy density are finite values but pressure diverges. Big Freeze singularity (type III) is observed when effective energy density and pressure diverge at a finite value of the scale factor while Big Brake (type IV) for the finite scale factor, effective energy density, and pressure with a divergence in the time derivative of the pressure or change of energy density rate. In that context our singularities belongs to the type III.
In our approach effective energy density and pressure (as well as the coefficient of the equation of state w eff ) can be simply expressed in terms of the potential, namely In our case potential as well as the effective pressure diverges. This singularity is the singularity of acceleration as the derivative of the potential goes to plus infinity on the left of this point while on the right side one has minus infinity. At this singularity da/dt also diverge while the scale factor is finite. In the diagram of the scale factor as a function of time one can observe how the function a(t) changes an inflection along the vertical line t = t fsing . Such types of singularities appear in the context of Loop Quantum Cosmology [37,[41][42][43] where are called a hyper-inflation state [44]. There exists a close relation between Palatini gravity and an effective action of Loop Quantum Gravity that reproduces the dynamics [45,46] of considered f (R) models of polytropic spheres in Palatini formalism f (R) = R±λR 2 (λ has order of squared Planck length). The freeze type of singularity in the model under consideration is a solution of the algebraic equation where K ∈ [0, 3). If α = 0 then equation (4.5) simplifies to The solution of the above equation is (4.7) From equation (1.15) and (4.7) one finds an expression for a value of the scale factor for the degenerated freeze singularity 8) or in the term of redshift This singularity, which is the horizontal inflection singularity point a = a fsing , corresponds to the diagram of the scale factor of t ( fig. 13). Note that at this singularity the potentialṼ diverges.
For the special case of the Chaplygin gas (α = 0) the exact formulas or value of a fsing can be obtained [1]. The diagram of function b + d 2 = f (A s , Ω ch,0 , Ω γ , α) is presented in fig. 14. From our numerical analysis (see fig. 15) we have obtained that a single freeze type singularity is generic property of dynamics for the broad range of model parameters (A s , α, Ω ch,0 , Ω γ ).

Singularities and astronomical observations
In this section, we discuss the status of singularities appearing in the model under the consideration. If we compare this model with the ΛCDM model which formally can be obtained after putting α = 0 then one can conclude that the latter inherits these singularities. Basing on the estimation of the model parameters performed in our previous paper [1] one can calculate the value of the scale factor (or redshift) corresponding to this event in the history of the universe when singularities appear. This value depends on the value of density parameters Ω γ , A s , α, H 0 . From numerical simulations we obtain that this value is only sensitive on Ω γ parameter and the dependence on the parameter α is very weak. This means that corresponding values of redshifts obtained for these singularities in the case α ∈ (0, 1) do not differ. In figures 12, 15 we illustrate how values of redshifts (the scale factor marked in the figures) depend on the crucial value of density parameter Ω γ . While for positive Ω γ this function is a growing function of the scale factor, but for negative values of Ω γ it is a decreasing function of the scale factor. In both cases it is a monotonic relation. Of course for the case of Ω γ = 0 the correspondence with the ΛCDM model is achieved and both sudden and freeze singularities vanish. Therefore the presence of singularities of the model under consideration should be treated as a property which is strictly related with the Palatini formalism.
For the best fitted values of the model parameters obtained in our previous paper [1] one could estimate simply the value of redshift corresponding of the sudden singularity at the bounce. This value is 1103.67. In the case of the positive Ω γ , this parameter lies on boundary which unable us formulation of an analogous conclusion.
In the monograph by Capozziello and Faraoni [6, p. 73], section 3.4.2, some problems with the Palatini formalism were addressed. The authors remarked that the f (R) gravity suffers from two serious problems. The first problem is the presence of curvature singularities at the surface of stars and the second one is incompatibility with the Standard Model of particle physics. Let us consider our model in the context of singularities. Our remark in the context of cosmology is that the Palatini formalism rather generates singularities than suffers from the singularities like in the case of stars.
Capozziello and Faraoni also noted that different formulation of gravity can give rise to new physical effects. If we assume that freeze singularities can appear before the recombination epoch and sudden singularities being before the nucleosynthesis then such requirements should be treated as minimal conditions which guarantee that physics does not change in comparison to the ΛCDM model. In consequence if we assume that the nucleosynthesis epoch was for redshift z = 3 × 10 8 (see e.g. [47]) then the value of Ω γ parameter should belong to the interval (−10 −25 , 0) for the case the negative Ω γ parameter. If we analyze the likelihood functions for values of model parameters [1] for the negative Ω γ parameter then we get that the probability of appearing the sudden singularity after the nucleosynthesis epoch is 1 − 10 −16 while the probability of appearing one before the nucleosynthesis epoch is 10 −16 . In consequence one can reject the case with the negative values of Ω γ even it is favored by the data. If the recombination epoch takes place after the freeze singularity then it is required that the values of the positive Ω γ parameter belong to the interval (0, 10 −9 ).

Conclusions
We have classified all evolutionary paths of the cosmological model f (R) =R + γR 2 the Palatini formulation due to reduction of model dynamics to the dynamical system of a Newtonian type. We have found a phase space structure of dynamics organized through the two saddle points representing the static Einstein universe and center.
Moreover, the localization of the freeze singularity during the cosmic evolution of the early universe was presented. We investigated in details how this localization depends on density parameters for contribution originating from the presence of R 2 term in the Lagrangian. The model possesses four phases: a first matter dominating deceleration phase, an intermediate phase of inflation, a second matter dominating deceleration phase, and a late accelerating phase of evolution of the current universe. This evolutionary scenarios becomes in agreement with the ΛCDM model only if the freeze singularity is shifted to the matter dominating initial singularity.
While in our case acceleration (and in consequence pressure) is undefined (note that all four types assume that acceleration is a well defined function for which we calculate the limits). It is similar to a singularity of type III (finite scale factor singularity) at which the scale factor remains finite, but both ρ and p diverge (as well as the Hubble parameter H). A particular example of such a singularity is a 'big-freeze' singularity (both in the past and the future). They are characterized by the generalized Chaplygin equation of state [42]. In the model under consideration one deals with the situation in which freeze singularities in the past and in the future are glued. This type of degeneration is beyond the standard classification. We will call this type of weak singularity a degenerate freeze singularity. Note that trajectories in the phase space pass through this singularity and continue their evolutions.
We would like to draw the reader's attention to the conclusion that the dynamics presented in terms of dynamical systems of a Newtonian type enable us to classify all evolutionary paths. That it, the dynamics is reduced to the 2-dimensional sewn dynamical system. We have constructed the phase portrait which represents the global dynamics. From the construction, one gets all evolutionary paths admissible for all initial conditions. The set of sewn trajectories is a critical point located at the infinity (ȧ = ∞, a = a fsing ). The trajectories pass through this critical point. On the phase portrait the trajectory of the flat model (k = 0) divides trajectories in the phase space on the domains occupied by closed (k = +1) and open (k = −1) models. From the phase portrait one can derive conclusion that while all open models possess freeze singularity there are also presented bouncing models and closed ones with initial and final singularities without the freeze singularity. The latter class of models without freeze singularity is also generic. Note that there exists also a class of generic models without the initial singularity.
We have also classified all solutions for the negative Ω γ . This case is favored by the astronomical data [1]. The classification is performed in both configurational and phase spaces. From the classification one can conclude that that there is a generic class of cosmological models in which the Big Bang singularity is replaced by the bounce.
For the case of negative γ the phase space is in the form of two disjoint domains without any causal communication between them. The simple method of the avoidance the domain on the left-side of line {b = 0} is to assume that f (R) > 0. For our model this means that b > 0 at the very beginning. The domain of the phase space for a negative value of b can be interpreted as a region occupied by trajectories with a negative constant coupling of matter and scalar field.
The final conclusion is that the extended cosmology with γR 2 term in the Palatini formalism has status sewn dynamical systems from the point of view dynamical systems theory. A point of sewing represents the freeze singularity for γ > 0 or sudden singularity in the opposite case. From the cosmological point of view a sudden singularity is represented by a bounce while a freeze one is represented by an inflationary phase.
The type of singularities in the Starobinsky model in the Palatini formalism crucially depends on the sign of the parameter γ. If γ is positive then we obtain new singularity apart from the Big-Bang singularity. This type singularity is similar formally to the singularity of type III but has a more complex nature, because is composed de facto with two finite scale factor singularities of type III (in a past and a future). The behavior of the scale factor near this type singularity can be obtained from the expansion of function t = t(a) near an inflection point.
Therefore, a − a sing ∝ −(t sing − t) 1/2 for t → t − sing +(t − t sing ) 1/2 for t → t + sing . (6. 2) The acceleration p =ä goes to +∞ as t → t − sing and −∞ as t → t + sing . In the consequence p is undefinite because of the lack of continuity in this point of gluing. Note that the system under consideration is an example of a piecewise-smooth dynamical systems [48].
Our general conclusion is that the presence of a sudden singularity in the model falsifies the the negative case of γ in the Palatini cosmology. The agreements of the physics which imply evolutional scenario of the universe with the observations requires that the value of Ω γ parameter should be extremely small beyond the possibilities of contemporary observational cosmology. From the statistical analysis of astronomical observations, we deduce that the case of negative values of Ω γ can be rejected.
There are in principle two interpretation of obtained results. In the first interpretation, we treated singularities as artifacts of the Palatini variational principle which in consequence limits its application. However there is also another interpretation. We in consideration of MGT take the simplest, quadratic correction of general relativity. Of course we do not know the exact form of f (R) and our theory plays the role of some kind effective theory. It is possible that if additional terms in a Taylor expansion be included these singularities disappear in a natural way. It seems to be interesting to build cosmology in the Palatini formalism under higher order terms (with respect to the Ricci scalar R or its inverse R −1 ) in a Taylor expansion for checking whether above mentioned singularities can appear during the cosmic evolution.