Bouncing solutions from generalized EoS

TL;DR

This paper derives an exact analytical bouncing cosmological solution in a closed universe with a single exotic fluid obeying a generalized equation of state, demonstrating stability and scalar field interpretation.

Contribution

It introduces a new exact bouncing solution with a specific generalized equation of state and analyzes its stability and scalar field potential.

Findings

Solution obeys Null Energy Condition under certain initial conditions

Scalar field potential approximates a Gaussian function

Bouncing solution is stable under small parameter variations

Abstract

We present an exact analytical bouncing solution for a closed universe filled with only one exotic fluid with negative pressure, obeying a Generalized Equations of State (GEoS) of the form $P(\rho)=A\rho+B\rho^{\lambda}$, where $A$, $B$ and $\lambda$ are constants. In our solution $A=-1/3$ and $\lambda=1/2$ and $B<0$ is kept as a free parameter. For particular values of the initial conditions, we obtain that our solution obeys Null Energy Condition (NEC), which allows us to reinterpret the matter source as that of a real scalar field, $\phi$, with a positive kinetic energy and a potential $V(\phi)$. We compute numerically the scalar field as a function of time as well as its potential $V(\phi)$, and find an analytical function for the potential that fits very accurately with the numerical results 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(\phi)$. We further show that the bouncing scenario is structurally stable under small variations of the parameter $A$, such that a family of bouncing solutions can be find numerically, in a small vicinity of the value $A=-1/3$.

Equations

$$P(\rho )=A\rho +B{\rho }^{\lambda },$$

$$\rho =3{\left(\frac{\dot{a}}{a}\right)}^2+\frac{3}{a^2},$$

$$\dot{\rho }+3H(\rho +p)=0.$$

$$p(\rho )=-\frac{1}{3}\rho -\frac{2}{s\sqrt{3}}{\rho }^{1/2}.$$

$$a(t)=s\left[\cosh \left(\frac{t-t_0}{s}\right)-c\right],$$

$$H(u)=\frac{\sinh (u)}{s(\cosh (u)-c)},\text{vspace}2mm$$

$$\dot{H}(u)=\frac{1-c\cosh \left(u\right)}{{\left[s(\cosh \left(u\right)-c)\right]}^2},\text{vspace}2mm$$

$$\ddot{H}(u)=-\frac{\sinh (u)\left(2-c^2\right)+\sinh (2u)\left(\frac{r}{3}+\frac{1}{s}\right)}{{\left[s(\cosh (u)-c)\right]}^3},$$

$$u=\frac{t-t_0}{s}.$$

$$\rho =\frac{3{\cosh }^2(u)}{{\left[s(\cosh (u)-c)\right]}^2}.$$

$$\begin{array}[]{c}a_{0}=s(1-c),\quad\ \ 3=a_{0}^{2}\rho_{0}.\end{array}$$

$$\begin{array}[]{c}a_{0}=s(1-c)\dfrac{v^{2}}{R_{0}\sqrt{8\pi G}},\quad\ \ \dfrac{v^{2}}{R_{0}^{2}}=\dfrac{8\pi G}{3}\rho_{0}a_{0}^{2},\end{array}$$

$$H(t)=\sqrt{\frac{\lambda }{3}}\tanh \left(\sqrt{\frac{\lambda }{3}}(t-t_0)\right),$$

$$H(t)=\frac{1}{s}\tanh \left(\frac{t-t_0}{s}\right).$$

$$\rho (a)=3{\left(\frac{1}{s}+\frac{c}{a}\right)}^2.$$

$$p(a)=-{\left(\frac{1}{s}+\frac{c}{a}\right)}^2-\frac{2}{s}\left(\frac{1}{s}+\frac{c}{a}\right).$$

$$\rho (a)=\frac{3}{s^2}+\frac{6c}{sa}+\frac{3c^2}{a^2}={\rho }_1+{\rho }_2+{\rho }_3,\text{vspace}2mm$$

$$p(a)=-\frac{3}{s^2}-\frac{4c}{sa}-\frac{c^2}{a^2}=p_1+p_2+p_3.$$

$$\omega =-1+\frac{2c}{3\cosh \left(\frac{t-t_0}{s}\right)}.$$

$$\mathcal{L}=\frac{1}{2}{\varphi }_{,\mu }{\varphi }^{,\mu }-V(\varphi ).$$

$$\begin{array}[]{c}\rho=\dfrac{1}{2}\dot{\phi}^{2}+V(\phi),\quad p=\dfrac{1}{2}\dot{\phi}^{2}-V(\phi).\\ \\ \end{array}$$

$$\frac{d\varphi }{du}=\pm \frac{\sqrt{2c\cosh (u)}}{\cosh (u)-c},\text{vspace}2mm$$

$$V(u)=\frac{\cosh (u)(2\cosh (u)-c)}{{\left[s(\cosh (u)-c)\right]}^2},$$

$$\varphi (\text{u})=-i\sqrt{8c}\left[F\left(\frac{i\text{u}}{2},\sqrt{2}\right)+\frac{c}{1-c}\Pi \left(\frac{i\text{u}}{2},\frac{2}{1-c},\sqrt{2}\right)\right],$$

$$\varphi (u)={\theta }_1\cdot \tanh ({\theta }_2\cdot u),$$

$$V(\varphi )={\sigma }_1\cdot \exp \left(-{\sigma }_2\cdot {\varphi }^2\right)+{\sigma }_3,$$

$$V(\varphi )=(V_{max}-V_{min})\exp \left[-{\sigma }_1{\varphi }^2-\frac{{\sigma }_2{\varphi }^2}{{\varphi }_{max}^2-{\varphi }^2}\right]+V_{min}$$

$$\begin{array}[]{l}\dot{\phi}(a)=\dfrac{2c}{sa}+\dfrac{2c^{2}}{a^{2}},\quad\ \ V(a)=\dfrac{3}{s^{2}}+\dfrac{5c}{sa}+\dfrac{2c^{2}}{a^{2}}.\end{array}$$

$$\begin{array}[]{l}V_{max}=\dfrac{3-c}{s^{2}(1-c)^{2}},\quad\quad V_{min}=\dfrac{3}{s^{2}}.\end{array}$$

$$p=\left(-\frac{1}{3}+ϵ\right)\rho -\frac{2}{s\sqrt{3}}{\rho }^{1/2},$$

$$\begin{array}[]{l}a(u)=a_{(0)}(u)+\epsilon a_{(1)}(u)+O(\epsilon^{2}),\vspace{2mm}\\ \rho(u)=\rho_{(0)}(u)+\epsilon\rho_{(1)}(u)+O(\epsilon^{2}),\vspace{2mm}\\ p(u)=p_{(0)}(u)+\epsilon p_{(1)}(u)+O(\epsilon^{2}),\end{array}$$

$$\begin{array}[]{l}a_{(0)}(u)=s\left(\cosh(u)-c\right),\vspace{2mm}\\ \rho_{(0)}(u)=\dfrac{3(\cosh(u))^{2}}{s^{2}(\cosh(u)-c)^{2}},\vspace{2mm}\\ p_{(0)}(u)=-\dfrac{1}{3}\rho_{(0)}-\dfrac{2}{s\sqrt{3}}\rho_{(0)}^{1/2}.\end{array}$$

Full text

I Introduction

Non singular cosmologies such as the described by an emergent or bouncing universe have been studied during the last decades as alternative scenarios to the inflationary paradigm, which is the most accepted one to describe the early universe Lemoine , Martin . Nevertheless, in inflation the problem of the initial singularity still remains Borde . 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 an scale-invariant spectrum can also be obtained.

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 includes among others, higher order theories of gravity, scalar-tensor theories and braneworlds. See Novello for a detailed discussions of different approaches to obtain bouncing solutions.

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. A wide variety of cosmological models have been investigated considering a GEoS of the form

[TABLE]

where $A$, $B$ and $\lambda$ 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  Barrow assumed a GEoS with $A=-1$ and $B>0$, which corresponds to the standard EoS of a perfect fluid $p=(B-1)\rho$ when $\lambda=1$. A non singular flat universe was found for the case $\lambda=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 Barrow1 . This is a consequence of the inclusion of bulk viscosity in the Eckart s theory, which leads to a viscous pressure $\Pi$ of the form $-3\xi H$, where $\xi$ is assumed usually in the form $\xi=\xi_{0}\rho^{\lambda}$. Other emergent flat solutions were found by Mukherjee et al Mukherjee 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\rho$ and a polytropic EoS with the polytropic exponent $\lambda=n/(n+1)$, where $n$ is the polytropic index. Non singular inflationary scenarios were investigated in Chavanis taken particular values for $A$, $B$, and $n>0$.

In the study of late time evolution of the universe, it has been also assumed GEoS of the type given by Eq.(1), motivated by the fact that the constraints from the observational data implies $\omega\approx-1$ for the EoS of the dark energy component, if it is ruled by a barotropic EoS. Nevertheless, the values $\omega<-1$, corresponding to a phantom fluid, or $\omega>-1$, corresponding to quintessence can not be discarded. Within a phenomenological approach to phantom fluids, a GEoS of the form $p=-\rho-f(\rho)$, with $f(\rho)>0$, was proposed in Odintsov . To overcome the hydrodynamic instability of a fluid with an EoS $p=w\rho$, with $w=const<0$, a general linear EoS of the form $p=A(\rho-\rho_{0})$ was postulated in Babichev , being $A$ and $\rho_{0}$ constant and free parameters. This EoS corresponds to the particular choice $\lambda=0$ and $B=A\rho_{0}$ and was investigated as a dark fluid filling the universe. A bouncing solution was obtained when $1+A<0$ and $A\rho_{0}<0$. For a Bianchi-I cosmology, the inclusion of a perfect fluid obeying a GEoS with $\lambda=2$ leads to a great suppression on the anisotropies in the contracting phase of a bouncing cosmology Bozza .

The case with $A=-1$ and $\lambda=1/2$ was considered in  Odintsov and Stefancic . In both works 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  Paul ,  Paul1 , where the allowed values of the parameters $A$ and $\gamma$ were constrained using H(z)-z data, a model independent BAO peak parameter and cosmic parameter (WMAP7 data).

Also theoretical studies like the so called running vacuum energy in QFT (see  Shapiro ) gives rise to a cosmological constant with a dynamical evolution during the cosmic time, which allows to conclude that GEoS of the type of Eq.(1) could also effectively represent these scenarios under some specific assumptions.

In this work we use a rather conservative setup introducing a positive curvature and the particular values $A=-1/3$ and $\lambda=1/2$, letting $B<0$ as a free parameter of the model. With this election the strong energy condition is violated, which is a condition to have bouncing solutions, but NEC holds, and thus our particular GEoS has a parameter $\omega$ that evolves with the cosmic time, but lies in the range of quintessence fluids, for some choice of the initial conditions, except for $t\rightarrow\pm\infty$, where the fluid behaves like a cosmlogical constant. These particular values of $A$ and $\lambda$ allow to find an exact analytical bouncing solution for the scale factor.

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 fit 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 $\epsilon$ and the first order correction allows 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 as well to conclude that the properties of the scalar potential (shape and minimum) proposed as source for the effective equation of state are stable, provided $\epsilon$ remains small enough.

This paper is organized as follows. In section II we present the particular considered GEoS and show the analytical bouncing solution found and their main properties. In particular, we present the evolution of the parameter $\omega$ with the cosmic time, discussing its quintessential behavior and how this allow to describe the matter content by a usual real scalar field with a potential. In section III 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 section IV we investigate the stability of the bouncing solution when the GEoS is modified disturbing the parameter $A$ by an small quantity. So in this case, we make a study of structural stability under small variations of the parameter $A$. Finally, in section V we discuss some features and their further possible applications to suitable bouncing models.

II Exact bouncing solution from GEoS

In what follows we will discuss an analytical solution for a closed universe found in Contreras , for the case in which the parameter $A$ takes the value $-1/3$ and $\lambda=1/2$ in Eq.(1). As we discuss bellow the $\omega=p/\rho$ parameter can represent phantom and quintessence fluids, depending on the initial conditions. This exact solution describe a bouncing universe, assuming that one fluids with a GEoS is present in the early universe.

For a universe with positive curvature ($k=1$), the equation of constraint of the Friedmann equations is given by

[TABLE]

and the equation of continuity by

[TABLE]

Using the change of variable $s=-\dfrac{2}{B\sqrt{3}}$ the GEoS of Eq. (1) can be rewritten as

[TABLE]

Solving the above equations one finds the following solution

[TABLE]

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\in(-\infty,\infty)$. The scale factor takes a minimum value $a(t=t_{0})=s(1-c)$. The positivity of the scale factor constraints $c$ to be in the following range $(-\infty,1)$. Before to express $c$ in terms of the initial energy density we evaluate $H$, $\dot{H}$ and $\ddot{H}$ using the Eq.(5). Their expressions are the following:

[TABLE]

where

[TABLE]

For $c<0$ the Hubble parameter is a strictly increasing function, so there are no critical points and we have $H(t\rightarrow-\infty)=-\frac{1}{s}$ and $H(t\rightarrow\infty)=\frac{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.(2) and the expression for $a(t)$ and $H(t)$ given by Eq.(5) and Eq.(6), respectively. The expression for the energy density is then given by

[TABLE]

Using the initial conditions $a(t_{0})=a_{0}$ in (5) and $\rho(t_{0})=\rho_{0}$ in (10) we obtain that

[TABLE]

One dimensional restoration lead to the Eq.(11) takes the following form:

[TABLE]

where $v$ is the speed of light, $R_{0}$ is the radius of curvature, $G$ is the gravitational contant and $\rho$ is the energy density. Because the model considers the universe with curvature positive, the radius the curvature $R_{0}$ also is a free parameter.

A very special situation occurs for $c=0$ or equivalently $\rho_{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 expand with acceleration but the energy density remains constant, like in de Sitter solution for a closed universe. Note that the Hubble parameter for the de Sitter solution is given by

[TABLE]

and in our case the Eq.(6) with $c=0$ takes the similar form:

[TABLE]

In the next subsection we give an interpretation of the used GEoS in terms of well known fluids.

II.1 Fluid sources of the bouncing solution

We can obtain the energy density as a function of the scale factor if we replace the Eq.(5) in the Eq.(10)

[TABLE]

Introducing Eq.(15) in Eq.(4) we obtain the fluid pressure as a function of the scale factor

[TABLE]

Expanding the terms of the above both expressions yields:

[TABLE]

Comparing each terms of the expansions our fluid can be seen as the sum of three fluids with the EoS given by $\omega_{1}=p_{1}/\rho_{1}=-1$, $\omega_{2}=p_{2}/\rho_{2}=-2/3$ and $\omega_{3}=p_{3}/\rho_{3}=-1/3$, respectively. So the first fluid corresponds to a cosmological constant, the second one is a quintessence and the last corresponds to a fluid which drives an expanding universe with zero acceleration. Notice that in the above descomposition the EoS of each fluid is constant. Therefore, each $\omega_{i}$ is independent of the parameters $s$, $t_{0}$ and $c$.

Lets us evalue the EoS, $\omega=p/\rho$, for this fluid, which in terms of the cosmic time takes the expression

[TABLE]

The lower plot in Fig. 1. depicted the behavior of this parameter $\omega$.

Note that the EoS becomes like a cosmological constant $\omega\rightarrow-1$ for $t\rightarrow\pm\infty$. For the case with $c=0.5$ the fluid ruled by the EoS given in Eq.(19) behaves like quintessence for the lapse associated at the time of bouncing. For the case $c=-0.5$ the EoS 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 real scalar field with a lagrangian minimally coupled to gravity given by

[TABLE]

In order to reinterpret the matter source as that of a scalar field, we will evaluate the scalar field $\phi$ and the potential $V(\phi)$ in the next section.

III 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 term of the scalar field are given by

[TABLE]

Using the Eq.(4) and the Eq.(10) in the Eq.(21) we obtain:

[TABLE]

where $u$ is given by the Eq.(9). Integrating the Eq.(22), the function $\phi$ is obtained as:

[TABLE]

where $i$ is the imaginary unit, $F$ is the elliptic integral of the first kind and $\Pi$ is the elliptic integral of the third kind. Both are defined in Gradshteyn .

The function $d\phi\over du$ in Eq.(22) is a continuous function. Therefore, this function have a real primitive function, but this can’t be represented by elementary functions. Thus the imaginary value in $\phi(u)$ in the Eq.(24) is only a artifice of the representation of the function.

In order to numerically obtain $\phi(u)$ from the above equations, we have used standard integration subroutines from Matlab, imposing for consistence the initial condition $\phi(0)=0$. The result is displayed in Fig. 2, where for comparison we have also plotted its Maclaurin expansion up to order 14th, whose coefficients are obtained in Appendix A. In addition, in this appendix the convergence radius of this series shown to be $\arccos(c)$. A remarkable agreement among both results is found, within the common range of $u$, which represents a severe test of accuracy to the numerical solution.

Moreover, we have compared the results for $\phi(u)$ obtained by the Maclaurin series with the one obtained by the numerical integration by computing the Pearson’s coefficient $r=0.999985$, which allows to use the numerical solution beyond the convergence radius of the series expansion.

We have also computed numerically $V(\phi)$ from Eq.(23) as well as its Maclaurin expansion using the expression deduced in Appendix A (see the Eq.(A21)). Similarly to the analysis explained above, we have measured the degree of agreement among both methods by computing the Pearson’s coefficient, which in this case is 1 ($r=1$). Both results are displayed in Fig. 3.

In the next subsecction we will find analytical expressions for the field $\phi$ and its potencial $V$ by performing high accuracy fits.

III.1 ** Analytical representation of $\phi(u)$ and $V(\phi)$**

In order to characterize analytically the shape of the scalar potential, and eventually to compare it with other quintaessence potentials, we perform a fit of the numerical data for $\phi(u)$ using the function $\tanh(x)$,

[TABLE]

where $\theta_{1}$ and $\theta_{2}$ are paremeters positives of fit. The fit of field $\phi$ can be observed in the Fig. 4.

The fit represented by Eq.(25) is of high quality as the corresponding coefficient of determination is $r^{2}=0.99981$, for the fit parameters $\theta_{1}\approx 3.5892$ and $\theta_{2}\approx 0.42631$.

We have also fitted the numerical data for the field $\phi(u)$ to the function $\theta_{1}\arctan(\theta_{2}u)$, where $\theta_{1}$ and $\theta_{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.(25).

Because of the shape of the scalar potential $V(\phi)$, we have used a Gaussian as a trial function :

[TABLE]

where $\sigma_{1}$, $\sigma_{2}$ and $\sigma_{3}$ are fit parameters. The result of this fit is displayed in Fig. 5.

As it is shown in Fig. 5, the Gaussian function fits very accurately the numerical data, moreover the coefficient of determination $r^{2}=0.99965$ for the fit parameters $\sigma_{1}\approx 7.2164$, $\sigma_{2}\approx 0.32953$ and $\sigma_{3}\approx 2.896$.

Due to the bounded domain of the potential $V(\phi)$, one has to modify the Gaussian such that it falls off to zero as the field $\phi$ approaches the limits $\pm\phi_{max}$. We have fulfilled this constraint by introducing the modified trial function:

[TABLE]

where $\sigma_{1}$ and $\sigma_{2}$ are the fit parameters. $\phi_{max}$ turns out to be $\phi_{max}\approx 3.6009$. $V_{max}$ and $V_{min}$ are the maximum and minimum values of the potential $V(\phi)$. They can explicitly be obtained as follows: inserting Eqs. (15) and (16) into Eq.(21) one obtains:

[TABLE]

Since the constants $c$ and $s$ are positive, the maximum (minimum) of $V$ is obtained when $a$ is a minimum (maximum). But from Eq.(5) the maximum and minimum values of $a(t)$ are $\infty$ and $s(1-c)$ respectively. Inserting these values into Eq.(28) we obtain:

[TABLE]

The result of the fit of $V(\phi)$ using the function defined in Eq.(27) is shown in Fig. 6.

The modified Gaussian function fits remarkable well the numerical data for $V(\phi)$ as the coefficient of determination is $r^{2}=0.99999$ for the values $\sigma_{1}=0.29243$ and $\sigma_{2}=0.32674$. Considering the constraint on the domain of $\phi$ and the higher accuracy of this modified Gaussian function, we conclude that the expression given by Eq.(27) is the faithfuliest representation of $V(\phi)$.

III.2 Analysis of the modified Gaussian potential

As a consistency check it is possible to start with the expression of Eq.(27) for the potential $V(\phi)$ and solve numerically the set of equations (28). By using standard integration subroutines this numerical strategy allows to obtain the functions $\phi(t)$ and $a(t)$, whose results are displayed in Figs. 7 and 8, where for comparison, we have included the exact solutions given by Eqs. (24) and (5) respectively.

Both figures show a remarkable agreement among the numerical solutions and the exact expressions, which is quantified by the coefficients of determination $r^{2}=0.99999$ for the scalar field $\phi(t)$, and $r^{2}=0.99962$ for the scale factor $a(t)$.

IV Structural stability of the GEoS and prevalence of the bouncing solution

Since and exact bouncing solution is obtained for the very particular value $A=-1/3$, it is an open question whether variations of the parameter $A$ leads to structural stability of the field equations, or in other words whether the system still have bouncing solutions when the GEoS is modified by disturbing the parameter $A$ by $A\rightarrow A+\epsilon$ in the Eq.(1). We will explore this issue taking into account a GEoS of the form

[TABLE]

where $\epsilon\ll|1/3|$. We study first perturbations on the exact solution found, driven by the initial GEoS of Eq.(4). We shall consider the following perturbations on the parameters $a$, $\rho$ and $P$

[TABLE]

where for the case $\epsilon=0$ the solutions are given by the Eq.(5) in the following form

[TABLE]

Then using Eq.(2) and the continuity equation

[TABLE]

we can obtain a system to first order in $\epsilon$, which allows to find $a_{(1)}$ and $\rho_{(1)}$. The system is:

[TABLE]

where the “prime” ( {}_{{}^{\prime}}\) denotes the derivative with respect to $u$, and the cofficients $A(u)$, $B(u)$, $C(u)$, $D(u)$, $E(u)$, $F(u)$ are obtain in the Appendix $B$ (see the Eqs.(B4) and (B7)).

Solving the system using the given initial condition,

[TABLE]

we obtain that the first order contribution for the scale factor, the energy density and the pressure are given by

[TABLE]

Now, we analyze the functions $a_{1}$ and $\rho_{1}$ of above equation. Let us consider the first two nonzero terms of the Maclaurin expansion of Eq.(32), which represents $a_{0}$, and the first nonzero nonzero of the Maclaurin expansion of the expression for $a_{1}$, which is given by Eq.(36). We obtain that

[TABLE]

Note that $a_{1}(0)=0$ and that the sign of $a_{1}$ is always negative contrary to $a_{0}$ that is positive. For this reason for $\epsilon>0$ the scale factor of the bouncing universe in a neighborhood of $u=0$ grows lesser than the original solution and for $\epsilon<0$ the scale factor grows fasther than the original one. This behavior can be seen from Eq.(30), since for $\epsilon>0$ the GEoS is the same that the original GEoS plus the term $\epsilon\rho$. Therefore, the expected behavior of the new scale factor should be less pronunced that the scale factor of the original solution, because the new quintessence fluid have one state parameter $\omega$ greater than original. Finally, if we capare magnitudes of the order $u^{2}$ for $a_{0}$ and $a_{1}$, and if in addition we include the term $\epsilon$ in $a_{1}$, we obtain that

[TABLE]

Therefore, we will have a bouncing behavior for $u\approx 0$ for any couple of costants $c$ and $\epsilon$ that satisfy the Eq.(38). Note that for any $c$ between $(0,1)$ always exist one $\epsilon>0$ and $|\epsilon|\ll{1\over 3}$ which satisfy the above equation.

Now, let’s analyze $\rho_{1}$ of Eq.(36). Evaluating the first nonzero terms of Maclaurin series of $\rho_{0}$ and $\rho_{1}$, we obtain

[TABLE]

Note that $a_{1}(0)=\rho_{1}(0)=0$ and also that the sign of $\rho_{1}$ is negative for $c<0.5$ and positive for $c>0.5$. In the case $c=0.5$ the first $\rho_{1}=0$. Thus, for any $|\epsilon|\ll{1\over 3}$, $\rho_{0}$ will dominate over $\rho_{1}$ for $u\approx 0$. Now, if we compare magnitudes of the order $u^{2}$ for $\rho_{0}$ and $\rho_{1}$, and if in addition we include the term $\epsilon$ in $\rho_{1}$, we obtain that for $c\neq 0.5$

[TABLE]

This equation tell us what values of $c$ and $\epsilon$ lead $\rho_{0}$ dominates over $\rho_{1}$ for $u\approx 0$.

The behavior of original scale factor and the perturbed ones are shown in Fig. 9. The $\epsilon$ that we consider in the graphics corresponds to a $45\%$ of $1/3$ and $c=0.5$. With these $c$ and $\epsilon$ we obtain that the Eq.(38) is satisfied and therefore, there is a bouncing behavior.

The coefficient of determination is $r^{2}=0.94759$, for $a_{0}$ and $a_{0}+\epsilon a_{1}$, and $r^{2}=0.99875$, for $a_{0}$ and $a_{0}-\epsilon a_{1}$. The energy density $\rho_{1}$ given by Eq.(36) is shown in Fig. 10.

The coefficient of determination is $r^{2}=0.99980$, for $\rho_{0}$ and $\rho_{0}+\epsilon\rho_{1}$ and $r^{2}=0.98566$, for $\rho_{0}$ and $\rho_{0}-\epsilon\rho_{1}$. We can note that despite of the high value of $\epsilon$, the behavior of scale factor and of the energy density perturbed to first order are quite similar to the 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 $\phi$, which we denoted by $\phi_{1}(u)$. Its behavior is showed in the Fig. 11.

Like the scale factor and the energy density, we obtain that the first order perturbed field $\phi_{1}$ preserve the shape of unperturbed solution $\phi_{0}$. The coefficient of determination is $r^{2}=0.99999$, for $\phi_{0}$ and $\phi_{0}+\epsilon\phi_{1}$, and $r^{2}=0.99998$ for $\rho_{0}$ and $\rho_{0}-\epsilon\rho_{1}$.

To obtain the potential $V_{1}$ we exapand the function $V(\phi)$ in $\epsilon$. From this expansion we obtain the following expression

[TABLE]

The potential $V(\phi)$ is plotted in the Fig. 12.

We obtain a $r^{2}=0.99975$ for $V_{0}$ and $V_{0}+\epsilon V_{1}$ and a $r^{2}=0.99975$, for $V_{0}$ and $V_{0}-\epsilon V_{1}$.

The perturbation of the GEoS in the parameter $A$ around the va ue $-1/3$ leads also to bouncing solutions whose behavior in the scale factor and in the energy density are quite similar than the exact analytical solution found in this work. This also applied for the scenario in which a scalar field describe the matter content of the universe.

V 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 and exact analytical bouncing solution for a closed universe filled with one fluid which obeys a GEoS of the form $p=-1/3\rho+B\rho^{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 $\omega$ to evolve with the cosmic time in the domain of quintessence. For $t\rightarrow\pm\infty$, $\omega\rightarrow-1$ 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 of our result is the possibi-lity to interpret the fluid ruled by the GEoS of Eq.(1) in terms of known fluids. In fact, expanding the expressions for the pressure, $p$, and the energy density, $\rho$, in terms of the scale factor and inserting them into the GEoS, we obtained that the matter content can be seen as the sum of the contributions coming from three fluids: a cosmological constant, quintessence ($\omega=-2/3$) and the corresponding fluid which arises from the particular value $\omega=-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 $\phi$, minimally coupled to gravity with a positive kinetic term and a potential $V(\phi)$.

We have solved for $\phi$ and $V(\phi)$ the set of coupled equations (21) by using Maclaurin expansions up to 14th order and, in paralell as an accuracy test, we have computed them numerically from the implicit Eqs.(22) and (23). A remarkable agreement was found by comparing the scalar field $\phi$ and the potential $V(\phi)$ obtained by both methods, whose precision is characterized by coefficients of determination $r^{2}=0.99997$ and $r^{2}=1$ respectively, (see Figs. (2) and (3)), which holds 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 one up to $10^{-4}$ and $10^{-5}$ respectively (see Eqs. (25) and (27). 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 $\phi_{min}<\phi<\phi_{max}$. With this exact analytical expression for the scalar potential we evaluated numerically the scale factor, founding a curve very close to the exact analytical bouncing solution already found. This represents a very astringent test on the accuracy of the numerical method used to compute the scalar potential $V(\phi)$.

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.(30)). For sufficiently small $\epsilon$, the scenario predicted by the analytic solution is still preserved as the scale factor $a$ and the density $\rho$ behave quite similar to the unperturbed solutions (see Eqs. (38) and (40)). 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 $\epsilon$ 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.

The analytic scalar potential found in this article can further be used to study other interesting issues associated, like the consequences of considering perturbations in the background metric in the trivial minimum and absence of spontaneous symmetry breaking .

Acknowledgements

This work was supported by CONICYT through Grant FONDECYT N0 1140238 (NC). GP acknowledges partial funding by Dicyt-Usach Grant N0 041531PA.

Appendix A: Maclaurin coefficients of $\phi$ and $V(\phi)$

We will find the coefficients of Maclaurin series of the functions $\phi(u)$ and $V(\phi)$ and we will compute their convergence radius. These formal expansions up to order 14th are required in Figs. 2 and 3. We will firstly find the Maclaurin series of $\phi(u)$, and then we will invert this series to obtain the expansion for $u(\phi)$. These expressions will allow us to obtain the Maclaurin series of $V(\phi)$. In fact, starting with the definitions of the coefficients $\alpha_{n}$, $\beta_{n}$, and $c_{n}$ through the formal expansions

[TABLE]

we will proceed finding the $\beta_{n}$ coefficients in three steps:

• •

First, we find the coefficients $\alpha_{n}$ of the Maclaurin series of $\phi(u)$, using its exact analytical expression.

• •

Second, we find the coefficients $c_{n}$ of the Maclaurin series of $V(u)$ by using Eq.(23).

• •

third, we find the coefficients $\beta_{n}$ of the Maclaurin series of $V(\phi)$ solving the implicit expression:

[TABLE]

where $u(\phi)$ should be computed by inverting the expansion for $\phi(u)$ (see Eq. (A1)). We begin by finding $\alpha_{n}$. To find the Maclaurin series of $\phi(u)$ we integrate Eq.(22), which gives

[TABLE]

with $s_{1}$ being an integration constant. Now we need the Maclaurin series of the integrand ${\sqrt{\cosh(u)}\over\cosh(u)-c}$, which can be derived by using the identity

[TABLE]

where:

$U_{k}(u)=\sum_{i=1}^{k}\dfrac{(-1)^{i+1}}{(i-1)!}y^{i-1}(u)\left[\left(y^{k+1-i}(u)\right)^{(n)}(u)\right]\times$ 

$\left.\quad\quad\quad\ \right.\left[\Pi_{j=1}^{i-1}(k+1-j)\right]$, and $f(u)=(F\circ y)(u)$.

Note that each $U_{k}$ also depends of $n$. If we consider $F(y)={\sqrt{y}\over y-c}$ and $\psi(u)=\cosh(u)$, we obtain that $f(u)$ which is precisely the integrand of Eq.(A3).

Now, reemplacing $y=\cosh(u)$ into Eq.(A4) and evaluating $u$ at $u=0$, we obtain

[TABLE]

where:  

$U_{k}(0)=\sum_{i=1}^{k}{(-1)^{i+1}\over(i-1)!}\left[\left(\cosh^{k+1-i}(u)\right)^{(n)}(0)\right]\times$  

$\left.\quad\quad\quad\right.\left[\prod_{j=1}^{i-1}(k-j+1)\right]$.

Using Leibniz’s rule for the derivative of a product

[TABLE]

with $A(y)=\sqrt{y}$ and $B(y)=(y-c)^{-1}$, and the expressions

[TABLE]

we obtain

[TABLE]

In order to obtain the functions $U_{k}(0)$ we use the identities

[TABLE]

From which we obtain for $l\geq 1$ and $n\geq 1$

[TABLE]

Moreover, using Eq.(A8), Eq.(A10) and Eq.(A5), we find

[TABLE]

where  $\delta_{k}=\sum_{i=0}^{k}{(2k-2i-3)!!\over(k-i)!2^{k-i}}(1-c)^{-i-1}$, which leads to the expression

[TABLE]

where $\ n\geq 1$, and

$U_{2k}(0)=\sum_{i=1}^{k}\left\{{\eta_{i,k}\over(2i-2)!}\left[\prod_{j=1}^{2i-2}(2k-j+1)\right]\right.$ 

$\left.\quad\quad\quad\quad-{\gamma_{i,k}\over(2i-1)!}\left[\prod_{j=1}^{2i-1}(2k-j+1)\right]\right\}$,  

$U_{2k-1}(0)=(2k-1)+\sum_{i=1}^{k-1}\left\{{\gamma_{i,k}\over(2i-2)!}\left[\prod_{j=1}^{2i-2}(2k-j)\right]\right.$ 

$\left.\quad\quad\quad\quad\quad-{\eta_{i+1,k}\over(2i-1)!}\left[\prod_{j=1}^{2i-1}(2k-j)\right]\right\}$, 

$\eta_{i,k}=2^{2i-2k-1}\sum_{j=0}^{k-i}{2(k+1-i)\choose j}(2(k+1-i-j))^{2n}$,  

$\gamma_{i,k}=2^{2i-2k}\sum_{j=0}^{k-i}{2k+1-2i\choose j}(2k+1-2i-2j)^{2n}$.

Finally, using the above results for $f(u)$ we obtain

[TABLE]

Now, inserting Eq.(A13) into Eq.(A3) it follows

[TABLE]

From the above expression, one can compute the first coefficients of the Maclaurin series of $\phi(u)$

[TABLE]

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 closest singularity and that $0<c<1$, we conclude that the convergence radius of Maclaurin series of $\phi(u)$ of Eq.(A14) is equal to $\arccos(c)$.

Now, we will compute the coefficients $\beta_{n}$ of the Maclaurin series of $V(\phi)$. To this aim, we use the coefficients $c_{n}$, which can be obtained in a similar way to how the coefficients $\alpha_{n}$ were computed since we can rewrite Eq.(23) as

[TABLE]

A straightforward computation yields

[TABLE]

where:

$Y^{(2n)}={2c\over 5}\sum_{k=1}^{n-1}{2n\choose 2k}g_{2(n-k)}g_{2k}+\left(1+{4c\over 5(1-c)}\right)g_{2n}$, with $Y^{(2n-1)}=0$ for $n\geq 1$, and constants $g_{2n}$ given by

[TABLE]

with

$U_{2k}=\sum_{i=1}^{k}\left\{{\eta_{i,k}\over(2i-2)!}\left[\prod_{j=1}^{2i-2}(2k-j+1)\right]\right.$

$\left.\quad\quad\quad-{\gamma_{i,k}\over(2i-1)!}\left[\prod_{j=1}^{2i-1}(2k-j+1)\right]\right\}$,

$U_{2k-1}=(2k-1)+\sum_{i=1}^{k-1}\left\{{\gamma_{i,k}\over(2i-2)!}\left[\prod_{j=1}^{2i-2}(2k-j)\right]\right.$,

$\left.\quad\quad\quad\quad-{\eta_{i+1,k}\over(2i-1)!}\left[\prod_{j=1}^{2i-1}(2k-j)\right]\right\}$,

$\eta_{i,k}=2^{2i-2k-1}\sum_{j=0}^{k-i}{2(k+1-i)\choose j}(2(k+1-i-j))^{2n}$, $\gamma_{i,k}=2^{2i-2k}\sum_{j=0}^{k-i}{2k+1-2i\choose j}(2k+1-2i-2j)^{2n}$,

$\delta_{k}=(1-c)^{-1-k}$.

The expansion of Eq.(A17) has a convergence radius $r_{c}=\arccos(c)$. Now we are equipped with the required relations to finally obtain the coefficients $\beta_{n}$ of the expansion of $V(\phi)$. We use Eq.(A2) together with the following identity

[TABLE]

where $c_{(0,n)}=a_{0}^{n}$, and $c_{(m,n)}={1\over ma_{0}}\sum_{k=1}^{m}(kn-m+k)a_{k}c_{(m-k,n)}$ for $m\geq 1$. For the particular case $s_{1}=0$, which leads to a well defined and unique solution for the parameters $\beta_{n}$, we obtain

[TABLE]

To obtain the coefficient $\beta_{k}$, we used that $\phi(u)$ is an odd function while $V(u)$ is an even function of the argument. Thus solving Eq.(A20) for $\beta_{k}$ we obtain

[TABLE]

The convergence radius of $V(\phi)$ is $r_{c}=\phi(\arccos(c))$. Finally, we evaluate explicitly the first $\beta_{n}$-coefficients from the relations of Eq.(A21)

[TABLE]

**Appendix B: Differential equations for $a$ and $\rho$ up to first order **

We will compute the coefficients $A$, $B$, $C$, $D$, $E$ and $F$ that appear in Eq.(34). We first compute the function $p_{(1)}$ appearing in Eqs.(30) and(31):

[TABLE]

which corresponds to Eq.(36). Now, in order to obtain the coefficients A, B, C, D, E and F, we insert Eqs.(31) and (32) into Eq.(2), which leads to the following expression

[TABLE]

where the “prime” ( {}_{{}^{\prime}}\) denotes the derivative with respect to $u$. Now, inserting $a_{0}(u)$, $\rho_{0}(u)$ and $a^{\prime}_{0}(u)$ in the above equation we obtain

[TABLE]

where

[TABLE]

Moreover, if we use Eqs.(31) and (33), we obtain

[TABLE]

When the functions $a_{0}(u)$, $\rho_{0}(u)$, $p_{0}(u)$ and $a^{\prime}_{0}(u)$ are inserted into Eq.(B5), $\rho^{\prime}_{1}(u)$ can explicitly be expressed as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]