Bouncing solutions from generalized EoS

We present an exact analytical bouncing solution for a closed universe filled with only one exotic fluid with negative pressure, obeying a generalized equation of state (GEoS) of the form p(ρ)=Aρ+Bρλ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\rho )=A\rho +B\rho ^{\lambda }$$\end{document}, where A, B and λ\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} are constants. In our solution A=-1/3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=-1/3$$\end{document}, λ=1/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda =1/2$$\end{document}, and B<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B<0$$\end{document} is kept as a free parameter. For particular values of the initial conditions, we find that our solution obeys the null energy condition (NEC), which allows us to reinterpret the matter source as that of a real scalar field, ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document}, with a positive kinetic energy and a potential V(ϕ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(\phi )$$\end{document}. We numerically compute the scalar field as a function of time as well as its potential V(ϕ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(\phi )$$\end{document}, and we find an analytical function for the potential that fits very accurately with the numerical data obtained. The shape of this potential can be well described by a Gaussian-type of function, and hence there is no spontaneous symmetry minimum of V(ϕ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V(\phi )$$\end{document}. We show numerically that the bouncing scenario is structurally stable in a small vicinity of the value A=-1/3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=-1/3$$\end{document}. We also include the study of the evolution of the linear fluctuations due to linear perturbations in the metric. These perturbations show an oscillatory behavior near the bouncing and approach a constant at large scales.


Introduction
Non-singular cosmologies such as described by an emergent or bouncing universe have been studied during the last decades as alternative scenarios to the inflationary paradigm, which is the one most widely accepted to describe the early universe [1,2]. Nevertheless, in inflation the problem of the initial singularity still remains [3]. On the other hand, the scale-invariant spectrum of cosmological perturbations can be obtained in most inflationary models and one natural question is if in these non-singular scenarios a scale-invariant spectrum can also be obtained. a e-mail: norman.cruz@usach.cl b e-mail: guillermo.palma@usach.cl In the case of bouncing models the universe has emerged from a cosmological bounce, where the scale factor takes a non-zero minimum value so there is no initial singularity. Bouncing universes have been investigated in a wide variety of frameworks, which include higher order theories of gravity, scalar-tensor theories and braneworlds. See [4] for a detailed discussions of different approaches to obtaining bouncing solutions. A more recent review of bouncing models driven by quantum gravity, canonical quantum cosmology, ekpyrotic scenarios, string theory, galileons, massive gravity, Horǎva-Lifshitz, among others, is presented in [5].
In this paper our aim was to find bouncing solutions for a universe only filled with one exotic fluid with negative pressure, obeying a GEoS, in the framework of general relativity. It is worthwhile mentioning that in the case of Ref. [4], examples with GEoS in the framework of general relativity are scarcely mentioned, except in the discussion of a viscous fluid as a source, where the non-singular cosmological solutions found are unstable and display non-causal propagations, in the context of the Eckart approach for non-equilibrium thermodynamics. An inhomogeneous time-dependent EoS for the dark energy has also been discussed as the corresponding fluid-like description of modified theories of gravity (see [6]). In particular, inhomogeneous viscous fluids [7][8][9] were investigated in order to obtain bouncing solutions [10][11][12], leading to a bouncing cosmology. A common method of these articles consists in proposing different Ansatzes for the scale factor a(t)-displaying a bouncing behavior-which leads as a consequence to particular equations of state.
A wide variety of cosmological models have been investigated considering a GEoS of the form where ρ is the energy density and p is the pressure of the fluid. The above factors A, B and λ are constants. In the framework of general relativity the inclusion of Eq. (1) has been used to describe the behavior of the cosmic fluid components at early and late times, as well as the possible present phantom epoch. For example, at early times and aiming to extend the range of known inflationary behaviors, Barrow [13] assumed a GEoS with A = −1 and B > 0, which corresponds to the standard EoS of a perfect fluid p = (B − 1)ρ when λ = 1. A non-singular flat universe was found for the case λ = 1/2 and B < 0, representing an emergent cosmological solution. It is interesting to mention that the doubled exponential behavior of this solution was previously found for a bulk viscous source in the presence of an effective cosmological constant [14]. This is a consequence of the inclusion of bulk viscosity in Eckart's theory, which leads to a viscous pressure of the form −3ξ H , where ξ is usually assumed to be of the form ξ = ξ 0 ρ λ . Other emergent flat solutions were found by Mukherjee et al. [15] for A > −1 and B > 0. The GEoS represented in Eq. (1) can also be seen as the sum of the standard linear EoS p = Aρ and a polytropic EoS with the polytropic exponent λ = n/(n + 1), where n is the polytropic index. Non-singular inflationary scenarios were investigated in [16] taking particular values for A, B, and n > 0.
In the study of late time evolution of the universe, there has also been assumed a GEoS of the type given by Eq. (1), motivated by the fact that the constraints from the observational data imply ω ≈ −1 for the EoS of the dark energy component, if it is described by a barotropic EoS. Nevertheless, the values ω < −1, corresponding to a phantom fluid, or ω > −1, corresponding to quintessence, cannot be discarded. Within a phenomenological approach to phantom fluids, a GEoS of the form p = −ρ − f (ρ), with f (ρ) > 0, was proposed in [8,17]. To overcome the hydrodynamic instability of a fluid with an EoS p = wρ, with w = const < 0, a general linear EoS of the form p = A(ρ − ρ 0 ) was postulated in [19], A and ρ 0 being constant and free parameters. This EoS corresponds to the particular choice λ = 0 and B = Aρ 0 and was investigated as a dark fluid filling the universe. A bouncing solution was obtained when 1 + A < 0 and Aρ 0 < 0. For a Bianchi-I cosmology, the inclusion of a perfect fluid obeying a GEoS with λ = 2 leads to a great suppression on the anisotropies in the contracting phase of a bouncing cosmology [20,21].
The case with A = −1 and λ = 1/2 was considered in [8,17,18]. In both references the cosmological solutions of dark energy models with this fluid was analyzed, focusing in the future expansion of the universe. A late time behavior of a universe filled with a dark energy component with an EoS given by Eq. (1) has been investigated in [22,23], where the allowed values of the parameters A and γ were constrained using H(z) − z data, a model-independent BAO peak parameter and a cosmic parameter (WMAP7 data).
Also theoretical studies like the so-called running vacuum energy in QFT (see [24][25][26]) give rise to a cosmological constant with a dynamical evolution during the cosmic time, which allows one to conclude that GEoS of the type of Eq. (1) could also effectively represent these scenarios under some specific assumptions.
In this paper we use a rather conservative setup, introducing a positive curvature and the particular values A = −1/3 and λ = 1/2, letting B < 0 be a free parameter of the model. With this choice the strong energy condition is violated, which is a condition to have bouncing solutions, but the NEC holds, and thus our particular GEoS has a parameter ω that evolves with cosmic time, but it lies in the range of quintessence fluids, for some choice of the initial conditions, except for t → ±∞, where the fluid behaves like a cosmological constant. These particular values of A and λ allow one to find an exact analytical bouncing solution for the scale factor, which for some initial conditions coincides with the usual de Sitter solution for a closed universe. This exact solution has the positive feature that it allows one to interpret the fluid involved by the GEoS of Eq. (1) in terms of known fluids. The importance of this result lies, contrarily to many investigations in which the scale factor a(t) is proposed as an Ansatz, from which different equations of state are derived, in the fact that we are proposing a physically well-motivated GEoS, which leads to a bouncing solution, which in addition happens to be analytic.
It has to be pointed out that the presence of a positive curvature leads to bouncing solutions with a matter content which has a strictly positive energy density. For bouncing models with a vanishing spatial curvature, the inclusion of negative energy fields is required. In the framework of general relativity it has been shown that, for a bouncing model dominated by a single exotic hydrodynamical perfect fluid, the adiabatic perturbations grow without bounds either at the bouncing point or at an exact time, when the NEC is violated or restored. Nevertheless, since NEC is not violated in our case, the solution is free of instabilities that may result, for example, from the introduction of ghost fields [27].
Reinterpreting the matter source in terms of a real scalar field, we can compute numerically the scalar field and its potential. We also found an analytical expression for this potential that fits very accurately the numerical solutions, with a coefficient of determination (r 2 ) equal to r 2 = 0.99999.
We also study the robustness of the bouncing solution when the GEoS is modified by including a perturbative term in the standard linear coefficient A. We find that under reasonable constraints on the perturbative parameter, the solution is analytic in and the first order correction allows one to extend the behavior of the bouncing solution beyond the value −1/3, for which an explicit analytic solution was found. The perturbative expansion leads us as well to conclude that the properties of the scalar potential (shape and minimum) proposed as the source for the effective equation of state are stable, provided remains small enough. This paper is organized as follows. In Sect. 2 we present the particular considered GEoS and show the analytical bouncing solution found and the main properties. In particular, we present the evolution of the parameter ω with cosmic time, discussing its quintessential behavior and how this allows one to describe the matter content by the usual real scalar field with a potential. At the end of this section we study the evolution of linear scalar perturbations to the metric. In Sect. 3 we evaluate numerically the scalar field and its potential associated to our exact solution. We also find an analytical expression for this potential that fits the numerical results found. In Sect. 4 we investigate the stability of the bouncing solution when the GEoS is modified on disturbing the parameter A by a small quantity. So in this case, we perform a study of the structural stability under small variations of the parameter A. Finally, in Sect. 5 we discuss some features and their further possible applications to suitable bouncing models.

Exact bouncing solution from GEoS
In what follows we will discuss an analytical solution for a closed universe as found in [28], for the case in which the parameter A takes the value − 1/3 and λ = 1/2 in Eq. (1). As we discuss below, the ω = p/ρ parameter can represent phantom and quintessence fluid, depending on the initial conditions. This exact solution describes a bouncing universe, assuming that one fluid with a GEoS is present in the early universe.
We consider the Friedmann-Robertson-Walker metric, and the energy-momentum tensor for a perfect fluid, with energy density ρ and pressure p, where we have chosen 8π G = 1 = c. For a universe with positive curvature (K = 1), the equation of constraint of the Friedmann equations is given by and the equation of continuity bẏ Using the change of variable s = − 2 the GEoS of Eq. (1) can be rewritten as Solving the above equations one finds the following solution: where t 0 and c are integration constants. The bouncing solution is obtained when s > 0. This solution represents a universe expanding exponentially for t ∈ (−∞, ∞). The scale factor takes a minimum value a(t = t 0 ) = s(1 − c). The positivity of the scale factor constrains c to be in the range (−∞, 1). Before we express c in terms of the initial energy density we evaluate H ,Ḣ andḦ using Eq. (7). Their expressions are the following: where For c < 0 the Hubble parameter is a strictly increasing function, so there are no critical points and we have H (t → − ∞) = − 1 s and H (t → ∞) = 1 s , then for late times this solution behaves like a de Sitter universe.
It is straightforward to evaluate the energy density as a function of the cosmic time using Eq. (4) and the expression for a(t) and H (t) given by Eq. (7) and Eq. (8), respectively. The expression for the energy density is then given by Using the initial conditions a(t 0 ) = a 0 in (7) and ρ(t 0 ) = ρ 0 in (12) we find that A dimensional restoration as regards Eq. (13) takes the following form: where v is the speed of light, R 0 is the radius of curvature, G is the gravitational constant and ρ is the energy density. Because the model considers the universe with positive curvature, the radius of the curvature R 0 also is a free parameter. A very special situation occurs for c = 0 or equivalently ρ 0 = 3/s 2 and a 0 = s, because the energy density preserves this constant value during all the cosmic evolution. It means that, for a closed universe with the EoS that we are considering, the universe expands with acceleration but the energy density remains constant, like in a de Sitter solution for a closed universe. Note that the Hubble parameter for the de Sitter solution is given by and in our case Eq. (8) with c = 0 takes a similar form: In the next subsection we give an interpretation of the used GEoS in terms of well-known fluids.

Fluid sources of the bouncing solution
We can obtain the energy density as a function of the scale factor if we substitute Eq. (7) in Eq. (12), Introducing Eq. (17) in Eq. (6) we obtain the fluid pressure as a function of the scale factor, Expanding the terms of the above expressions yields Comparing each term of the expansions our fluid can be seen as the sum of three fluids with the EoS given by So the first fluid corresponds to a cosmological constant, the second one is a quintessence term and the last corresponds to a fluid which drives an expanding universe with zero acceleration. Notice that in the above decomposition the EoS of each fluid is constant. Therefore, each ω i is independent of the parameters s, t 0 and c.
The lower plot in Fig. 1 depicts the behavior of this parameter ω.
Note that the GEoS becomes like a cosmological constant ω → −1 for t → ±∞. For the case with c = 0.5 the fluid described by the GEoS given in Eq. (21) behaves like quintessence for the lapse associated with the time of bouncing. For the case c = −0.5 the GEoS behaves like a phantom fluid within the period associated to the bouncing. We will focus on the particular interval 1 > c > 0, because in this case the GEoS leads to a quintessence-type of behavior. With c within this range, the matter content of the universe can be described by a canonical scalar field with a lagrangian minimally coupled to gravity given by where we have used the (−, +, +, +) signature. In order to reinterpret the matter source as that of a scalar field, we will evaluate the scalar field φ and the potential V (φ) in the next section.

Linear perturbations
In this section we consider the evolution of the scalar linear perturbations of the bouncing model. In what follows we use the metric notation inspired by [29,30], where the perturbed order variables α, β, ϕ and γ are scalartype perturbations; the transverse B (v) α and C (v) α are vectortype perturbations (rotation); and a transverse traceless C (t) αβ is tensor-type perturbation (gravitational wave). The energymomentum tensor is where the traceless entity π i j is the anisotropic stress tensor. Defining the variable as [31][32][33][34][35] the scalar-type perturbation of a fluid with vanishing anisotropic stress in the Einstein gravity is described by where [29,30]; ϕ v is the same as ϕ in the comoving gauge (v ≡ 0), and ϕ χ is the same as ϕ in the zero-shear gauge (a(β +aγ ) ≡ 0), etc. In a positive curvature model the wave number varies as k = (n 2 − 1)K where n is a natural not equal to 1 or 2 [37].
In the case of a minimally coupled scalar field, we have an additional nonvanishing entropic perturbation (the isotropic stress), e = (1 − c 2 s )δ v ρ). In the notation of Bardeen [36] we have We have ϕ v = R in [38]; ϕ v was also originally introduced by Lukash as − 1 3 q in [39,40]. From Eqs. (26) and (27) we can derive the following differential equation: where z = √ a/(ρ + p) and y = H/ √ (ρ + p)a. Given the explicit solutions given by Eq. (7), we solve the numerical equations (30), (28) and (26) for ϕ χ , δ v and , respectively, by using Wolfram Mathematica. The solutions are displayed  26) and (30), respectively. As no exponential grow/decay is displayed by the solutions and ϕ χ , the perturbations in this model show stability represented by an oscillatory type of behavior, close to the bouncing region. In the large scale (expanding region) and ϕ χ and δ v asymptotically approach a constant.

Computation of the scalar field and its potential
If we consider the matter content of the universe modeled by a perfect fluid, then the density and pressure in terms of the scalar field are given by Using Eqs. (6) and (12) in Eq. (31) we obtain where u is given by Eq. (11). Integrating Eq. (32), the function φ is obtained: where i is the imaginary unit, F is the elliptic integral of the first kind and is the elliptic integral of the third kind. Both are defined in [46]. The function dφ du in Eq. (32) is a continuous function. Therefore, this function has a real primitive function, but this cannot be represented by elementary functions. Thus the imaginary value in φ(u) in Eq. (34) is only a artifice of the representation of the function.
In order to numerically obtain φ(u) from the above equations, we have used standard integration subroutines from Matlab, imposing for consistence the initial condition φ(0) = 0. The result is displayed in Fig. 4, where for comparison we have also plotted its Maclaurin expansion up to order 14, whose coefficients are obtained in Appendix A. In addition, in this appendix the convergence radius of this series is shown to be arccos(c). A remarkable agreement between the two results is found, within the common range of u, which represents a severe test of the accuracy to the numerical solution.
Moreover, we have compared the results for φ(u) obtained by the Maclaurin series with the one obtained by the numerical integration by computing the Pearson coefficient r = 0.999985, which allows one to use the numerical solution beyond the convergence radius of the series expansion. We have also computed numerically V (φ) from Eq. (33) as well as its Maclaurin expansion using the expression deduced in Appendix A (see Eq. (A21)). Similarly to the analysis presented above, we have measured the degree of agreement among both methods by computing the Pearson coefficient, which in this case is 1 (r = 1). Both results are displayed in Fig. 5.
In the next subsection we will find analytical expressions for the field φ and its potential V by performing high accuracy fits.

Analytical representation of φ(u) and V (φ)
In order to characterize analytically the shape of the scalar potential and eventually to compare it with other quintessence potentials, we perform a fit of the numerical data for φ(u) using the function tanh(x), where θ 1 and θ 2 are positive parameters of the fit. The fit of field φ can be observed in Fig. 6. The fit represented by Eq. (35) is of high quality, as the corresponding coefficient of determination is r 2 = 0.99981, for the fit parameters θ 1 ≈ 3.5892 and θ 2 ≈ 0.42631.
We have also fitted the numerical data for the field φ(u) to the function θ 1 arctan(θ 2 u), where θ 1 and θ 2 are the fit parameters, but the quality of the fit was not quite comparable to the one obtained by using the analytic form given by Eq. (35).
Because of the shape of the scalar potential V (φ), we have used a Gaussian as a trial function: where σ 1 , σ 2 and σ 3 are fit parameters. The result of this fit is displayed in Fig. 7. As shown in Fig. 7, the Gaussian function fits very accurately the numerical data; moreover, the coefficient of determination r 2 = 0.99965 for the fit parameters σ 1 ≈ 7.2164, σ 2 ≈ 0.32953 and σ 3 ≈ 2.896.
Due to the bounded domain of the potential V (φ), one has to modify the Gaussian such that it falls off to zero as the field φ approaches the limits ±φ max . We have fulfilled this constraint by introducing the modified trial function: where σ 1 and σ 2 are the fit parameters. φ max turns out to be φ max ≈ 3.6009 for c = 0.5. V max and V min are the maximum and minimum values of the potential V (φ). They can explicitly be seen as follows: inserting Eqs. (17) and (18) into Eq.
Since the constants c and s are positive, the maximum (minimum) of V is obtained when a is a minimum (maximum). But from Eq. (7) the maximum and minimum values The result of the fit of V (φ) using the function defined in Eq. (37) is shown in Fig. 8.
The modified Gaussian function fits remarkably well the numerical data for V (φ) as the coefficient of determination is r 2 = 0.99999 for the values σ 1 = 0.29243 and σ 2 = 0.32674. Considering the constraint on the domain of φ and the higher accuracy of this modified Gaussian function, we conclude that the expression given by Eq. (37) is the most faithful representation of V (φ).

Analysis of the modified Gaussian potential
As a consistency check it is possible to start with the expression of Eq. (37) for the potential V (φ) and solve numerically the set of equations (38). By using standard integration subroutines this numerical strategy allows one to obtain the functions φ(t) and a(t), whose results are displayed in Figs. 9 and 10. For comparison, we have included the exact solutions given by Eqs. (34) and (7), respectively.
The two figures show a remarkable agreement between the numerical solutions and the exact expressions, which is quantified by the coefficients of determination r 2 = 0.99999 for the scalar field φ(t), and r 2 = 0.99962 for the scale factor a(t).

Structural stability of the GEoS and prevalence of the bouncing solution
Since (1). We will explore this issue taking into account a GEoS of the form where |1/3|. We study first perturbations of the exact solution found, driven by the initial GEoS of Eq. (6). We shall consider the following perturbations on the parameters a, ρ and P: where for the case = 0 the solutions are given by the Eq. (7) in the following form: Then using Eq. (4) and the continuity equatioṅ we can obtain a system to first order in , which allows one to find a (1) and ρ (1) . The system is where the prime ( ) denotes the derivative with respect to u, and the coefficients A(u), B(u), C(u), D(u), E(u), F(u) are obtained in Appendix B (see Eqs. (B4) and (B7)).
Solving the system using the given initial condition, we find that the first order contribution for the scale factor, the energy density and the pressure are given by , Now, we analyze the functions a 1 and ρ 1 of the above equation. Let us consider the first two non-zero terms of the Maclaurin expansion of Eq. (42), which represent a 0 , and the first non-zero of the Maclaurin expansion of the expression for a 1 , which is given by Eq. (46). We find that Note that a 1 (0) = 0 and that the sign of a 1 is always negative, contrary to a 0 , which is positive. For this reason, for > 0 the scale factor of the bouncing universe in a neighborhood of u = 0 grows smaller than the original solution and for < 0 the scale factor grows faster than the original one. This behavior can be seen from Eq. (40), since for > 0 the GEoS is the same that the original GEoS plus the term ρ. Therefore, the expected behavior of the new scale factor should be less pronounced than the scale factor of the original solution, because the new quintessence fluid have one state parameter ω greater than original. Finally, if we compare magnitudes of the order u 2 for a 0 and a 1 , and if in addition we include the term in a 1 , we obtain Therefore, we will have a bouncing behavior for u ≈ 0 for any couple of constants c and that satisfy Eq. (48). Note that for any c in (0, 1) there always exist > 0 and | | 1 3 which satisfy the above equation.
Now, let us analyze ρ 1 of Eq. (46). Evaluating the first nonzero terms of the Maclaurin series of ρ 0 and ρ 1 , we obtain Note that a 1 (0) = ρ 1 (0) = 0 and also that the sign of ρ 1 is negative for c < 0.5 and positive for c > 0.5. In the case c = 0.5 the first term, ρ 1 = 0. Thus, for any | | 1 3 , ρ 0 will dominate over ρ 1 for u ≈ 0. Now, if we compare magnitudes of the order u 2 for ρ 0 and ρ 1 , and if in addition we include the term in ρ 1 , we find that for c = 0.5 3c This equation tell us what values of c and lead ρ 0 to dominate over ρ 1 for u ≈ 0.
The behavior of the original scale factor and the perturbed ones are shown in Fig. 11. The that we consider in the graphics corresponds to a 45% of 1/3 and c = 0.5. With these c and we find that Eq. (48) is satisfied and, therefore, there is a bouncing behavior.
The coefficient of determination is r 2 = 0.94759, for a 0 and a 0 + a 1 , and r 2 = 0.99875, for a 0 and a 0 − a 1 . The energy density ρ 1 given by Eq. (46) is shown in Fig. 12.
The coefficient of determination is r 2 = 0.99980, for ρ 0 and ρ 0 + ρ 1 and r 2 = 0.98566, for ρ 0 and ρ 0 − ρ 1 . We may note that despite the high value of , the behavior of the scale factor and of the energy density perturbed to first order are quite similar to the one obtained with the exact solution. Moreover, if we change the parameters c and s in their respective domain we will get universes keeping the same shape of the bouncing but with different growths.
We evaluate numerically the first order perturbation of the field φ, which we denoted by φ 1 (u). Its behavior is showed in Fig. 13. Like the scale factor and the energy density, we find that the first order perturbed field φ 1 preserves the shape of the unperturbed solution φ 0 . The coefficient of determination is r 2 = 0.99999, for φ 0 and φ 0 + φ 1 , and r 2 = 0.99998 for φ 0 and φ 0 − φ 1 .
To obtain the potential V 1 we expand the function V (φ) in . From this expansion we obtain the following expression: The potential V (φ) is plotted in Fig. 14.
We obtain a r 2 = 0.99975 for V 0 and V 0 + V 1 and a r 2 = 0.99975, for V 0 and V 0 − V 1 .
The perturbation of the GEoS in the parameter A around the value −1/3 leads also to bouncing solutions whose behavior in the scale factor and in the energy density are quite similar to the exact analytical solution found in this work. This also applied for the scenario in which a scalar field describes the matter content of the universe.

Conclusions
The study of GEoS has been very important in the exploration of new scenarios in the very early phases of the universe like inflation, or bouncing universe theories in which there are no initial singularities.
In this work we have found an exact analytical bouncing solution for a closed universe filled with one fluid which obeys a GEoS of the form p = −1/3ρ + Bρ 1/2 , where B < 0 is a free parameter. We have chosen the initial conditions that allow no violation of NEC, which leads the parameter ω to evolve with the cosmic time in the domain of quintessence. For t → ±∞, ω → −1, and the fluid behaves like a cosmological constant. We have also shown that the well-known de Sitter solution with positive curvature is obtained as a particular case of our exact analytical bouncing solution.
Another interesting feature, derived from the fact that we have an exact solution, is the possibility to interpret the fluid described by the GEoS of Eq. (1) in terms of known fluids. In fact, expanding the expressions for the pressure, p, and the energy density, ρ, in terms of the scale factor and inserting them into the GEoS, we found that the matter content can be seen as the sum of the contributions coming from three fluids: a cosmological constant, quintessence (ω = − 2/3) and the corresponding fluid which arises from the particular value ω = − 1/3.
Since the investigated fluid behaves effectively like quintessence, it is possible to reinterpret the matter source in terms of an ordinary scalar field φ, minimally coupled to gravity with a positive kinetic term and a potential V (φ). Moreover, since NEC is not violated, the solution is free of instabilities at the perturbative level.
We have solved for φ and V (φ) the set of coupled equations (31) by using Maclaurin expansions up to 14th order and, in parallel, as an accuracy test, we have computed them numerically from the implicit equations (32) and (33). Remarkable agreement was found by comparing the scalar field φ and the potential V (φ) obtained by the two methods, whose precision is characterized by coefficients of determination r 2 = 0.99997 and r 2 = 1, respectively (see Figs. 4 and 5), holding in the common range.
Performing a high accuracy fit we have found an analytical expression for the field and its scalar potential with coefficients of determination typically of the order of 1 up to 10 −4 and 10 −5 , respectively (see Eqs. (35) and (37). The shape of this potential can be very precisely described by a Gaussian-type of function that has a bounded domain given by the condition φ min < φ < φ max . With this exact analytical expression for the scalar potential we evaluated numerically the scale factor, finding a curve very close to the exact analytical bouncing solution already found. This represents a very stringent test on the accuracy of the numerical method used to compute the scalar potential V (φ).
We have also studied the structural stability of the analytical bouncing solution when the GEoS is modified by including a perturbative term in the standard linear coefficient A (see Eq. 40). For sufficiently small , the scenario predicted by the analytic solution is still preserved as the scale factor a and the density ρ behave quite similar to the unperturbed solutions (see Figs. 11 and 12). The shape of the scalar potential -introduced as a possible source to generate the effective GeoS -also confirms the unperturbed scenario: the absence of a spontaneous symmetry breaking minimum, for small enough within the validity range of the first order approximation.
In summary, the exact analytical bouncing solution can be extended to a vicinity of A = −1/3, confirming a bouncing scenario beyond the particular value required for the exact solution. Moreover, an analytical quintessence potential has been found by using a high accuracy fit to the numerical data. A scalar field theory minimally coupled to gravity and ruled by this potential leads to bouncing solutions for closed universes, which does not present spontaneous symmetry breaking.
We have studied the time evolution of the scalar-type perturbations in terms of the variables , ϕ χ and δ v . We have found that an oscillatory behavior characterizes both and ϕ χ in a region close to the bouncing, while a constant behavior is displayed for large scales. In the case of δ v , it keeps essentially constant on the whole scale. As a consequence, the scalar perturbations in this model do not lead to instability.
An analytic scalar potential was found in this article, and it can further be used to study perturbations in the background metric in the trivial minimum and in the absence of spontaneous symmetry breaking.
From these we obtain for l ≥ 1 and n ≥ 1 Moreover, using Eqs. (A5), (A8) and (A10), we find where where n ≥ 1, and Finally, using the above results for f (u) we obtain (1−c) 2 , α 4 = 0. (A15) Using the fact that the convergence radius of a power series of an analytic function is the distance from the center of the power series to the closest singularity and that 0 < c < 1, we conclude that the convergence radius of Maclaurin series of φ(u) of Eq. (A14) is equal to arccos(c). Now, we will compute the coefficients β n of the Maclaurin series of V (φ). To this aim, we use the coefficients c n , which can be obtained in a similar way to how the coefficients α n were computed, since we can rewrite Eq. (33) as A straightforward computation yields where Y (2n) = 2c 5 n−1 k=1 2n 2k g 2(n−k) g 2k + 1 + 4c 5(1−c) g 2n , with Y (2n−1) = 0 for n ≥ 1, and the constants g 2n given by with