Do sewn up singularities falsify the Palatini cosmology?

We investigate further (cf. Borowiec et al. JCAP 1601(01):040, 2016) the Starobinsky cosmological model R+γR2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R+\gamma R^2$$\end{document} in the Palatini formalism with a Chaplygin gas and baryonic matter as a source in the context of singularities. The dynamics reduces to the 2D sewn dynamical system of a Newtonian type (a piece-wise-smooth dynamical system). We demonstrate that the presence of a sewn up freeze singularity (glued freeze type singularities) for the positive γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} is, in this case, a generic feature of the early evolution of the universe. It is demonstrated that γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} equal zero is a bifurcation parameter and the dynamics qualitatively changes as the γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} sign is changing. On the other side for the case of negative γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} instead of the big bang the sudden bounce singularity of a finite scale factor does appear and there is a generic class of bouncing solutions. While the Ωγ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{\gamma } > 0$$\end{document} is favored by data only very small values of Ωγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{\gamma }$$\end{document} parameter are allowed if we require agreement with the Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}CDM model. From the statistical analysis of astronomical observations, we deduce that the case of only very small negative values of Ωγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _\gamma $$\end{document} cannot be rejected. Therefore, observation data favor the universe without the ghost states (f′(R^)>0\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})>0$$\end{document}) and tachyons (f′′(R^)>0\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})>0$$\end{document}).

general relativity (GR), is not able to explain. Issues like the nature of dark matter and dark energy, the origin and source of inflation or large scale structure 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 have enriched with the so-called Starobinsky term γR 2 [2] keeping in mind that the parameter γ is small and estimated by observational data and analytical analysis. We have used the Palatini formalism [3][4][5], which treats the metric and the symmetric connection as independent dynamical variables. The Palatini formalism modifies the Starobinsky model by providing the other field equations and degrees of freedom in comparison to the original purely metric Starobinsky model. The torsion free connection is used to construct Riemann and Ricci tensors, while the contraction with the metric provides a (generalized) Ricci scalar (see e.g. [4,6]).
Let us notice that the Palatini approach we are using for the description of the gravitational field was originally introduced by Einstein himself (see e.g. [7, p. 415] and references therein) but historical misunderstanding decided on its name in this context [7, p. 191, 485]. The main idea of this formalism is to treat the symmetric (i.e. torsion free) connection in the definition of the Ricci tensor as a variable independent of the spacetime metric g. That is, the connection is not necessarily a Levi-Civita connection of the metric g, since the dependence on g arises from solutions of field equations. E.g., there is no special reason to apply the Palatini variational principle in GR where the two formalisms are equivalent. The situation changes if one considers the extended theory of gravity (ETG) because in these cases both metric and Palatini variational principles satisfy different field equations, which lead to different physics [7, p. 486, 769]. The applications of the Palatini approach in the context of cosmological investi-gations have been the subject of many papers [7, p. 211, 721, 722, 843, 1129] and references therein. A Newtonian potential can be obtained if we considered the weak-field limit of ETG and its relations with a conformal factor [7, p. 794]. In particular, it has been shown [8,9] that vacuum solutions of Palatini gravity differ from GR by the presence of a cosmological constant. It appears that the values of the cosmological constant which are admitted by solar system tests are many orders of magnitude bigger than the values obtained from cosmological estimations [10]. The Palatini approach is very important in cosmology as the equations of motion are second order differential equations, to be compared to the standard metric approach in which the field equations give rise to fourth order ones (see e.g. [1] and references therein).
The second modification in our approach is performed to the matter part of the modified Einstein equations. Instead of considering a perfect fluid with the barotropic equation of state, p = ωρ, we study the so-called generalized Chaplygin gas, which has also, similarly to cosmological constant, a negative pressure. Chaplygin gas has recently gained a lot of attention in cosmology [11][12][13][14][15][16][17][18][19][20][21] as it merges 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 [22,23]. For this reason it seems to be an important and interesting object to study.
We would like to mention that in [1] we have obtained, under the same assumptions, a model which reaches very good agreement 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 turns out that this singularity is of type III and provides an intermediate ephemeral inflation era to the evolution of our universe. More exactly, the model possesses four phases of cosmic evolution while the big bang singularity is preserved: the decelerating phase dominated by matter, an intermediate ephemeral 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. The present paper extends and intensifies the analysis presented in [1]. Now on, for convenience of the reader, let us briefly recall the main properties of the formalism which is used in the paper.
The general action of the Palatini f (R)-gravity is written in the standard way, where f (R) is a function of the generalized Ricci scalar R = g μνR μν (ˆ ). One should notice that the Palatini scalar R is constructed with both objects, that is, the metric and the torsion free 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 model with a generalized barotropic matter satisfying equation of state p = p(ρ), it can be postulated similar to GR to have the following form [24]: The variation of the total action with respect to the metric gives rise to the following field equations: where a prime denotes differentiation with respect toR, while the energy-momentum tensor T μν is obtained by variation of the matter action with respect to g μν . We have also used the geometric units 8π G = c = 1. By taking the g-trace of (3) one obtains a structural equation in the form where T is the g-trace of the energy-momentum tensor. We shall discuss these equations later.
The variation with respect to the independent symmetric 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 metric conformally related 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 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 (4) for different choices of the function f (R) leads to the algebraic type equation depending onR. Such an 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 GR, then we obtain a 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 based on the analogy to GR. 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 [6]. Hence, the continuity equation is given by where˙≡ d dt denotes the differentiation with respect to the cosmological time, and ρ 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 usual, the Hubble parameter. For completeness, the form of the equation of state p = p(ρ) should be postulated in order to obtain the dependence on the scale factor a(t) 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 [12]. The equation of state is taken as where A is positive constant and 0 ≤ α ≤ 1. Note that a 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 We parametrize the ρ(a) dependence through the physical parameters ρ ch,0 and assume 0 ≤ A s ≤ 1 (and B > 0) 1 [11]. In this parametrization the square of the speed of sound is αρ 1+α ch,0 A s ρ −(α+1) , today c 2 s = α A s < 1 as a consequence of 0 < α ≤ 1. Note that ρ(a) has a lower limit, that is, ρ(a) ≥ ρ ch,0 (A s ) 1 1+α . Moreover, the Chaplygin gas does not violate the null energy condition ρ + p ≥ 0 because ch,0 (1 − A s ) > 0 in the case under consideration. Therefore, for small values of the scale factor, the density ρ(a) behaves like dust matter, ρ(a) ≈ a −3 , while for the large scale factor one gets ρ(a) = A 1/(1+α) , i.e. the effect of the cosmological constant. We use the idea of a Chaplygin gas [26] because such an equation of state interpolates between the matter dominating phase and the quintessence epoch at which dominates. There is also an intermediate phase mimicking the Zeldovich stiff matter domination. Moreover, for α = 0 the 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 1 In some models of fluid inflation negative values are also allowed [25]. 2 One should notice that the best fit obtained in [1] corresponds to a small value of α = 0.0194. 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 an exact form of the 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 fields and their inverses are proposed. It enables us to study how different problems of contemporary cosmology like the dark energy and dark matter issues can be solved [1,6,12,[27][28][29][30]. 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 (11) for which one deals with the relation obtained from (4), Finally, after substituting formulas (9) and (10) we obtain the following a(t)-dependence for the Palatini scalar: In this case the generalized Friedmann equation can be rewritten in the form H 2 = H 2 (a) as was done in [1]. Following that result, the relation H 2 where Here, the quantity ρ ch,0 denotes the present value of ρ ch , k = −1, 0, +1 is the space curvature and A s = A ρ 1+α ch,0 [1,27,28].
This model was examined by us [1] where we demonstrated how the Palatini formulation modifies evolutionary scenarios with respect to the CDM model. Henceforth, it can be considered as a natural extension of the standard cosmological model. In this paper we focus our attention on an investigation of the dynamics of this model provided by dynamical systems methods. We show that the dynamics of the model can be treated as a two-dimensional dynamical system of a Newtonian type [28,29]. It turns out that the phase space structure is more complicated than for the standard dynamical system because of the presence of the sewn up singularity which belongs to type III [31,32]. This singularity has an intermediate character [33] and divides evolutionary paths into 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' [34]. Following this approach the full trajectories are sewn up trajectories of two cuts of half-trajectories along the singularity. We have found that a weak sewn up 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 [35]. 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 detail 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.
2 Dynamical system approach in study of evolution of the universe As already mentioned above, 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 [28,29]. The Lagrangian has a form of one for standard mechanical systems. Therefore, the Hamiltonian H has a kinetic term quadratic in the momenta and the potential as a function of the 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 onedimensional potential. This enables us to classify all evolutionary paths in the configuration space and in the phase too. Dynamical systems of a Newtonian type are special because they describe the evolution of conservative systems as in classical mechanics. For such systems the construction of a phase portrait can be obtained directly from the functional form of the potential.
In cosmological applications the scale factor a ≥ 0 plays the role of positional variable x, while the localization of the critical points as well as their type are determined by a shape of the potential V (x). Let us recall 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 say that the point (x 0 , 0) is a critical point of a Newtonian system if it is a critical point of the function of the potential V (x), it means: is the total energy of the system. Spatially flat models admit the case y =ẋ; E = 0, while the ones with the spatial curvature k = 0 (constant) have From the above classification it should be clear that the shape of the potential function determines critical points and their stability. The integral of the energy levels defines the algebraic curves in the phase space (x, y) which imitate the evolution of the system with time. 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 other critical points apart from the mentioned ones. The centers are structurally unstable [36], opposite to saddles which represent structurally stable critical points.

Classification of the possible evolutional scenarios in the configuration 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 the Lagrangian L =ȧ where from now on we will use dot as˙≡ d dτ and τ is a rescaled cosmological time, that is, H 0 t = τ . The potential V a Fig. 1 The diagram presents the potentialṼ (a) for A s = 0.7264, α = 0.0194, and γ = 10 −9 . The shaded region represents a nonphysical domain forbidden for motion of a classical system for whicḣ As already mentioned, Eq. (16) has the simple mechanical interpretation of the evolution of the universe in terms of the positional type variable a(t). Hence the dynamics of the model is the same as the dynamics of a particle of unit mass moving in the potentialṼ along a constant energy level. Such an interpretation enables us to examine admissible trajectories and their classification in the configuration and phase space. All information as regards the dynamics is encoded in the geometry of the potential function.
The diagram of the potential function (17) for typical values of model parameters is presented in Fig. 1.
We also define function V (a) in the form The motion of the 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 with the boundary Note that the domain E − V < 0 beyond this boundary is forbidden for classical motion. Let us classify all possible 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; 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 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 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 are situated above the maximum of the potentialṼ . Moreover, one deals with a singularity at which the accelerationä is ill defined. The 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 together two singularities: a type III singularity with time running forward and a type III singularity with reverse time. In other words, this special type of singularity goes beyond the classification of four types of finite-time singularities [37][38][39].
In the paper we have considered the cosmological model with a new kind of singularities which is beyond the standard classification of cosmological singularities. The standard classification of singularities concerns only one-side future (big crunch type) or past (big bang type) singularities which are isolated. Notice that big bang and big crunch are mirror images of each other. In the model under consideration the compound singularity appears. It consists of two singularities sewn together: i.e. a future singularity on the evolution path for a < a sing sewn up with the past singularity on the evolution path in the domain a > a sing . The two components of this singularity are of the same type. In the domain a < a sing the universe expands toward a sudden future singularity at a = a min . Reversing the time in the domain a > a sing one sees that the universe collapses to a sudden past singularity and reaches a = a min . At this point the two singularities meet and the new singularity is created. It is the sewn up singularity of the future sudden singularity and past sudden singularity. We called it the sewn up sudden singularity. For the case of the sewn up freeze singularity the future and past freeze singularities meet at the point a = a sing . Summing up, we have found that the existence of sewn up (i.e. double-sided) nonisolated singularities is a generic feature of cosmological models in the Palatini formalism. Note that in such singularities the pressure is not well defined because its value on the left-hand side (future singularity) has different sign from the value taken at the right-hand side (past singularity).
Recently several dark energy models have been postulated as a hypothetical explanation of the conundrum of acceleration. The appearance of such a phase in evolution requires the violation of the strong energy condition, which is related with the presence of singularities in the future evolution of the universe [37].
Among them is the big freeze singularity [40][41][42], which takes place at a finite scale factor, in a (flat) Friedmann-Robertson-Walker (FRW) universe. For example the future big freeze singularity can appear in a FRW universe with a (phantom) generalized Chaplygin gas (GCG) [40].
The big freeze singularity (or type-III singularity in the nomenclature of [37]) takes place at a finite scale factor and a finite cosmic time in the (flat) FRW universe. At this singularity, both the Hubble rate and its cosmic time derivative blow up.
It is interesting that a big freeze singularity might be simply evidence of the change from a Lorentzian to a Euclidean signature in the brane cosmology [43].

Phase portrait from a potential
Due to the Hamiltonian formulation of the dynamics it is also possible to classify all evolution paths in the phase space by constructing the phase portrait of the system, on the phase plane (a, x) with the constraint a 2 = −2V (a), It would be useful to look at the dynamical problem from the point of view of sewn dynamical systems [44,45]. Our strategy is as follows. For the construction of a global phase portrait one divides the dynamics in two parts, that is, one for the scale factor a < a fsing and the 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 the corresponding dynamical system assumes the following form for a < a fsing : where is a function with respect to the new time, σ , where η(a) denotes the Heaviside function.
In an analogous way for the domain {a : a > a fsing } in configuration space we havė where V 2 = V η(a − a s ) and η 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 [31,32,46]. The classification of all possible evolutionary paths in the phase space completed our previously classification in the configuration space.
Note that the trajectory representing the evolution of our universe is located in a close neighborhood of the trajectory of the spatially flat model, i.e. the universe is expanding and starting from the initial singularity, going through a freeze singularity, and after an accelerating phase is going toward the de Sitter attractor.

Classification of 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 for the classification of the evolutionary paths, which takes into consider- where ch = ch,0 or  Fig. 4 Diagram of γ (a) dependence used for a classification of all evolutional paths of the spatially flat models with a negative γ . The type of evolution we obtain after consideration of levels = γ . We can simply discover oscillating models without the initial and final singularity, oscillating models with the initial and final singularity, models evolving to infinity with the initial singularity and bouncing models. From this analysis we conclude that for the physical case of sufficiently small values of γ all models possess a bounce instead of the initial singularity. The blue line represents b = 0 (Eq. (28)). Above the blue line is the region for b > 0 and this region corresponds to the physical domain. The region below the blue line is for b < 0 and represents the non-physical domain. The red line represents the function (29). Between the red line and the blue line H 2 < 0 (light gray domain). The gray domain also represents a non-physical region The diagram of the γ (a) function derived from (27) is shown in Fig. 4. The functions (28) and (29)  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 the 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 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 −1 . Because K ∈ [0, 3), the above equation has one solution which is or, in terms of the scale factor a sing ,

Equation (36) has a real solution for
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 the Hubble parameter remains 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 Fig. 8. For this case the dynamical system is expressed by 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 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 changes sign, the corresponding dynamical system possesses a discontinuity point on the line {b = 0}.
All properties of the model dynamics under consideration are summarized in the phase portrait, which is a picture of the global dynamics (see Fig. 7). The phase space is the union of the disjoint domains A = {a : a < a sud.sing } and B = {a : a > a sud.sing }. In both regions trajectories repre- That is, it is a decreasing function of its argument so the universe is a 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 singularities appearing in the cosmological models with the Chaplygin gas. Figure 11 shows the scale factor for the flat models as a function of time with 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 corresponds tȯ a = 0. In consequence, the Hubble parameter is finite. For the FRW model with initial singularity this is a divergence of conformal time dt a(t) because dt a(t) > 1 const dt t diverges as t → 0. The diagram of the relation between a sud.sing. ( γ ) is presented in Fig. 12.
The dynamical systems under consideration belongs to the class of dynamical systems for which the r.h.s. are in general piece-wise smooth functions of state variables (a, da/dτ ). In the phase space there are present vertical lines given by the equation a = a sing which on the r.h.s. are undefined at a finite domain. These lines are called lines of discontinuity or lines of singularity because the system breaks down in this domain. The lines of singularity divide the phase space into disjoint domains in which the trajectories are determined by smooth dynamical systems.  The lines of discontinuity correspond to the poles of the potential function. In any part of the diagram apart from the singularity points the potential is a smooth function and defines a part of smooth dynamical systems. In each smooth piece of the potential the dynamics is modeled by a smooth dynamical system, therefore the Cauchy problem is well defined. The global dynamics of the universe is obtained through the construction of sewing up solutions from the separate regions.
In the phase space the sewn up sudden singularity or the sewn up freeze singularity does appear. They are singularities of a finite scale factor type. Trajectories of the sudden-like singularity are sewn up by identification of opposite points at infinity (a = a sing , x = ∞) with (a = a sing , x = −∞).
For the case of the freeze like singularity the trajectories are connected from both separate regions at the point at infinity (a = a sing , x = ∞). This construction guarantees us the continuity of a solution which is passing through the singularity. The Cauchy problem is then well defined.
It was shown in [47] that, as in general relativity, the initial value problem is well formulated when the Palatini f (R) gravity is considered as a scalar-tensor theory with w = −3/2. The difference between the Cauchy problem in general relativity and Palatini f (R) theory is the appearance of the secondary constraint, which turns out to be the field equation for the scalar field f (R). It has an algebraic form relating the trace of the energy-momentum tensor and the scalar field. Moreover, the detailed discussion in [48] suggests that the Cauchy problem is also well posed but the analysis should be done case by case, since it cannot be performed for a general function f (R).
The dynamics of the model is governed by the first order differential equation called the Friedmann equation. Therefore the solution of the Cauchy problem requires specifying initial conditions. In the model that we are considering, the phase space is divided on two disjoint domains separated by the line a = a sing . Let us try initial conditions in the domain a < a sing . Then the solutions are determined in this domain breaking their evolution at a sing . Beyond this domain the Cauchy problem is unsolved. A similar construction can be performed in the domain a > a sing . If we choose the initial condition in this domain and inverse time, the Friedmann equation determines the evolution of the system, which breaks at the singularity point a sing .
This singularity point is of course a discontinuity point, but one can glue any trajectory at the left domain with the corresponding trajectory in the right domain at the singularity point using the symmetry of the a(t) function with respect to the singularity point. The behavior of the scale factor in the neighborhood of a singularity has the universal character because a = a sing is an inflectional point. We have a symmetry of a(t) with respect to a sing . Because of this symmetry any trajectory in the left domain is prolonged uniquely along the trajectory in the right domain. In this sense the Cauchy problem is solved and the relation a(t) is continuous (lim t→t − sing a(t) = lim t→t + sing a(t) = a sing ). The sudden type singularities (type II) are weak in Krolak's classification and they can be passed through. An example for the continuation of geodesics through the weak singularities can be found in [49]. In the case of the sudden singularity, trajectories pass through the bounce because in this singularity is the reflectional symmetryȧ → −ȧ (σ → −σ ). Therefore there is a natural identification of these points in the singularity of this type.

Sewn up singularity of type III as a model of endogenous intermediate ephemeral inflation
From Eq. (14) it is possible to detect the singularity of type III, also called 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 FRW models, which is filled with a perfect fluid with effective energy density ρ eff and pressure p eff then all singularities can be classified on the four groups [38,39]. The first class is a big rip (type I) singularity, where the energy density, pressure, and the scale factor diverge. A sudden singularity (type II) happens when the scale factor and effective energy density are finite values while pressure diverges. A big freeze singularity (type III) is observed when the effective energy density and pressure diverge at a finite value of the scale factor. A big brake (type IV) holds for the finite scale factor; the effective energy density and pressure exhibit divergences in the time derivative of the pressure or changes of the energy density rate. In that context our singularities belong to the type III.
In cosmology a singularity can occur not only at vanishing of the scale factor but even at an infinite value of the scale factor (the big rip singularity), or at a finite value of the scale factor (the sudden or big freeze singularities) [50]. They can appear in the past or in the future history of the universe. Further, singularities can be classified as being strong or weak [51][52][53][54]. The strong singularities are identified by the tidal forces which are infinite. In the case of the weak singularities we can have a divergence in the spacetime curvature but they are not strong enough to destroy arbitrarily strong detectors. For example, big bang and big rip as well as big freeze singularities are examples of the strong singularities beyond which the geodesics cannot be extended.
Different types of singularities in the FRW spacetime are considered for matter satisfying the barotropic equation of state in the form p = p(ρ) (of course other types of singularities are possible if we do not specify a form of equation of state. We do not considered them here.). The pressure in this classification is understood as an effective pressure which is a function of the effective energy density. Unlike big bang/big crunch singularities, these singularities need not occur at a vanishing scale factor, but may even occur at infinite volume. Such singularities in general relativity can be classified using the scale factor a, the energy density ρ, and the pressure p in the following types: -Type-I singularities: In this case, the scale factor diverges in finite time. The energy density and pressure diverge, causing curvature invariants to become infinite. This type of singularities is known as big rip singularities. -Type-II singularities: These events, also known as sudden singularities [55], occur at a finite value of the scale factor and at a finite time. In this case the energy density is finite, but the pressure diverges, which causes the divergence in the curvature invariants. -Type-III singularities: As in the case of type-II singularities, these singularities occur at a finite value of the scale factor, but they are accompanied by the divergence in the energy density as well as the pressure. These singularities are also known as big freeze singularities [40]. -Type-IV singularities: Unlike the case of big bang/crunch, and other singularities in this case the curvature invariants remain finite at these singularities. But their derivatives blow up [56]. In this sense, such singularities can be regarded as soft singularities. The type-IV singularities are the weakest, because they do involve any divergence of the curvature invariants. The type-I and type-III singularities, on the other hand, share the properties of big bang or big crunch singularities, in the sense that they are strong and geodesically inextensible events.
In our approach the 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 The potential as well as the effective pressure diverges. This singularity is a 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 diverges, 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 [38,40,57,58] where they are called a hyper-inflation state [59]. There exists a close relation between Palatini gravity and an effective action of Loop Quantum Gravity that reproduces the dynamics [60,61] of the considered f (R) models of polytropic spheres in the Palatini formalism f (R) = R ± λR 2 (λ is of the order of a 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 Eq. (45) simplifies to − 3K − K γ ch,0 (3A s ) The solution of the above equation is From Eqs. (15) and (47) one finds the expression for the value of the scale factor for the sewn up freeze singularity, or in the term of the redshift, This singularity, which is the horizontal inflection singularity point a = a fsing , is depicted on the diagram of the scale factor of t (Fig. 13). Note that at this singularity the potential V diverges.
For the special case of a Chaplygin gas (α = 0) the exact formulas or value of a fsing can be obtained [1]. The diagram  Fig. 14. From our numerical analysis (see Fig. 15) we have found that a single freeze type singularity is a generic prop-erty 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 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.
Taking into account the estimation of the model parameters performed in our previous paper [1] one can compute 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 see that this value is only sensitive on the γ 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. Figures 12 and 15 illustrate how values of redshifts (the scale factor marked in the figures) depend on the crucial value of the density parameter γ . While for a positive γ this is a growing function of the scale factor, 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 to 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 the redshift corresponding to the sudden singularity at the bounce. This value is 1103.67. In the case of positive γ , this parameter lies on the boundary, which disables us to formulate an analogous conclusion.
In the monograph by Capozziello and Faraoni [7, p. 73], section 3.4.2, some problems with the Palatini formalism were addressed. The authors remarked that 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 the 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 a 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. [62]) 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 the model parameters [1] for a negative γ parameter, then we find that the probability of the sudden singularity appearing after the nucleosynthesis epoch is 1-10 −16 , while the probability of one appearing before the nucleosynthesis epoch is 10 −16 . In consequence, one can reject the case with the negative values of γ even if 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 in the Palatini formulation due to reduction of the model to a dynamical system of a Newtonian type. We have found that its phase space structure is organized through the existence of two saddle points representing the static Einstein universe and the center.
Moreover, the localization of a freeze singularity during the cosmic evolution of early universe was discussed. We investigated in detail how this localization depends on a density parameters and a contribution originating from the presence of a R 2 term in the Lagrangian. The model possesses four phases: a first matter dominated deceleration phase, an intermediate accelerated phase of ephemeral inflation, a second matter dominated decelerated phase, and, finally, a second accelerating phase of current universe evolution. This evolutionary scenario becomes in agreement with the CDM model only if the freeze singularity is shifted back to the matter dominated initial singularity.
The obtained singularity 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 parameters H andḢ ). 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 [57]. In the model under consideration one deals with the situation in which freeze singularities in the past and in the future are glued together. This type of degeneration is beyond the standard classification, which requires that the acceleration parameterä is a well-defined function for which one can calculate the limits. In our case the acceleration (and in consequence the pressure) is ill defined. We will call this type of weak singularity a sewn up freeze singularity. Note that trajectories in the phase (as well as in the configuration) space pass through this singularity and continue their evolutions.
We would like to draw the reader's attention to the conclusion that the Newtonian type dynamics enables us to classify all evolutionary paths when the dynamics is reduced to the two-dimensional sewn up dynamics. 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 up 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 the conclusion that, while all open models possess a freeze singularity, there are also present bouncing models and closed ones with initial and final singularities without the freeze singularity. The latter class of models is also generic. Note that there also exists 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 configuration and phase spaces. From the classification one can conclude 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 corresponding to negative and positive values of the parameter b. There is no causal communication between them. A simple way to avoid one of the domains is by assumption 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 between matter and scalar field.
The final conclusion is that the extended cosmology with γ R 2 term in the Palatini formalism has the status of sewn dynamical systems. A point of sewing represents the freeze singularity for γ > 0 or a 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 ephemeral 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 a new singularity apart from the big bang singularity. This type of singularity is formally similar 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 of singularity can be obtained from the expansion of the function t = t (a) near an inflection point. We have t − t s ± 1 2 Therefore, a − a sing ∝ −(t sing − t) 1/2 for t → t − sing , +(t − t sing ) 1/2 for t → t + sing .
The acceleration p =ä goes to +∞ as t → t − sing and −∞ as t → t + sing . In consequence p is indefinite because of the lack of continuity in this point of gluing. Note that the system under consideration is an example of a piece-wise-smooth dynamical system [63].
Our general conclusion is that the presence of a sudden singularity in the model falsifies the negative case of γ in the Palatini cosmology. The agreement with observations which implies the evolutionary scenario for the universe requires that the value of the γ 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.
For avoiding ghost states in ETG it is sufficient to satisfy the condition f (R) > 0. In turn, to avoid the negative mass squared of a scalar-field degree of freedom (tachyon) f (R) > 0 [4] is required. From estimations of our model one can conclude that the universe without the ghost and tachyon is favored.
There are in principle two interpretations of the results obtained here. 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: considering ETG we took the simplest, quadratic correction to general relativity. Of course, we do not know the exact form of f (R) and our theory plays the role of an effective one. It is possible that after including additional higher order terms in the Taylor series expansion the singularities disappear in a natural way. Therefore, it seems to be interesting to study cosmological models with higher order terms (with respect to the Ricci scalar R or its inverse R −1 ) in the Taylor expansion for checking whether the above-mentioned singularities can appear during the cosmic evolution. It may also happen that after adding more terms into Lagrangian one can find that more consecutive freeze type singularities are allowed in the early universe, enforcing inflationary effects.
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .