Stability of a viable non-minimal bounce

The main difficulties in constructing a viable early Universe bouncing model are: to bypass the observational and theoretical \emph{no-go} theorem, to construct a stable non-singular bouncing phase and perhaps, the major concern of it is to construct a stable attractor solution which can evade the BKL instability as well. In this article, in the homogeneous and isotropic background, we extensively study the stability analysis of the recently announced viable non-minimal bouncing theory in the presence of an additional barotropic fluid and show that, the bouncing solution remains stable and can evade BKL instability for a wide range of the model parameter. We provide the expressions that explain the behavior of the Universe in the vicinity of the required fixed point i.e., the bouncing solution and compare our results with the minimal theory and show that ekpyrosis is the most stable solution in any scenario.


I. INTRODUCTION
In solving non-linear differential equations, initial conditions play a crucial role in determining the dynamical behavior of the system as different initial conditions may lead to different trajectories which may behave completely differently. Since gravity is highly nonlinear in nature, solving early Universe evolution depends highly on the initial conditions as well. However, the already highly acceptable paradigm of the early Universe -the inflationary paradigm [1][2][3][4][5][6][7][8][9][10][11][12][13], perhaps succeeds based on the fact that, not only it is in line with the recent tighter observational constraints [14,15] but it also provides an attractor solution [16,17]. This implies that the early Universe inflationary solution does not depend on the initial condition and because of this reason, even if the Universe contains additional matter(s), inflation makes them decays exponentially and the Universe solely dominates by the inflaton field responsible for inflation.
However, the greatest achievement comes from the fact that the inflationary solution can also evade Belinsky-Khalatnikov-Lifshitz (BKL) instability: the anisotropic energy density can grow faster than the energy density responsible for the evolution of the Universe [18] which can cause highly unstable system. However, even with this much of success, inflationary paradigm also suffers many crucial difficulties as well. First, despite the ever-tightening observational constraints, there seems to exist many inflationary models that continue remain consistent with the data [19][20][21][22], even leading to the concern whether inflation can be falsified at all [23]. Also, inflation being insensitive to initial conditions is still debatable.
For instance, in Ref. [24], using non-perturbative simlulations, authors showed that the inflationary expansion starts under very specific circumstances. In Ref. [25], authors pointed out the fact that small field potential fails to start the inflationary expansion under most initial conditions. Also, it suffers from the trans-Planckian problem: the cosmological scales that we nowadays observe may correspond to length scales smaller than the Planck length at the onset of inflation, which in contradiction to the inherent assumption that the scales are classical even at the initial stage of inflation [26].
These issues leads to in search for alternatives to inflationary paradigm and the most popular one is the classical non-singular bouncing scenario where, the Universe undergoes a phase of contraction until the scale factor reaches a minimum value, before it enters the expanding phase [27][28][29][30][31][32]. However, the problem with the most bouncing models is that the models are extremely difficult to construct. While, only known attractors are the ekpyrotic models [33], in fact, in Ref. [34], authors showed that the ekpyrotic contraction is a 'super-smoother', i.e., it is robust to a very wide range of initial conditions and avoid Kasner/mixmaster chaos and a similar statement could never be proven about any other primordial scenario, it fails to be in line with the observations. On the other hand, while matter bounce models may provide correct theoretical predictions consistent with the observations, the solution fails to be stable and thus, any tiny initial deviation can grow exponentially to cause an instability in the system. Another difficulty is, while constructing the contracting phase is rather easy to achieve in, obtaining the non-singular bouncing phase is extremely difficult as it requires to violate the null energy condition: often called as the theoretical no-go theorem [35][36][37][38][39][40][41][42][43][44]. And, most importantly, even if one may construct a model evading the stability and theoretical no-go theorem, perturbations suffer from many difficulties such as gradient instability, and these models fail to be in line with the observational constraints: a small tensor-to-scalar ratio (r 0.002 0.06) and simultaneously, very small scalar non-Gaussianity parameter (f NL ∼ O(1)): referred to as the observational no-go theorem [45,46]. Apart from these main three issues, in general, bouncing models possess another, rather weak, difficulty as well: a natural exit mechanism from the bouncing phase to enter into the conventional reheating phase.
In Refs. [55,56], we have shown that, the issue of stability can easily be resolved by using a simple non-minimal coupling. This has been achieved with the help of conformal transformation and we also have shown that, it may lead to viable tensor spectra as well [57]. Recently, we extended the work to construct the first viable non-minimal bouncing model which evades all the issues and is consistent with the constraints [54]. This model is constructed in such a way that the scalar sector of the action is conformal to that of the slow-roll inflationary model. Since, under conformal transformation, perturbations remain invariant and the constraints are obtained from the correlation functions of the scalar and tensor perturbations, we achieved to bypass observational no-go theorem. The choice of the coupling function also helps us to achieve the non-singular bouncing phase without any in-stability and therefore, it evades the theoretical no-go theorem and leads to the conventional reheating scenario. Also, it has been pointed out that since the work is similar to works in Refs. [55,56], it can also resolve the issue of stability, mainly the BKL instability as well. However, the model contains a free model parameter α (which needs to be greater than zero to achieve the bouncing solution) and there is no bound on it. Also, it is not clear how the system behaves in the presence of an additional matter.
At this moment, we need to stress that, when a model is studied for perturbations and is compared with the observations, it is assumed that that, either the initial conditions are chosen with extreme measure, or the model solution is stable enough so that there is no need to fine-tune for choosing the initial conditions. Otherwise, even without the external matter, the perturbations itself can cause the instability as it also possesses a tiny energy. This is the main reason why stable is solution is always preferred (one may even argue that only the stable solutions are allowed). Therefore, even our model is conformal to inflationary solutions and obeys all relevant observational constraints, one needs to verify the stability of the model, especially in the presence of an external matter. This leads to the main motivation of this article.
In this work, we analyze the stability of the model in presence of an additional barotropic matter. The aim is to find the general attractor behavior of the bouncing solution that is consistent with the observations. There is mainly one free model parameter α (for minimal theory, α = −1, since for this value, the coupling function becomes unity) along with the equation of state of the additional barotropic fluid w m . In addition to that, since we consider conformal single field model that leads to slow-roll inflation, we also consider the potential slow-roll parameter in the minimal Einstein frame to be (nearly) constant (in case of exponential potential is strictly constant, as examined in Refs. [16,17]). With the help of these model parameter, we find the condition for stability for the bouncing solution and show that, to avoid the BKL instability leads to α < 2. Since bouncing solution requires α > 0, the bound becomes 0 < α < 2. Also, in order to understand the behavior of the stability, we study the phase space in the vicinity of the fixed point corresponding to the bouncing solution and find that the system acts like the Universe contains three different types of matter and due to the attractor behavior only the leading bouncing matter component remains intact, whereas, other two decay very fast. This helps to realize how the external barotropic fluid influence the system and how the attracting nature of the bouncing solution get rid of the initial deviations.
A few words on our conventions and notations are in order at this stage of our discussion.
In this work, we work with the natural units such that = c = 1, and we define the Planck mass to be M Pl = (8 π G) −1/2 . We adopt the metric signature of (−, +, +, +). Also, we should mention that, while the Greek indices are contracted with the metric tensor g µν , the Latin indices contracted with the Kronecker delta δ ij . Moreover, we shall denote the partial and the covariant derivatives as ∂ and ∇. The overdots and overprimes, as usual, denote derivatives with respect to the cosmic time t and the conformal time η associated with the Friedmann-Lemaître-Robertson-Walker (FLRW) line-element, respectively. The sub(super)script 'I' and b denote the quantity in the minimal Einstein theory and nonminimal theory, respectively.

II. A BRIEF INTRODUCTION OF THE BOUNCING MODEL
The non-minimal bouncing model is constructed in such a way that the scalar sector of the non-minimal action transforms conformally to that of a minimal canonical Einstein model that leads to slow-roll inflation. The Slow-roll inflationary model can easily be constructed by using a single canonical scalar field φ minimally coupled to the gravity, i.e., the Einstein gravity as R I is the Ricci scalar for the metric g I µν , and V I (φ) is the potential responsible for the inflationary solution. The part of the action S I M (g µν , Ψ M ) indicates the presence of an additional fluid. As mentioned above, since inflation is an attractor, the evolution of the field do not depend on the additional fluid and solely governed by the scalar sector of the system. Assuming the above theory is responsible for the slow-roll inflationary dynamics, the scale factor solution, during inflation, can approximately be written as a function of the scalar field φ as with V I,φ ≡ ∂V I ∂φ . Now we shall construct a model which is conformal to the scalar part of the above inflationary action (1) in such a way that the new scale factor behaves as a bouncing solution. The transformation can be written as g b µν and a b (η) are the required bouncing metric and scale factor solutions, respectively, along with f (φ) being the coupling function. Using the above conformal transformation along with the action (1), we can construct the action responsible for such bouncing solution as ω(φ) and the potential V b (φ) depend on the coupling function f (φ) and the inflationary where f ,φ ≡ ∂f ∂φ . S b M is the external fluid present in the system in the non-minimal theory. Note that, while the scalar part of the above action (4) is conformal to that of (1), these two actions in reality are not conformal due to the presence of the additional fluid. However, if the scalar field dominated solution in the non-minimal theory also becomes stable, similar to minimal theory, then, for a given inflationary model (1) with its solution (2), one can obtain the bouncing scale factor solution a b (η) in the non-minimal theory by setting f (φ) in a certain manner. For simplicity, if we desire the bouncing solution of the form where η is the comoving time, then by using Eq. (2) and (3), one can obtain the required solution of the coupling function f (φ) as f 0 is normalized in such a way that at the minima of the potential, f (φ) becomes unity.
The non-minimal theory (4) is the desired bouncing model: given an inflationary potential V I (φ), we can completely determine the coupling function f (φ), which in turn, provides non-minimal theory (4) responsible for the bouncing scale factor solution (6). Note that, since Eq. (2) is an approximated form, by using the above expressions (7) and (3), the obtained scale factor solution a b (η) is also not exact, but close to (6) (we will evaluate it later). Also, we need to stress that the form of the scale factor (6) is valid only in the contracting phase, and during and after the bounce, it is extremely difficult to obtain analytical solution. Therefore, the bouncing model, in our case, is also asymmetrical. Also, the underlying assumption of the theory is that, either the bouncing solution is stable or we must choose precise initial condition that leads to the bounce. Last of all, while α > 0 provides bouncing solution, there are actually no bound on α as it can take any value and negative values of it will lead to different Universe with different scale factor solution (6).
For instance, α = −1 brings us the conventional minimal Einstein scenario as, for this value, the coupling function becomes unity. In this work, we will concentrate only on the bouncing solutions where α > 0 and obtain the condition for stability and the behavior of the system in the vicinity of the contracting solution, which we will explore in the next section.
Few things to note before we move in the next section. In case of scalar perturbation, in the Einstein frame (1), the perturbed action is where, ζ is the curvature perturbation and ∇ is the divergence operator. In the non-minimal frame (4), the expression of the action for the scalar perturbation changes by replacing z I (η) with z b (η). However, z b (η) is simply the conformally transformed z I (η), i.e., . This implies that, z b (η) = z I (η) which leads to the action being identical in both frames and hence, ζ I = ζ b at linear order. Also, as the speed of sound is conformally invariant, c s is unity in both minimal inflationary as well as in non-minimal bouncing scenarios. Similarly, one can evaluate the perturbed interaction Hamiltonian (for detailed evaluation, see Refs. [58][59][60][61][62]) at any order and show that, it also remains invariant in all conformal frames.
This implies that the curvature perturbation at any order in both the frames are identical, which is not surprising, as we know, under conformal transformation, curvature perturbation remains invariant. The argument extends to tensor perturbation as well. Therefore, given a viable inflationary model which is consistent with the observations, in our bouncing model, there appear no divergences or instabilities (e.g., gradient instability) in the perturbations anywhere, even at the bounce and, in addition to that, the model also satisfies all observational constraints: thus evading the no-go theorem. Also, the model leads to stable non-singular bounce as well as, shortly after the bounce, it leads to conventional reheating phase (for detailed information, see Fig. 1 and Ref. [54] as well as [63], to compare reheating in minimal theory with the same in non-minimal theory).

III. STABILITY ANALYSIS
Let us now analyze the stability condition for the above system (4). Assuming the additional fluid is barotropic, in nature, with the equation of state w m , the governing equations become where, T (M ) µν is the energy-momentum tensor of the additional barotropic fluid in the nonminimal theory (4), satisfying We have not used any subscript b since, from now on, all concerning equations will correspond only to the non-minimal theory. These equations, using Friedmann-Robertson-Walker (FRW) homogeneous and isotropic Universe with the line element turn into where ρ M ≡ −T (M )0 0 and T i j ≡ P M δ i j = w m ρ M δ i j are the energy density and the pressure of the barotropic fluid with w m being the equation of state. To study the stability analysis, instead of using quantities with dimensions, we can express and re-write these field equations in terms of dimensionless quantities. These are defined as Using these dimensionless variables, one can quickly express the evolution equation of them along with the constraint equation N is defined as the e-folding number in the non-minimal frame ≡ ln (a/a 0 ) , a 0 is the initial value of the scale factor and V I (φ) is the potential in the minimal Einstein frame and is related to V (φ) through Eq. (5). Note that, the model parameter γ e actually represents a parameter in the Einstein frame which is completely arbitrary and can be fixed from the observations. For example, while near the pivot scale, for chaotic inflation, γ e ∼ 0.01, for Starobinsky model, at the same scale, it becomes ∼ 0.003. However, chaotic inflation has already been ruled out from the observations while, the Starobinsky model still remains consistent with it and one of the most important model in the paradigm. Also, we should stress that the phase during the stability is analyzed is the contraction in a bouncing model and therefore, instead of e-N-fold time convention, here we are using e-fold time convention.
One can also express the slow-roll parameter b in terms of these dimensionless parameter Using the above expression, one can define the effective equation of state as Although γ e is a function of φ and thus, it is dynamical, the minimal theory that we are considering leads to slow-roll inflation and thus one can consider it to be nearly constant, similar to exponential potential.
In this system, as one can see, there exist three model parameters: {α, w m , γ e }. w m is arbitrary as it describes the nature of the additional field. For example, in case of anisotropic stress fluid that corresponds to the stiff matter, w m = 1. This field results in the BKL instability that, later, we will investigate. Also, as mentioned above, γ e is nearly constant and hence, we are left with one model parameter α that can be constrained. In the next section, we will try to constrain it using the lyapunov exponents.

A. Fixed points and the bouncing solution
Using the evolution equation, one can find the fixed/critical point of the system and these points correspond to the solutions of the system, in this case, the Universe depicted by the non-minimal theory (4). These can be found by setting Eqs. (18) and (19) to be equal to zero, i.e., the velocities of x and y vanishes at these points. There are seven such fixed points: 2.
The sign of the y * signify whether the Universe is expanding or contracting. Since, H appears in the denominator in the expression of y in (17), negative sign of y leads to contraction (bouncing), while positive y corresponds to expanding Universe. Therefore, {x * 1 , y * 1 } leads to contraction and {x * 2 , y * 2 } implies the expanding Universe. In fact, {x * 1 , y * 1 } is our desired solution and the scalar field dominated solution as the fractional energy density of the additional fluid and the effective equation of state (see Eq. (22)) of corresponds to this point are which corresponds to the the scale factor solution as As you can see, the exponent is not exactly α but β, which is contrary to the equation (6).
However, as mentioned before, γ e 1 for slow-roll models and only by using such model, we achieved such non-minimal bounce, and therefore, β α. Therefore, {x * 1 , y * 1 } corresponds to our bouncing solution that we will investigate thoroughly in the next section. Other fixed points are 6.
Since our desired fixed point is {x * 1 , y * 1 }, in this work, we will focus only in this point. In the next subsection, we will explore the Lyapunov exponents that determine the stability of the solution and then later, we will focus on the behavior of the Universe close to the fixed point.

B. Lyapunov exponent and the stability of the bouncing solution
Now, in order to study the stability of these fixed points, we need to linearize the equations (18) and (19) as where, A(x, y) and B(x, y) are the right-hand side of (18)  In other words, both eigenvalues being negative (positive) in an expanding (contracting) Universe implies that the fixed point is stable and the corresponding solution is an attractor solution.
Let us now calculate the Lyapunov exponents for our desired fixed point (23). It takes the form The above expression bears one of the main result of this work. As one can see properly, since the desired fixed point is the bouncing solution, we require both of the λ's to be positive.
Few things to note: since α = −1 represents the minimal solution as mentioned before (coupling function (7) becomes unity), one can easily obtain the results in Refs. [16,17] where a is the scale factor. Therefore, since BKL instability is caused by it, the condition to evade the BKL instability can easily be obtained by substituting w m = 1: Therefore, for exponential potential, γ e is constant and positive domain of α can only be obtained if γ e < 2. In case of slow-roll, γ e 1 and therefore, for slow-roll inflationary models which are conformal to the non-minimal bounce, the condition becomes as bouncing solution can only be obtained for α > 0. This is the main result of this article as, now, α cannot possess arbitrary positive values for stable viable solution. In the next section, in order to understand how the deviations impact the system near the fixed point, we will study the phase phase in the concerned region.

IV. EVOLUTION OF THE DEVIATIONS
Using the eigenvalues and the eigenvectors, one can easily solve (33). For the fixed point (23), it becomes δy (1) where, λ 1 and λ 2 are given in (34). The A's can be evaluated as where, Once we obtain the solution expression for the slow-roll parameter, we can solve the energy density equation as (N ) = −d N ln H(N ) and by performing the integral, one can easily solve for H. It becomes, H 2 (N ) = H 2 0 exp −2 dN (N ) , H 0 is the initial value of the Hubble parameter. Using (37) and (43), and by using the fact that, δx 0 and δy 0 are tiny and act like perturbations, we obtain the evolution of the Hubble parameter as where, One can easily obtain the the above expressions in terms of model parameters {α, w m , γ e } by using relations (39), (40), (41), (42), (44) and (45). It is clear from the above equation From the above expression, it becomes obvious that, for bouncing solution, if λ's are positive, the energy fraction vanishes after a while. In order to understand the correct nature of it, below, we presented three different scenarios with different choices of model parameters that try to help us realize the system.

A. Radiation bounce
As we know from (36), α must be positive and less than 2 in order to satisfy bouncing solution as well as to evade the BKL instability. Therefore, in the first example, we consider radiation bounce where α = 1, i.e., the scale factor solution approximately becomes a b (η) ∝ (−η). Also, we consider the additional fluid as the anisotropic stress fluid with stiff equation of state w m = 1. Assuming typical conformal slow-roll model (e.g., chaotic inflation with with γ e ≈ 0.01, the Hubble solution takes the form with Ω m1 Ω m2 = −1.125 δx 0 + 0.002 δy 0 δx 0 + 0.002 δy 0 .
The equation represents the attractor nature of the solution as the leading order solution, i.e., the effective radiation matter solution dominates over the additional fluid evolution.
Also, the actual energy density fraction of the additional fluid (50) takes the form For bouncing solution, as it is obvious, it quickly vanishes after a while.

B. Comparing with the minimal scenario
As it has been mentioned before, negative α represents expanding solution and α = −1 reduces to the minimal Einstein inflationary scenario as for this value, the coupling function becomes unity. In this subsection, we will compare the above result with the minimal inflationary case. Again using the similar value of the additional fluid equation of state as well the γ e , i.e., w m = 1, γ e ≈ 0.01, the Hubble solution becomes with Ω m1 Ω m2 = −0.05 + 122.47 δy 0 δx 0 .
Notice that, in this case, the additional fluids with Ω m1 and Ω m2 initial mass fractions decay faster compared to the radiation bounce scenario as the Lyapunov exponents, i.e., the λ's are having high values. This can even be seen from the expression of energy fraction of the additional fluid (50): It vanishes faster than that of the radiation bounce case. Therefore, one may consider, by judging the capability of stability solution, minimal inflationary models are preferred compared to the non-minimal bouncing models. In the next subsection, we will show that, it is clearly not the case.

C. Comparing with the viable ekpyrosis
In Ref. [34], it has been studied extensively and showed that the ekpyrosis is a superattractor (although it has been studied for minimal models, one can easily extend it for non-minimal solution as well. See Ref. [55] in this context). It can also been seen from (34) as with Ω m1 Ω m2 = −0.95 δx 0 − 0.003 δy 0 δx 0 + 0.003 δy 0 .
Again, compared to the leading solution, the externals fluids decay much faster than the minimal case, shown above. This is because of the Lyapunov exponents are, in this case, 38 and 30, respectively. For instance, the energy fraction of the additional fluid becomes We also showed that the ekpurosis is the best attarctor among them and thus is always preferred.
While we studied the system for homogeneous and isotropic background, one can easily extend the work for homogeneous and anisotropic system as well and may arrive at the similar condition. One can even go beyond the homogeneous system and study the stability using numerical relayivity. However, this is beyond the scope of the this article. Also, we did not consider the negative values of α as well as α = 0 case. The last one signifies the emergent gravity scenario. Apart from that, we also did not study the behavior of the other fixed points (28), (29), (30), (31) and (32). While in the minimal scenario, most of them are unstable, in our case, it may not be the same and may impact the system heavily. We reserve these works for our future prospect.