Global phase space analysis for a class of single scalar field bouncing solutions in general relativity

Abstract

We carry out a compact phase space analysis of a non-canonical scalar field theory whose Lagrangian is of the form $F(X)-V(\phi )$ within general relativity. In particular, we focus on a kinetic term of the form $F(X)=\beta X^m$ ($m\ne 1/2$) with power-law potential $V_0 \phi ^n$ and exponential potential $V_0 e^{-\lambda \phi /M_{Pl}}$ of the scalar field. The Cuscuton case $m=1/2$ where the scalar field is non-dynamical is left out of consideration. The main aim of this work is to investigate the genericity of nonsingular bounce in these models and to investigate the cosmic future of the bouncing cosmologies when they are generic. A global dynamical system formulation that is particularly suitable for investigating nonsingular bouncing cosmologies is used to carry out the analysis. We show that when $F(X)=\beta X^m$ ($\beta <0$), nonsingular bounce is generic for a power law potential $V(\phi ) = V_0 \phi ^n$ only within the parameter range $\left\{ \frac{1}{2}<m<1,\,n<\frac{2\,m}{m-1}\right\}$ and for an exponential potential $V(\phi ) = V_0 e^{-\lambda \phi /M_{Pl}}$ only within the parameter range $\left\{ \frac{1}{2}<m\le 1\right\}$. Except in these cases, nonsingular bounce in these models is not generic due to the non-existence of global past or future attractors. Our analysis serves to show the importance of a global phase space analysis to address important questions about nonsingular bouncing solutions, an idea that may and must be adopted for such solutions even in other theories.

1 Introduction

The nonsingular bouncing paradigm has been under study for quite a long time now as an alternative to the inflationary paradigm. The main problem with the inflationary paradigm is that it is inherently singular; timelike and null geodesics cannot be extended indefinitely towards the past in a nonsingular way [1]. Studies on the nonsingular bouncing paradigm were fuelled mainly by the theoretical necessity of avoiding this singularity problem. Nonsingular bouncing cosmologies consist of a pre-bounce contracting phase ($\dot{a}(t)<0$, a(t) being the scale factor of the Universe) connecting smoothly through a nonsingular bounce ($\dot{a}=0,\,\ddot{a}>0$ or $H=0,\,\dot{H}>0$, H(t) being the Hubble parameter) to the post bounce expanding phase ($\dot{a}(t)>0$) that we live in. Various aspects of the nonsingular bouncing solutions have been studied in the literature, for example, the anisotropy problem during the contraction phase and the evolution of cosmological perturbations through a bounce. For a critical review of the bouncing paradigm, one can see Refs. [2, 3]. An aspect that is rarely addressed in relevant literature is, even if a theory admits nonsingular bouncing solutions, how generic these bouncing solutions are i.e. how sensitive are these solutions to small changes in the initial conditions of the universe? This is precisely the question that we will address in this article.

In this article, we will be dealing with a single scalar field model whose Lagrangian density can be written as $\mathcal {L}=F(X)-V(\phi )$. With a phantom field, this model can produce nonsingular bouncing solutions [4, 5]. It is known that if one tries to achieve nonsingular bouncing solutions within general relativity (GR), one needs to allow for the violation of the null energy condition (NEC) and/or the strong energy condition (SEC) [6]. NEC violation is necessary for the spatially flat case we consider here. When the dynamics of inhomogeneous cosmological perturbations are considered, NEC-violating scalar fields may exhibit ghost instability [7,8,9] and gradient instability [10]. We must mention here that attempts to construct stable nonsingular bouncing models despite the fact that NEC violation has given birth to more sophisticated and involved models that involve, instead of a usual scalar field, either ghost condensates [11, 12] and Galileons [13,14,15]. However, since the main focus of this paper is not an analysis of the dynamics of inhomogeneous perturbations but to address the question of stability of a bouncing solution with slight perturbation in the initial condition, we stick to the more straightforward $F(X)-V(\phi )$ model here to explain our method. The analysis can, in principle, of course, be generalized to the more sophisticated bouncing models mentioned above.

Considering a simple form of the kinetic term F(X) and two particular examples for the potential $V(\phi )$, we will attempt to address the question of the generality of nonsingular bouncing solutions. The canonical field, often referred to as the quintessence field with potential, is the most basic type of scalar field [16]. However, there are several complex cosmological dynamics in the universe that the canonical scalar field is unable to fully explain. For instance, the crossing of the phantom divide line and the bouncing solution cannot be explained by the quintessence field model. This leads to a broader explanation of the non-canonical scalar field, a type of scalar field. The coincidence problem may be solved in the non-canonical scenario without causing any fine-tuning problems, which is another benefit over the canonical setup. In addition, the tensor-to-scalar ratio in non-canonical models is smaller than in canonical models, resulting in better agreement with CMB measurements. We are inspired to learn more about the cosmic dynamics of non-canonical scalar fields by their intriguing properties [17].

Our approach to tackling the problem of genericity of bouncing solutions will be to carry out a dynamical system analysis of the models and investigate the nonsingular bouncing solutions within the phase space. The dynamical systems technique has been a very useful tool in understanding the qualitative behaviors of cosmological models without analytically solving the system of differential equations [18, 19]. With this technique, one can reframe cosmological dynamics as the phase flow in a suitably defined phase space. Although this technique is now extensively used in investigating inflationary and dark energy models [20], it is rarely used to investigate nonsingular bouncing cosmologies. The main reason seems to be that the Hubble normalized dimensionless dynamical variables that are mostly used for the cosmological phase space analysis diverge at a nonsingular bounce. On the other hand, the dynamical system formulation seems to be perfectly fitted to answer qualitative questions like the genericity of a cosmological solution. By genericity of the bouncing solutions in theory, we mean how sensitive they are to the initial conditions. Suppose one sets some initial conditions at a random moment during the contracting phase of the evolution, numerically evolves the field equations and obtains a bounce. The question is whether one would also get a bounce if (s)he applies a random slight change in the initial condition. To investigate nonsingular bouncing solutions in the phase space picture, one needs to either define an alternative set of dimensionless dynamical variables or compactify the infinite phase space into a finite domain using a specific compactification prescription. Apart from the question of genericity, the phase space picture also provides information about the past and future asymptotics of a nonsingular bouncing cosmology. The future asymptotic is really important, as they determine the end states of the bouncing cosmologies. Do they end up being asymptotically De-Sitter like in $\Lambda$CDM, or lead to a big-rip? These are our main motivations behind this work.

Cosmological phase space of $F(X)-V(\phi )$ models have been investigated in Ref. [21], where the authors also attempt to touch upon the bouncing solutions from a phase space point of view by defining an alternative set of dynamical variables. Ref. [22] extends the same analysis to the case of Bianchi-I. However, none of the above works carries out a compact phase space analysis. The system, in general, may contain interesting cosmological scenarios described by fixed points hidden at infinity; in that case, compact analysis is required to extract the global dynamics of the system. For example, in Ref. [21], even though the authors could show phase trajectories that represent nonsingular bounce, one cannot really answer the question of genericity and future asymptotics. In the present work, we employ the same dynamical system formulation as was used in Refs. [21, 22] but complement these earlier works by carrying out a compact phase space analysis and hence answering the aforementioned questions. For the functions F(X) and $V(\phi )$ considered in this analysis, we explicitly identify the scenarios where a bounce occurs generically. Additionally, we determine when these scenarios exhibit asymptotically De-Sitter behavior or lead to a big-rip. It is crucial to note that the absence of generic bouncing solutions does not imply that bouncing solutions are unattainable altogether. Under certain specific initial conditions, it is still possible to obtain a bounce. However, unlike the case of generic bounces, a chosen initial condition that yields a bounce cannot be assured to retain this property under a random slight perturbation. Consequently, the cases where the bounce is generic are particularly intriguing when constructing bouncing models.

The paper is structured as follows: In Sect. 2, we present the basic equations for $F(X)-V(\phi )$ model considering a spatially flat homogeneous and isotropic Friedmann-Lamaitre-Robertson-Walker (FLRW) Universe. In Sect. 3, we present the dynamical system formulation suitable for investigating nonsingular bouncing cosmologies and carry out the phase space analysis for two particular models. In Sect. 4, we summarize the take-home results from our phase space analysis. Finally, we conclude in Sect. 5.

2 Basic cosmological equations for $\mathcal {L}=F(X)-V(\phi )$

The most general action of a minimally coupled scalar field theory is given by

$$\begin{aligned} S=\int d^4x \sqrt{-g}\left( \frac{M_{Pl}^2}{2} R+\mathcal {L}(\phi ,X)\right) +S_m, \end{aligned}$$

(1)

where $M_{Pl}$ is the reduced Planck mass, R is the Ricci scalar, g is the metric determinant, $\mathcal {L}(\phi , X)$ is the Lagrangian density of scalar field $\phi$ whose kinetic component is denoted by X (i.e., $X = -\frac{1}{2}\partial _{\mu } \phi \partial ^{\mu } \phi$) and the last term $S_m$ is the action corresponds to the matter component taken to be a perfect fluid.

Variation of (1) with respect to the metric $g_{\mu \nu }$ yields the Einstein field equations given by

$$\begin{aligned} G_{\mu \nu }= T_{\mu \nu }^{(\phi )} + T_{\mu \nu }^{(m)}, \end{aligned}$$

(2)

where $G_{\mu \nu }$ denotes the Einstein tensor, $T_{\mu \nu }^{(\phi )}$ is the scalar field energy-momentum tensor given by

$$\begin{aligned} T_{\mu \nu }^{(\phi )}=\frac{\partial \mathcal {L}}{\partial X} \partial _{\mu } \phi \partial _{\nu } \phi -g_{\mu \nu } \mathcal {L}, \end{aligned}$$

(3)

and the matter energy-momentum tensor $T_{\mu \nu }^{(m)}$ is given by

$$\begin{aligned} T_{\mu \nu }^{(m)}=(\rho _m+P_m) u_\mu u_\nu +P_m g_{\mu \nu }. \end{aligned}$$

(4)

Here $\rho _{m}$ and $P_m$ denote the energy density and pressure of the matter component, respectively, with the four-velocity vector $u_\mu$. In this paper, we consider a spatially flat FLRW cosmology described by the metric

$$\begin{aligned} ds^2=-dt^2+a(t)^2(dx^2+dy^2+dz^2), \end{aligned}$$

(5)

where a(t) is a scale factor, t is cosmic time and x, y, z are the Cartesian coordinates. We also focus on a scalar field model whose generic form of the Lagrangian is given by

$$\begin{aligned} \mathcal {L}(\phi ,X) = F(X) - V(\phi ), \end{aligned}$$

(6)

where $V(\phi )$ is a scalar field potential and F(X) is an arbitrary function of X.

The energy–momentum tensor (4) of a perfect fluid under the FLRW cosmology is

$$\begin{aligned} T^{\mu (m)}_{\nu }= \text {diag}(-\rho _m, P_m, P_m, P_m) \end{aligned}$$

(7)

The spatially flat FLRW space time transforms the above Einstein field equations to the following cosmological field equations,

$$\begin{aligned} H^{2}= & {} \frac{1}{3M_{Pl}^{2}}\left[ 2XF_X - F + V + \rho _{m}\right] , \end{aligned}$$

(8a)

$$\begin{aligned} \dot{H}= & {} -\frac{1}{2M_{Pl}^{2}}\left[ 2XF_X + (1+\omega _{m})\rho _{m}\right] , \end{aligned}$$

(8b)

where $H=\frac{\dot{a}}{a}$ is the Hubble parameter and an over dot denote derivative with respect to t, the subscript X denotes derivative with respect to X and $\omega _{m}$ is the perfect fluid equation of state parameter defined as $P_m=\omega _{m} \rho _{m}$. When there is no energy exchange between the fluid and the field, the fluid component scales according to the ordinary continuity equation

$$\begin{aligned} \dot{\rho }_m + 3H(1+\omega _m)\rho _m = 0 \quad \Rightarrow \quad \rho _{m}\propto a^{-3(1+\omega _{m})}, \end{aligned}$$

(9)

and the scalar field satisfies the generic Klein–Gordon equation

$$\begin{aligned} \frac{d}{dN}(2XF_X - F + V) + 6XF_X = 0, \end{aligned}$$

(10)

where $N=\ln a$.

In the next section, we shall construct the corresponding autonomous system for the above cosmological equations and then we investigate the bouncing scenarios via a global dynamical system analysis.

3 Dynamical system formulation of $\mathcal {L}=F(X)-V(\phi )$ models suitable for investigating nonsingular bounces

In order to obtain an autonomous system of equations that is suitable to investigate nonsingular bouncing solutions through a phase space point of view, we follow a dynamical system construction that was presented in Ref. [21]. We define the dynamical variables

$$\begin{aligned} \begin{aligned} x&= \frac{\sqrt{3}M_{Pl}H}{\sqrt{|\rho _k|}}, \quad y = \sqrt{\frac{|V|}{|\rho _k|}}{{\,\textrm{sgn}\,}}(V), \quad \Omega _{m} = \frac{\rho _{m}}{\vert \rho _{k}\vert },\\ \sigma&= -\frac{M_{Pl}V_{\phi }}{V}\sqrt{\frac{2X}{3|\rho _k|}}{{\,\textrm{sgn}\,}}(\dot{\phi }) = -\frac{M_{Pl}}{\sqrt{3|\rho _k|}}\frac{dlog{V}}{dt} \end{aligned} \end{aligned}$$

(11)

where we have denoted the kinetic part of the energy density by

$$\begin{aligned} \rho _{k} = 2XF_{X} - F. \end{aligned}$$

(12)

The kinetic part of the pressure is just

$$\begin{aligned} P_k = F. \end{aligned}$$

(13)

Motivated with this, one can define

$$\begin{aligned} \omega _{k} \equiv \frac{P_k}{\rho _k} = \frac{F}{2XF_{X}-F}, \end{aligned}$$

(14)

which can be interpreted as the equation of state parameter for the kinetic part of the Lagrangian. The equation of state of the scalar field can be obtained in terms of the new variable as,

$$\begin{aligned} \omega _{\phi }=\frac{p_{\phi }}{\rho _{\phi }}=\frac{\omega _{k}x^{2}-y^{2}}{x^{2}+y^{2}}. \end{aligned}$$

(15)

Next, we define two auxiliary variables

$$\begin{aligned} \Xi =\frac{X F_{XX}}{F_{X}}, \quad \Gamma =\frac{V V_{\phi \phi }}{V_{\phi }^{2}}, \end{aligned}$$

(16)

which will be required to write the dynamical system. Lastly, we define the phase space time variable [21]

$$\begin{aligned} d\tau = \sqrt{\frac{|\rho _{k}|}{3M_{Pl}^{2}}}dt. \end{aligned}$$

(17)

With respect to the dynamical variables and auxiliary variables defined in Eqs. (11) and (16), the Friedmann constraint and the dynamical equations become

$$\begin{aligned}{} & {} x^2 - y|y| - \Omega _{m} = 1\times {{\,\textrm{sgn}\,}}(\rho _{k}), \end{aligned}$$

(18a)

$$\begin{aligned}{} & {} \frac{dx}{d\tau } = \frac{3}{2}x\left[ (\omega _{k}+1)x-\sigma y|y|{{\,\textrm{sgn}\,}}(\rho _{k})\right] \nonumber \\{} & {} \quad - \frac{3}{2}\left[ (\omega _{k}-\omega _{m}){{\,\textrm{sgn}\,}}(\rho _{k})+(1+\omega _{m})(x^{2}-y|y|)\right] ,\nonumber \\ \end{aligned}$$

(18b)

$$\begin{aligned}{} & {} \frac{dy}{d\tau } = \frac{3}{2}y\left[ -\sigma +(\omega _{k}+1)x-\sigma y|y|{{\,\textrm{sgn}\,}}(\rho _{k})\right] ,\end{aligned}$$

(18c)

$$\begin{aligned}{} & {} \frac{d\sigma }{d\tau } = -3\sigma ^{2}(\Gamma -1)\nonumber \\{} & {} \quad + \frac{3\sigma [2\Xi (\omega _{k}+1)+\omega _{k}-1]}{2(4\Xi +1)(\omega _{k}+1)}\left( (\omega _{k}+1)x-\sigma y^{2}\right) .\nonumber \\ \end{aligned}$$

(18d)

Because of the definition of the dynamical variable y, y|y| can also be written as $y^{2}{{\,\textrm{sgn}\,}}(V)$. Therefore, the dynamical system (18) can also be written as

$$\begin{aligned}{} & {} x^2 - y^{2}{{\,\textrm{sgn}\,}}(V) - \Omega _{m} = 1\times {{\,\textrm{sgn}\,}}(\rho _{k}),\end{aligned}$$

(19a)

$$\begin{aligned}{} & {} \frac{dx}{d\tau } = \frac{3}{2}x\left[ (\omega _{k}+1)x-\sigma y^{2}{{\,\textrm{sgn}\,}}(V){{\,\textrm{sgn}\,}}(\rho _{k})\right] \nonumber \\{} & {} \quad - \frac{3}{2}\left[ (\omega _{k}-\omega _{m}){{\,\textrm{sgn}\,}}(\rho _{k})+(1+\omega _{m})(x^{2}-y^{2}{{\,\textrm{sgn}\,}}(V))\right] ,\nonumber \\ \end{aligned}$$

(19b)

$$\begin{aligned}{} & {} \frac{dy}{d\tau } = \frac{3}{2}y\left[ -\sigma +(\omega _{k}+1)x-\sigma y^{2}{{\,\textrm{sgn}\,}}(V){{\,\textrm{sgn}\,}}(\rho _{k})\right] ,\nonumber \\ \end{aligned}$$

(19c)

$$\begin{aligned}{} & {} \frac{d\sigma }{d\tau } = -3\sigma ^{2}(\Gamma -1)\nonumber \\{} & {} \quad + \frac{3\sigma [2\Xi (\omega _{k}+1)+\omega _{k}-1]}{2(2\Xi +1)(\omega _{k}+1)}\left( (\omega _{k}+1)x-\sigma y^{2}\right) . \end{aligned}$$

(19d)

Since $\Omega _m\ge 0$, the phase space of the system (19) is given by

$$\begin{aligned} \lbrace (x,y,\sigma ) \in \mathbb {R}^{3}: x^2-y^2{{\,\textrm{sgn}\,}}(V)-{{\,\textrm{sgn}\,}}(\rho _k)\ge 0\rbrace . \end{aligned}$$

(20)

It is important to express the quantity $\frac{\dot{H}}{H^2}$ in terms of the dynamical variables as this will help us to find the cosmological evolution corresponding to a fixed point

$$\begin{aligned} \frac{\dot{H}}{H^{2}}=-\frac{3}{2x^{2}}\left[ (1+\omega _{m})(x^{2}-y\vert y\vert )+(\omega _{k}-\omega _{m}){{\,\textrm{sgn}\,}}(\rho _{k})\right] . \end{aligned}$$

(21)

Additionally, one can introduce the effective equation of state $\omega _\textrm{eff}$ given by

$$\begin{aligned} \omega _\textrm{eff}&=-1-\frac{2}{3}\frac{\dot{H}}{H^2}\nonumber \\&=-1+\frac{1}{x^{2}}\left[ (1+\omega _{m})(x^{2}-y\vert y\vert )+(\omega _{k}-\omega _{m}){{\,\textrm{sgn}\,}}(\rho _{k})\right] . \end{aligned}$$

(22)

Below we consider some specific examples of $F(X)-V(\phi )$ models and look for nonsingular bouncing solutions in the phase space.

3.1 Specific case I: $F(X)=\beta X$

In this section we consider scalar field Lagrangians of the form

$$\begin{aligned} \mathcal {L}(\phi ,X) = \beta X - V(\phi ). \end{aligned}$$

(23)

Since $\rho _{\phi }+P_{\phi }=2XF_X={2\beta X}$ and $X\ge 0$, it reduces to a canonical scalar field when $\beta >0$ and to a phantom scalar field when $\beta <0$. One can notice from Eq. (8b) that for $\dot{H}>0$ near the bounce one must necessarily require $\beta <0$, i.e. a phantom scalar field.

3.1.1 Power law potential $V(\phi )=V_0 \phi ^n$

The first example we consider is the specific case given by $F(X)=\beta X,\,V(\phi )=V_0 \phi ^n$ where $\beta , V_0$ are constants with suitable dimension and n is a dimensionless constant. For this choice, we have

$$\begin{aligned} \Xi =0,\quad \Gamma =1-\frac{1}{n}, \quad \omega _{k}=1, \quad \rho _{k}=\beta X. \end{aligned}$$

(24)

Finite fixed point analysis

In this case, the system (19) reduces to

$$\begin{aligned}{} & {} \frac{dx}{d\tau } = \frac{3}{2}x\left[ 2x - \sigma y^2 {{\,\textrm{sgn}\,}}(V){{\,\textrm{sgn}\,}}(\beta )\right] \nonumber \\{} & {} \quad - \frac{3}{2}\left[ (1-\omega _{m}){{\,\textrm{sgn}\,}}(\beta ) + (1+\omega _{m})(x^{2}-y^2 {{\,\textrm{sgn}\,}}(V))\right] , \nonumber \\ \end{aligned}$$

(25a)

$$\begin{aligned}{} & {} \frac{dy}{d\tau } = \frac{3}{2}y\left[ -\sigma + 2x - \sigma y^2 {{\,\textrm{sgn}\,}}(V){{\,\textrm{sgn}\,}}(\beta )\right] ,\end{aligned}$$

(25b)

$$\begin{aligned}{} & {} \frac{d\sigma }{d\tau } = \frac{3}{n}\sigma ^{2}. \end{aligned}$$

(25c)

One can notice that the dynamical system (25) is symmetric under reflection against the $y=0$ plane, i.e. under the transformation $y\rightarrow -y$. This implies that it suffices to take into consideration only the $y>0$ region of the phase space, since the phase portrait in the $y<0$ region will just be a reflection of that against $y=0$. From the physical point of view, this means that the qualitative behaviour of the model is independent of the signature of the potential and we can confine our attention to ${{\,\textrm{sgn}\,}}(V)=1$. Also, the $y=0$ line acts as an invariant submanifold, meaning that dynamics taking place in the positive branch of the potential can never cross into the negative branch of the potential, and vice versa. Furthermore, for the physical viability of the model, one requires the condition $\Omega _{m}\ge 0$. The phase space of the system (25) is therefore constrained within the region given by

$$\begin{aligned} \lbrace (x,y,\sigma ) \in \mathbb {R}^3: y \ge 0, x^2-y^2-{{\,\textrm{sgn}\,}}(\beta )\ge 0\rbrace . \end{aligned}$$

(26)

Table 1 Existence and physical viability conditions for finite fixed points for a scalar field with kinetic term $F(X)=\beta X$ and potential $V(\phi )=V_0 \phi ^n$, calculated from the system (25)

Table 2 Stability condition of physically viable fixed points given in Table 1 along with their cosmological behavior. Here and throughout the paper, NH stands for nonhyperbolic

The system (25) contains two finite fixed points viz. $A_{1+}$ and $A_{1-}$ (see Table 1). Both points exist only in the case of the canonical scalar field, i.e., $\beta >0$. From the value of x-coordinate and the nature of scale factor a(t) (see Table 2), we see that while point $A_{1+}$ corresponds to a decelerated expansion of the universe, point $A_{1-}$ corresponds to a decelerated contraction of the universe. We note here that point $A_{1-}$ is non-hyperbolic with an empty unstable subspace near a point and therefore, we refer to the center manifold theory for further analysis. On performing the analysis, we obtained that the point $A_{1-}$ is a saddle point (see Appendix A). We also find that the trajectories flow from an unstable point $A_{1+}$ towards a saddle point $A_{1-}$. Therefore, from the finite analysis one can not extract any bouncing solution but it shows recollapsing solutions (see Fig. 1).

Fig. 1

Phase portrait of system (25) for $\omega _{m}=0$, $n=4$ and $\beta >0$

Fixed points at infinity

To get a global picture of the phase space, we introduce the following compact dynamical variables

$$\begin{aligned} \bar{x}=\frac{x}{\sqrt{1+x^{2}}}, \quad \bar{y}=\frac{y}{\sqrt{1+y^{2}}} , \quad \bar{\sigma }=\frac{\sigma }{\sqrt{1+\sigma ^{2}}}. \end{aligned}$$

(27)

The evolution Eq. (19) can be converted to the following system of equations

$$\begin{aligned} \frac{d\bar{x}}{d\tau }= & {} \frac{3}{2}(1-\bar{x}^{2})^{\frac{3}{2}}\left[ (\omega _{k}-\omega _{m}) \left( \frac{\bar{x}^{2}}{1-\bar{x}^{2}}-{{\,\textrm{sgn}\,}}(\beta )\right) \right. \nonumber \\{} & {} \left. +\left( 1+\omega _{m}-\frac{\bar{\sigma }\bar{x}{{\,\textrm{sgn}\,}}(\beta )}{\sqrt{1-\bar{\sigma }^2}\sqrt{1-\bar{x}^{2}}}\right) \frac{\bar{y}^{2}{{\,\textrm{sgn}\,}}(V)}{1-\bar{y}^{2}}\right] , \nonumber \\ \end{aligned}$$

(28a)

$$\begin{aligned} \frac{d\bar{y}}{d\tau }= & {} \frac{3}{2}\bar{y}(1-\bar{y}^{2})\left[ -\frac{\bar{\sigma }}{\sqrt{1-\bar{\sigma }^2}} +\frac{(\omega _{k}+1)\bar{x}}{\sqrt{1-\bar{x}^{2}}}\right. \nonumber \\{} & {} \left. -\frac{\bar{\sigma } \bar{y}^{2}{{\,\textrm{sgn}\,}}(V)}{\sqrt{1-\bar{\sigma }^2}(1-\bar{y}^{2})}{{\,\textrm{sgn}\,}}(\beta )\right] , \end{aligned}$$

(28b)

$$\begin{aligned} \frac{d\bar{\sigma }}{d\tau }= & {} \frac{3}{n} \bar{\sigma }^2\sqrt{1-\bar{\sigma }^2}. \end{aligned}$$

(28c)

Note that the variables $\bar{x}, \bar{y}, \bar{\sigma }$ are bounded between $-1$ and 1. The dynamical system in Eq. (28) is however not regular at the boundaries of the compact phase space $\bar{x}^2=1,\,\bar{y}^2=1$ and $\bar{\sigma }^2=1$, it has a pole of order $\frac{1}{2}$ at $\bar{x}^2=1$ and $\bar{\sigma }^2=1$ and a pole of order 1 at $\bar{y}^2=1$. This can be regularized following a prescription from Ref. [23]. The idea is to redefine the phase space time variable as

$$\begin{aligned} d\tau \rightarrow d\bar{\tau } = \frac{d\tau }{(1-\bar{x}^{2})^{\frac{1}{2}}(1-\bar{y}^{2})(1-\bar{\sigma }^2)^\frac{1}{2}}. \end{aligned}$$

(29)

With respect to this redefined time variable, the dynamical system (28) can be rewritten as

$$\begin{aligned} \frac{d\bar{x}}{d\bar{\tau }}= & {} \frac{3}{2}(\omega _{k}-\omega _{m})\left( \bar{x}^{2}-(1-\bar{x}^{2}){{\,\textrm{sgn}\,}}(\beta )\right) \nonumber \\{} & {} \times (1-\bar{x}^{2})(1-\bar{y}^{2}) \sqrt{1-\bar{\sigma }^2} \nonumber \\{} & {} + \frac{3}{2}\left( (1+\omega _{m})(1-\bar{x}^{2}) \sqrt{1-\bar{\sigma }^2}\right. \nonumber \\{} & {} \left. -\bar{\sigma }\bar{x}{{\,\textrm{sgn}\,}}(\beta )\sqrt{1-\bar{x}^{2}}\right) (1-\bar{x}^{2})\bar{y}^{2}, \end{aligned}$$

(30a)

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= & {} \frac{3}{2}\bar{y}(1-\bar{y}^{2})\left[ (-\bar{\sigma }\sqrt{1-\bar{x}^{2}}\right. \nonumber \\{} & {} \left. +(\omega _{k}+1)\bar{x}\sqrt{1-\bar{\sigma }^2})(1-\bar{y}^{2})\right. \nonumber \\{} & {} \left. -\bar{\sigma }\bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )\right] , \end{aligned}$$

(30b)

$$\begin{aligned} \frac{d\bar{\sigma }}{d\bar{\tau }}= & {} \frac{3}{n}\bar{\sigma }^2(1-\bar{\sigma }^2) \sqrt{1-\bar{x}^2}(1-\bar{y}^2). \end{aligned}$$

(30c)

Table 3 Existence and physical viability condition for fixed points at infinity for a scalar field with kinetic term $F(X)=\beta X$ and potential $V(\phi )=V_0 \phi ^n$ calculated from the system (30)

Table 4 Stability condition of physically viable fixed points given in Table 3 along with their cosmological behavior. Stability of the lines of fixed points $B_{1\pm },\,B_{2\pm }$ and $B_{3\pm }$ can be determined by investigating the stability of the invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=0,1$ and $\bar{\sigma }=\pm 1$ (see Appendix B)

In terms of the compact variables, the constraint equation (19a) can be rewritten for this model as

$$\begin{aligned} \frac{\bar{x}^{2}}{1-\bar{x}^{2}} - \frac{\bar{y}^{2}}{1-\bar{y}^{2}} - {{\,\textrm{sgn}\,}}(\beta ) = \Omega _{m} . \end{aligned}$$

(31)

Therefore, for a canonical scalar field where ${{\,\textrm{sgn}\,}}(\beta )=1$, the physical viability condition requires

$$\begin{aligned} \frac{\bar{x}^{2}}{1-\bar{x}^{2}}-\frac{\bar{y}^{2}}{1-\bar{y}^{2}}-1=\Omega _{m}\ge 0 . \end{aligned}$$

(32)

It can be checked that the necessary and sufficient conditions for physical viability of canonical scalar field are given respectively as follows

$$\begin{aligned} \bar{y}^{2}\le 2-\frac{1}{\bar{x}^{2}}, \quad \bar{x}^{2}\ge \frac{1}{2}. \end{aligned}$$

(33)

For a phantom scalar field where ${{\,\textrm{sgn}\,}}(\beta )=-1$, the physical viability condition requires

$$\begin{aligned} \frac{\bar{x}^{2}}{1-\bar{x}^{2}}-\frac{\bar{y}^{2}}{1-\bar{y}^{2}}+1=\Omega _{m}\ge 0. \end{aligned}$$

(34)

The necessary and sufficient condition for the above to be satisfied is

$$\begin{aligned} \bar{y}^{2}\le \frac{1}{2-\bar{x}^{2}}. \end{aligned}$$

(35)

The constraints in Eqs. (33) and (35) specify the physically viable region of the entire 3-dimensional compact phase space for canonical and phantom scalar fields respectively. One can notice that the line $\bar{x}=0$ is physical only for the case of a phantom scalar field ($\beta <0$), which is consistent with the fact that a phantom scalar field is necessarily required to achieve a nonsingular bounce.

The system (30) presents a total of six invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=0,1$ and $\bar{\sigma }=\pm 1$. Their stability is calculated in Appendix B. The fixed points for system (30) are presented in Table 3, along with their stability conditions in Table 4. Although from Table 3 it might appear that there are five different pairs of fixed points at the infinity of the phase space, the pair $B_{4\pm }$ is physically viable only for $\bar{\sigma }=0$, in which case they already lie on the line of fixed points $B_{1\pm }$. Therefore there are only four different pairs of fixed points at infinity. $B_{1\pm }$ are De-Sitter solutions, whereas $B_{2\pm }$ correspond to cosmological phases dominated by the hydrodynamic matter component. For $\omega _{m}=0$, the fixed point $B_{2+}$ can be interpreted as the matter-dominated epoch since this is a saddle and therefore, always represents an intermediate phase of evolution. Fixed points on the line $B_{3\pm }$ can exhibit various cosmological scenarios. For instance, they can describe a matter-dominated universe for $\bar{x}=1,\,\beta >0$, static universe for $\bar{x}=0,\,\beta <0$. Finally, the fixed points $B_{5\pm }$ correspond to a static universe.

In the left panel of Fig. 2, we present the phase portrait of the system (30) in the compact phase space that corresponds to physically viable nonsingular bouncing solutions. We also observe that there is a violation of the null energy condition (NEC) during bounce (see the right panel of Fig. 2), as expected. The plots also clearly show the matter-dominated fixed points $B_{2\pm }$ are saddles, i.e. intermediate epochs of evolution, as expected.

The phase trajectories in the upper left panel, which corresponds to $F(X)=\beta X\,(\beta <0),\,V(\phi )=V_{0}\phi ^4$, show nonsingular bouncing solutions connecting the contracting and expanding De-Sitter phases $B_{1-}$ and $B_{1+}$ respectively. However, one cannot say this is the generic behaviour of phase trajectories as $B_{1-}$ and $B_{1+}$ are not global repellers and attractors. Apart from $B_{1-}$ and $B_{1+}$ there is also another repeller $B_{3-}$ and another attractor $B_{3+}$. Therefore, there can exist heteroclinic trajectories connecting $B_{3-}$ to $B_{1+}$, $B_{1-}$ to $B_{3+}$ and $B_{3-}$ to $B_{3+}$, which may or may not correspond to a nonsingular bouncing cosmology, depending on how the trajectories evolve.

The phase trajectories in the lower left panel, which corresponds to $F(X)=\beta X\,(\beta <0),\,V(\phi )=V_{0}\phi ^{-4}$, show nonsingular bouncing solutions connecting two De-Sitter phases ($B_{1-} \rightarrow B_{1+}$) or a De-Sitter phase with a static universe ($B_{1-} \rightarrow B_{5-}$ or $B_{5+} \rightarrow B_{1+}$). However, again, one cannot say that such behaviour is generic. For $n=-4<0$, the fixed points $B_{3\pm }$ are saddles, but $B_{5\pm }$ becomes an attractor/repeller pair. Thus, there can exist heteroclinic trajectories connecting $B_{5-}$ to $B_{5+}$, which, again, may or may not correspond to a nonsingular bouncing cosmology, depending on how the trajectories evolve.

Fig. 2

Phase portrait of the system (30) with the shaded region representing the non-physical region of the phase space (upper and lower left panels). Here $\omega _{m}=0$, $\beta <0$ with $n=4$ in the upper left panel and $n=-4$ in the lower panel. A plot of variables $\bar{x}$ and effective equation of state $\omega _\textrm{eff}$ is shown in the upper and lower right panel for the blue trajectories of the upper and lower left panel, which represents two characteristic types of nonsingular bouncing solutions, namely, asymptotically De-Sitter and asymptotically static in future

3.1.2 Exponential potential $V(\phi )=V_{0}e^{-\lambda \phi /M_{Pl}}$:

The second example we consider is the specific case given by $F(X)=\beta X,\,V(\phi )=V_{0}e^{-\lambda \phi /M_{Pl}}$ where $\beta , V_0$ are constants with suitable dimension, $\lambda$ is a dimensionless constant. For this choice, we have

$$\begin{aligned} \Xi =0, \quad \Gamma =1, \quad \omega _{k}=1, \quad \rho _{k}=\beta X, \quad \sigma = \frac{\sqrt{{2}/{3}}\lambda }{\sqrt{|\beta |}}. \end{aligned}$$

(36)

Finite Fixed Point Analysis

In this case, since $\sigma$ is a constant, the system (19) reduces to a 2-dimensional dynamical system

$$\begin{aligned} \frac{dx}{d\tau }= & {} \frac{3}{2}x\left[ 2x-\sigma y^{2}{{\,\textrm{sgn}\,}}(V_{0}){{\,\textrm{sgn}\,}}(\beta )\right] \nonumber \\{} & {} - \frac{3}{2}\left[ (1-\omega _{m}){{\,\textrm{sgn}\,}}(\beta )+(1+\omega _{m})(x^{2}-y^{2}{{\,\textrm{sgn}\,}}(V_{0}))\right] ,\nonumber \\ \end{aligned}$$

(37a)

$$\begin{aligned} \frac{dy}{d\tau }= & {} \frac{3}{2}y\left[ 2x-\sigma \left( 1+ y^{2}{{\,\textrm{sgn}\,}}(V_{0}){{\,\textrm{sgn}\,}}(\beta )\right) \right] . \end{aligned}$$

(37b)

As in the power law case, one can notice that the dynamical system (37) is symmetric under reflection against the $y=0$ line, i.e. under the transformation $y\rightarrow -y$. This implies that it suffices to take into consideration only the $y>0$ region of the phase space, since the phase portrait in the $y<0$ region will just be a reflection of that against $y=0$. From the physical point of view, this means that the qualitative behavior of the model is independent of the signature of the potential and we can confine our attention to ${{\,\textrm{sgn}\,}}(V_0)=1$. Therefore, the interpretation of the $y=0$ invariant submanifold is very clear in this case. As before, for the physical viability of the model, one requires the condition $\Omega _{m}\ge 0$. The phase space of the system (25) is therefore constrained within the region given by

$$\begin{aligned}{} & {} \lbrace (x,y) \in \mathbb {R}^2: y \ge 0, x^2-y^2-{{\,\textrm{sgn}\,}}(\beta )\ge 0\rbrace . \end{aligned}$$

(38)

Table 5 Existence and physical viability conditions for finite fixed points for a scalar field with kinetic term $F(X)=\beta X$ and potential $V(\phi )=V_0 e^{-\lambda \phi /M_{Pl}}$, calculated from the system (37)

Table 6 Stability condition of physically viable fixed points given in Table 5 along with their cosmological behavior. Stability of $A_{2\pm },\,A_{3\pm }$ and $A_4$ are determined from the Jacobian eigenvalues whereas the stability of $A_{1\pm }$ are determined by examining the stability of the invariant submanifolds $y=0$ and $\Omega _m=0$ (see Appendix B)

The system (37) contains seven finite fixed points viz. $A_{1\pm },\,A_{2\pm },\,A_{3\pm }$ and $A_{4}$ (see Table 5). The stability and cosmological evolution corresponding to each fixed point are listed in Table 6. The fixed points $A_{1\pm }$ that exist only for the case of a canonical scalar field, correspond to decelerated expanding and contracting phases respectively, which we also obtained for the power law potential. The rest of the finite fixed points obtained for the exponential potential have no counterpart for the power law potential. For the case of a canonical scalar field, we get another pair of expanding and contracting solutions, $A_{3\pm }$, which can be accelerated or decelerated according to $\sigma ^2>\frac{4}{3}$ or $\sigma ^2<\frac{4}{3}$ respectively. $A_{3\pm }$ coincide with $A_{1\pm }$ in the limit $|\sigma |\rightarrow 2$. Apart from that, for the canonical case, we also have an expanding power law solution $A_4$, which is accelerated or decelerated according to $1+3\omega _{m}>0$ or $1+3\omega _{m}<0$, i.e. whether the matter component satisfies the SEC or not. For the case of a phantom scalar field, we get two solutions $A_{2\pm }$, which represent finite time singularities in the past and the future respectively. In fact, $A_{2+}$ is a phantom-dominated phase that ends, as expected, in a big-rip singularity.

Fixed points at infinity

To get a global picture of the 2-dimensional phase space we employ the compact dynamical variables $\bar{x},\,\bar{y}$ as in Eq. (27). For a constant $\sigma$ Eq. (28) reduces to

$$\begin{aligned} \frac{d\bar{x}}{d\tau }= & {} \frac{3}{2}(1-\bar{x}^{2})^{\frac{3}{2}}\left[ (\omega _{k}-\omega _{m}) \left( \frac{\bar{x}^{2}}{1-\bar{x}^{2}}-{{\,\textrm{sgn}\,}}(\beta )\right) \right. \nonumber \\{} & {} \left. +\left( 1+\omega _{m}-\frac{\sigma \bar{x}{{\,\textrm{sgn}\,}}(\beta )}{\sqrt{1-\bar{x}^{2}}}\right) \frac{\bar{y}^{2}}{1-\bar{y}^{2}}\right] ,\nonumber \\ \end{aligned}$$

(39a)

$$\begin{aligned} \frac{d\bar{y}}{d\tau }= & {} \frac{3}{2}\bar{y}(1-\bar{y}^{2})\left[ -\sigma +\frac{(\omega _{k}+1)\bar{x}}{\sqrt{1-\bar{x}^{2}}}-\frac{\sigma \bar{y}^{2}}{1-\bar{y}^{2}}{{\,\textrm{sgn}\,}}(\beta )\right] .\nonumber \\ \end{aligned}$$

(39b)

As in the previous case, to regularize the dynamical system Eq. (39) we redefine the phase space time variable as

$$\begin{aligned} d\tau \rightarrow d\bar{\tau } = \frac{d\tau }{(1-\bar{x}^{2})^{\frac{1}{2}}(1-\bar{y}^{2})}. \end{aligned}$$

(40)

With respect to this redefined time variable, the dynamical system (39) can be rewritten as

$$\begin{aligned}{} & {} \frac{d\bar{x}}{d\bar{\tau }} = \frac{3}{2}(\omega _{k}-\omega _{m}) \left( \bar{x}^{2}-(1-\bar{x}^{2}){{\,\textrm{sgn}\,}}(\beta )\right) \nonumber \\{} & {} \quad \times (1-\bar{x}^{2})(1-\bar{y}^{2})+ \frac{3}{2}\Big ((1+\omega _{m})(1-\bar{x}^{2})-\sigma \bar{x}\nonumber \\{} & {} \qquad \times {{\,\textrm{sgn}\,}}(\beta )\sqrt{1-\bar{x}^{2}}\Big )(1-\bar{x}^{2})\bar{y}^{2}, \end{aligned}$$

(41a)

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}{} & {} = \frac{3}{2}\bar{y}(1-\bar{y}^{2})\left[ (-\sigma \sqrt{1-\bar{x}^{2}} +(\omega _{k}+1)\bar{x})(1-\bar{y}^{2})\right. \nonumber \\{} & {} \quad \left. -\sigma \bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )\right] . \end{aligned}$$

(41b)

The constraints in Eqs. (33) and (35), which specified the physically viable region of the phase space for power law potential, remain valid also for exponential potential. This implies we can only focus on phantom scalar field ($\beta <0$) in whatever follows if we want to achieve nonsingular bounces.

Table 7 Existence and physical viability condition for fixed points at infinity for a scalar field with kinetic term $F(X)=\beta X$ and potential $V(\phi )=V_0 e^{-\lambda \phi /M_{Pl}}$, calculated from the system (41)

Table 8 Stability condition of physically viable fixed points given in Table 7 along with their cosmological behavior. The stability of these fixed points is determined by investigating the stability of the invariant submanifolds $\bar{x}=\pm 1$ and $\bar{y}=0,1$ (see Appendix B)

Fixed points for the system (41) are listed in Table 7, along with their stability conditions in Table 8. Although from Table 7 it might appear that there are three different pairs of fixed points at the infinity of the phase space, the pair $B_{3\pm }$ is physically viable only for $\sigma =0$, in which case it coincides with $B_{1\pm }$. Therefore there are only two different pairs of fixed points at infinity. $B_{1\pm }$ are De-Sitter solutions whereas $B_{2\pm }$ corresponds to cosmological phases dominated by the hydrodynamic matter component. For $\omega _{m}=0$ the fixed point $B_{2+}$ can be interpreted as the matter-dominated epoch since this is a saddle and therefore always represents an intermediate phase of evolution. It might also be noted that $A_{2\pm }$ and $A_{3\pm }$ merges with $B_{1\pm }$ at the limit $|\sigma |\rightarrow \infty$.

Below we provide compact 2-dimensional phase portraits for various cases. In all the plots we also specify the physically viable region as calculated by the constraints (33) or (35) and the regions where the NEC and SEC are satisfied. The NEC and SEC are determined in terms of the “effective” equation of state parameter defined in Eq. (22).

Fig. 3

The phase space portrait for $\omega _{m}=0$ with $\sigma =0$, $\beta <0$ and $V_{0}>0$. Satisfaction of the NEC is shown in light green, the satisfaction of the SEC in dark green and the red area denotes a non-physical part of the phase space

Fig. 4

The phase space portrait for $\sigma =0.8$ (upper left panel) $\sigma =-0.8$ (upper right panel) $\sigma =3$ (lower left panel) $\sigma =-3$ (lower right panel) with $\omega _{m}=0$, $\beta <0$ and $V_{0}>0$. Satisfaction of the NEC is shown in light green, the satisfaction of the SEC in dark green and the red area denotes a non-physical part of the phase space

Fig. 5

The phase space portrait for $\omega _{m}=0$ with $\sigma =0$, $\beta >0$ and $V_{0}>0$. The shaded region represents the non-physical viable region

From Figs. 3 and 4, it can be seen that for a phantom scalar field ($\beta <0$), all the phase trajectories within the physically viable region show a nonsingular bounce. Nonsingular bounce is a generic feature in this case. the figures also confirm the violation of both NEC and SEC to achieve a bounce. For the sake of completeness, we also show the phase portrait for a particular case for $\beta >0$ in Fig. 5, which shows that all the trajectories representing a nonsingular bounce fall within the physically non-viable region, as expected.

3.2 Specific case II: noncanonical scalar field $F(X)=\beta X^{m}$ ($m\ne 1$)

In this section, we generalize our scalar field Lagrangian to the noncanonical form

$$\begin{aligned} \mathcal {L}(\phi ,X) = \beta X^m - V(\phi ), \end{aligned}$$

(42)

Since $\rho _{\phi }+P_{\phi }=2XF_X=2m\beta X^m$ and $X\ge 0$, it corresponds to a non-phantom scalar field when $m\beta >0$ and to a phantom scalar field when $m\beta <0$. Again, one can notice from Eq. (8b) that for $\dot{H}>0$ near the bounce one must necessarily need $m\beta <0$, i.e. a phantom scalar field. We particularly concentrate on the case $m\ne 1$, as the case $m=1$ has been considered previously.

The analysis in this paper relies on the scalar field being a dynamical degree of freedom. The case $m=1/2$, i.e. $F(X)\propto \pm \sqrt{X}$ is a very special case for which the scalar field becomes non-dynamical in the homogeneous limit [24, 25]. In this case, the field is called Cuscuton. Since the field is nondynamical, violation of NEC during a bounce does not lead to a pathology [26]. This case is also left out of consideration in this paper. For more about Cuscustun bounce cosmology (See Refs. [27,28,29]).

One can calculate that, in this case,

$$\begin{aligned} \rho _{k}=(2m-1)\beta X^{m}. \end{aligned}$$

(43)

We keep on using the same compact and non-compact dynamical variables as introduced earlier. The physical viability of the model requires the condition $\Omega _{m}\ge 0$. In terms of the non-compact dynamical variables (x, y), it can be written as

$$\begin{aligned} \lbrace (x,y,\sigma ) \in \mathbb {R}^3: x^2-y^2-{{\,\textrm{sgn}\,}}((2m-1)\beta )\ge 0\rbrace . \end{aligned}$$

(44)

In terms of the compact dynamical variables $(\bar{x},\bar{y})$ defined in (27), one can write

$$\begin{aligned} \frac{\bar{x}^{2}}{1-\bar{x}^{2}}-\frac{\bar{y}^{2}}{1-\bar{y}^{2}}-{{\,\textrm{sgn}\,}}{((2m-1)\beta )}=\Omega _{m}. \end{aligned}$$

(45)

Note that the physical viability conditions for $m\ne 1$ differ from that for $m=1$ only in the fact that $sgn(\beta )$ is replaced by $sgn((2m-1)\beta )$. When ${{\,\textrm{sgn}\,}}((2m-1)\beta )=1$ one needs

$$\begin{aligned} \frac{\bar{x}^{2}}{1-\bar{x}^{2}}-\frac{\bar{y}^{2}}{1-\bar{y}^{2}}-1=\Omega _{m}\ge 0. \end{aligned}$$

(46)

This leads to the necessary and sufficient conditions for the physical viability respectively as

$$\begin{aligned} \bar{y}^{2}\le 2-\frac{1}{\bar{x}^{2}}, \qquad \bar{x}^{2}\ge \frac{1}{2}, \end{aligned}$$

(47)

which are the same conditions obtained in Eq. (33). When ${{\,\textrm{sgn}\,}}((2m-1)\beta )=-1$ one needs

$$\begin{aligned} \frac{\bar{x}^{2}}{1-\bar{x}^{2}}-\frac{\bar{y}^{2}}{1-\bar{y}^{2}}+1=\Omega _{m}\ge 0. \end{aligned}$$

(48)

This leads to the necessary and sufficient conditions as

$$\begin{aligned} \bar{y}^{2}\le \frac{1}{2-\bar{x}^{2}}, \end{aligned}$$

(49)

which is the same condition as obtained in Eq. (35). Just like one could conclude for the case of $F(X)=\beta X$ that a physically viable nonsingular bounce requires $\beta <0$, one can conclude here that for the generic case $F(X)=\beta X^m$, a physically viable nonsingular bounce requires $(2m-1)\beta <0$. This implies the following parameter range

$$\begin{aligned} \{(m>1/2)\wedge (\beta<0)\} \vee \{(m<1/2)\wedge (\beta >0)\}. \end{aligned}$$

(50)

As we have already seen from (8b) that achieving a nonsingular bounce necessarily requires a phantom scalar field, i.e. $m\beta <0$, together with the condition $(2m-1)\beta <0$ this slightly constrain the parameter range

$$\begin{aligned} \{(m>1/2)\wedge (\beta<0)\} \vee \{(m<0)\wedge (\beta >0)\}. \end{aligned}$$

(51)

In summary, the parameter range $0\le m< 1/2$ is not allowed if we want to achieve a nonsingular bounce and the Cuscuton case $m=1/2$ is not taken into consideration.

3.2.1 Power law potential $V(\phi )=V_{0}\phi ^{n}$

For the kinetic term $F(X)=\beta X^{m}$ and the potential $V(\phi )=V_{0}\phi ^{n}$, we have

$$\begin{aligned}&\Xi =m-1, \quad \Gamma =1-\frac{1}{n},\nonumber \\ \omega _{k}&=\frac{1}{2m-1}, \quad \rho _{k}=(2m-1)\beta X^{m}. \end{aligned}$$

(52)

Finite fixed point analysis

In this case, the system (19) reduces to

$$\begin{aligned} \frac{dx}{d\tau }= & {} \frac{3}{2}x\left[ \left( \frac{2m}{2m-1}\right) x - \sigma y^2 {{\,\textrm{sgn}\,}}(V){{\,\textrm{sgn}\,}}[(2m-1)\beta ]\right] \nonumber \\{} & {} - \frac{3}{2}\left[ \left( \frac{1}{2m-1}-\omega _{m}\right) {{\,\textrm{sgn}\,}}[(2m-1)\beta ] \right. \nonumber \\{} & {} \left. + (1+\omega _{m})(x^{2}-y^2 {{\,\textrm{sgn}\,}}(V)) \right] , \end{aligned}$$

(53a)

$$\begin{aligned} \frac{dy}{d\tau }= & {} \frac{3}{2}y\left[ -\sigma + \left( \frac{2m}{2m-1}\right) x - \sigma y^2\right. \nonumber \\{} & {} \quad \left. {{\,\textrm{sgn}\,}}(V){{\,\textrm{sgn}\,}}[(2m-1)\beta ] \right] , \end{aligned}$$

(53b)

$$\begin{aligned} \frac{d\sigma }{d\tau }= & {} \frac{3}{n}\sigma ^{2}+3\sigma \frac{(2m-3)m+1}{(4m-2)m}\nonumber \\{} & {} \quad \times \left( \left( \frac{2m}{2m-1}\right) x-\sigma y^{2}{{\,\textrm{sgn}\,}}(V)\right) . \end{aligned}$$

(53c)

The above system reduces to the system (25) for $m=1$. As in the case of $m=1$, the system is symmetric under reflection around $y=0$, which happens to be an invariant submanifold. Therefore it suffices to consider only the part of the phase space given by

$$\begin{aligned} \lbrace (x,y,\sigma ) \in \mathbb {R}^3: y \ge 0, x^2-y^2-{{\,\textrm{sgn}\,}}((2m-1)\beta )\ge 0\rbrace . \end{aligned}$$

(54)

Table 9 Existence and physical viability conditions for finite fixed points for a non-canonical scalar field with kinetic term $F(X)=\beta X^m$ ($m\ne 1$) and potential $V(\phi )=V_0 \phi ^n$, calculated from the system (53)

Table 10 Stability condition of physically viable critical points given in Table 9 along with their cosmological behaviour

The system (53) contains two finite fixed points $A_{1\pm }$, which exist only for non-phantom fields (see Table 9). The stabilities and corresponding cosmologies are given in Table 10. Since these critical points exist only when $(2m-1)\beta >0$. When $(2m-1)\beta >0$, they cannot give rise to physically viable nonsingular bouncing trajectories.

Fixed points at infinity

In terms of the compact dynamical variables $(\bar{x},\bar{y},\bar{\sigma })$ defined in Eq. (27) and the redefined time variable defined in Eq. (29), the dynamical system for $F(X)=\beta X^{m},\,V=V_0\phi ^n$ can be rewritten in the following form

$$\begin{aligned} \frac{d\bar{x}}{d\bar{\tau }}= & {} \frac{3}{2}\left( \frac{1}{2m-1}-\omega _{m}\right) \left( \bar{x}^{2}-(1-\bar{x}^{2}){{\,\textrm{sgn}\,}}[(2m-1)\beta ]\right) \nonumber \\{} & {} \times (1-\bar{x}^{2})(1-\bar{y}^{2})\sqrt{1-\bar{\sigma }^2}\nonumber \\{} & {} +\frac{3}{2}\big ((1+\omega _{m}) (1-\bar{x}^{2}) \sqrt{1-\bar{\sigma }^2}\nonumber \\{} & {} -\bar{\sigma }\bar{x}{{\,\textrm{sgn}\,}}[(2m-1)\beta ]\sqrt{1-\bar{x}^{2}}\big )(1-\bar{x}^{2})\bar{y}^{2}, \end{aligned}$$

(55a)

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= & {} \frac{3}{2}\bar{y}(1-\bar{y}^{2})\left[ \left( -\bar{\sigma }\sqrt{1-\bar{x}^{2}} +\left( \frac{2m}{2m-1}\right) \bar{x}\sqrt{1-\bar{\sigma }^2}\right) \right. \nonumber \\{} & {} \times \left. (1-\bar{y}^{2})-\bar{\sigma }\bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}[(2m-1)\beta ] \right] , \end{aligned}$$

(55b)

$$\begin{aligned} \frac{d\bar{\sigma }}{d\bar{\tau }}= & {} (1-\bar{\sigma }^{2})\left[ \frac{3}{n}\bar{\sigma }^{2} \sqrt{1-\bar{x}^{2}}(1-\bar{y}^{2})\right. \nonumber \\{} & {} \left. +3\bar{\sigma }\frac{(2m-3)m+1}{(4m-2)m}\left( \left( \frac{2m}{2m-1}\right) \bar{x}\sqrt{1-\bar{\sigma }^{2}}(1-\bar{y}^{2})\right. \right. \nonumber \\{} & {} \left. \left. -\bar{\sigma }\bar{y}^{2}\sqrt{1-\bar{x}^{2}}\right) \right] . \end{aligned}$$

(55c)

Table 11 Existence and physical viability condition for fixed points at infinity for a noncanonical scalar field with kinetic term $F(X)=\beta X^m$ ($m\ne 1$) and potential $V(\phi )=V_0 \phi ^n$, calculated from the system (55)

The system (55) presents six invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=0,1$ and $\bar{\sigma }=\pm 1$. Their stability is calculated in Appendix C). The fixed points for the system (55) are presented in Table 11. In the presence of pressureless dust, their stability conditions and corresponding cosmologies are presented in Table 12. The lines of fixed points $B_{1\pm },\,B_{3\pm }$ and the isolated fixed points $B_{5\pm }$ are the same ones that we obtained earlier for the case $m=1$ (see Table 3). However, unlike in the case of $m=1$, the entire lines $B_{2\pm }$ and $B_{4\pm }$ are not lines of fixed points when $m\ne 1$. Instead, in this case, we only get two isolated fixed points $B^a_{2\pm }$ that lie on the line $B_{2\pm }$ respectively, and six isolated fixed points $B^{a,b,c}_{4\pm }$ that lie on the line $B_{4\pm }$ respectively. The points are listed in Table 11. Moreover, $B_{4\pm }^{a,b,c}$ are physically viable only for $\bar{\sigma }=0$, in which case they fall back into the lines of fixed points $B_{1\pm }$ respectively. Therefore at the infinity of the phase space there are only two lines of fixed points $B_{1\pm },\,B_{3\pm }$ and two pairs of isolated fixed points $B_{2\pm }^a$ and $B_{5\pm }$, whose nature of stability and the corresponding cosmology is listed in Table 12. The cosmologies are the same as we have earlier obtained for $m=1$, as expected. Figure 6 presents the 3D phase portrait of the system (55) in the compact space for two cases, showing different types of physically viable nonsingular bouncing trajectories.

Table 12 Stability condition and the cosmology of the fixed points given in Table 11 in presence of pressureless dust $(\omega _m=0)$. Stability of the lines of fixed points $B_{1\pm }$ and $B_{3\pm }$ can be determined by investigating the stability of the invariant submanifolds $\bar{x}=\pm 1$, $\bar{y}=0,1$ and $\bar{\sigma }\pm 1$ (see Appendix C)

The phase trajectories in the right panel of Fig. 6 shows nonsingular bouncing trajectories for the case $F(X)=\beta X^3\,(\beta <0),\,V(\phi )=V_{0}\phi ^4$. Several different types of bouncing trajectories are possible, e.g. trajectories connecting a contracting De-Sitter phase to an expanding De-Sitter phase ($B_{1-}\rightarrow B_{1+}$), trajectories connecting a static and a De-Sitter phase ($B_{1-}\rightarrow B_{5-}$ and $B_{1+}\rightarrow B_{5+}$) and trajectories connecting $B_{3\pm }$ with $B_{1\pm }$ or $B_{5\pm }$. However, one cannot say this is the generic behaviour of phase trajectories for $F(X)=\beta X^3\,(\beta <0),\,V(\phi )=V_{0}\phi ^4$, as $B_{1-}$ and $B_{1+}$ are not global repellers and attractors. For $m=3,\,n=4$, there are two other attractor/repeller pairs $B_{3\pm }$ and $B_{5\pm }$. Therefore, there can exist heteroclinic trajectories connecting them, which may or may not correspond to a nonsingular bouncing cosmology, depending on how the trajectories evolve.

The phase trajectories in the left panel of Fig. 6, which corresponds to $F(X)=\beta X^{2/3}\,(\beta <0),\,V(\phi )=V_{0}\phi ^{-5}$, show nonsingular bouncing trajectories connecting two De-Sitter phases ($B_{1-} \rightarrow B_{1+}$). For $m=\frac{2}{3},\,n=-5$, $B_{1+}$ and $B_{1-}$ are the only attractor/repeller pair possible, i.e., they are global attractors and repellers. The trajectories connecting them must necessarily undergo a nonsingular bounce and we have shown explicitly in the figure four such characteristic trajectories. In this case, one can confidently say that a nonsingular bounce is a generic behaviour of the phase trajectories.

Fig. 6

Phase portrait of the system (55) with shaded region represents a non-physical region of the phase space. The left plot is for $n=-5, m=\frac{2}{3}$ and the right plot is for $n=4, m=3$ moreover ${{\,\textrm{sgn}\,}}((2m-1)\beta )=-1$

3.2.2 Exponential potential $V(\phi )=V_0 e^{-\lambda \phi /M_{Pl}}$

For the kinetic term $F(X)=\beta X^m$ and the potential $V(\phi )=V_0 e^{-\lambda \phi /M_{Pl}}$, we have

$$\begin{aligned} \Xi =m-1, \quad \Gamma =1, \quad \omega _{k}=\frac{1}{2m-1}, \quad \rho _{k}=(2m-1)\beta X^{m}\nonumber \\ \end{aligned}$$

(56)

Finite fixed point analysis

In this case, the system (19) reduces to

$$\begin{aligned} \frac{dx}{d\tau }= & {} \frac{3x}{2}\left[ \left( \frac{2m}{2m-1}\right) x - \sigma y^2{{\,\textrm{sgn}\,}}[(2m-1)\beta ]\right] \nonumber \\{} & {} - \frac{3}{2}\left[ \left( \frac{1}{2m-1}-\omega _{m}\right) {{\,\textrm{sgn}\,}}[(2m-1)\beta ]\right. \nonumber \\{} & {} \left. + (1+\omega _{m})(x^{2}-y^2) \right] , \end{aligned}$$

(57a)

$$\begin{aligned} \frac{dy}{d\tau }= & {} \frac{3}{2}y\left[ -\sigma + \left( \frac{2m}{2m-1}\right) x - \sigma y^2{{\,\textrm{sgn}\,}}[(2m-1)\beta ]\right] , \end{aligned}$$

(57b)

$$\begin{aligned} \frac{d\sigma }{d\tau }= & {} 3\sigma \frac{(2m-3)m+1}{(4m-2)m}\left[ \left( \frac{2m}{2m-1}\right) x-\sigma y^{2}\right] . \end{aligned}$$

(57c)

For $m=1$ the dynamical system for the exponential potential was 2D ((37)). For $m\ne 1$, the dynamical system for the exponential potential becomes 3D. The dynamical system is independent of the parameter $\lambda$. The above system reduces to the system (37) for $m=1$. The system is symmetric under reflection around $y=0$, which happens to be an invariant submanifold. Therefore it suffices to consider only the part of the phase space given by

$$\begin{aligned} \lbrace (x,y,\sigma ) \in \mathbb {R}^3: y \ge 0, x^2-y^2-{{\,\textrm{sgn}\,}}((2m-1)\beta )\ge 0\rbrace , \end{aligned}$$

(58)

which corresponds to taking $V_0>0$.

Table 13 Existence and the physical viability condition for finite fixed points for a noncanonical scalar field with kinetic term $F(X)=\beta X^m$ ($m\ne 1$) and potential $V(\phi )=V_{0}e^{-\lambda \phi /M_{Pl}}$, calculated from the system (57)

Table 14 Stability and cosmological behavior of the physically viable fixed points given in Table 13

The system (57) contains two finite fixed points $A_{1\pm }$, which exist only for non-phantom fields (see Table 13). The stabilities and corresponding cosmologies are given in Table 14. Since these critical points exist only when $(2m-1)\beta >0$. When $(2m-1)\beta >0$, they cannot give rise to physically viable nonsingular bouncing trajectories. The finite fixed points $A_{2\pm },\,A_{3\pm }$ and $A_4$ that appeared in Table 5, are specific to the case $m=1$ and does not arise for $m\ne 1$.

Fixed points at infinity

In terms of the compact dynamical variables $(\bar{x},\bar{y},\bar{\sigma })$ defined in Eq. (27) and the redefined time variable defined in Eq. (29), the dynamical system for $F(X)=\beta X^{m},\,V=V_0 e^{-\lambda \phi /M_{Pl}}$ can be rewritten in the following form

$$\begin{aligned} \frac{d\bar{x}}{d\bar{\tau }}= & {} \frac{3}{2}\left( \frac{1}{2m-1}-\omega _{m}\right) \left( \bar{x}^{2}-(1-\bar{x}^{2}){{\,\textrm{sgn}\,}}[(2m-1)\beta ]\right) \nonumber \\{} & {} \times (1-\bar{x}^{2})(1-\bar{y}^{2})\sqrt{1-\bar{\sigma }^2}\nonumber \\{} & {} +\frac{3}{2}\big ((1+\omega _{m}) (1-\bar{x}^{2}) \sqrt{1-\bar{\sigma }^2}\nonumber \\{} & {} -\bar{\sigma }\bar{x}{{\,\textrm{sgn}\,}}[(2m-1)\beta ]\sqrt{1-\bar{x}^{2}}\big )(1-\bar{x}^{2})\bar{y}^{2}, \end{aligned}$$

(59a)

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= & {} \frac{3}{2}\bar{y}(1-\bar{y}^{2}) \left[ \left( -\bar{\sigma }\sqrt{1-\bar{x}^{2}} +\left( \frac{2m}{2m-1}\right) \bar{x}\sqrt{1-\bar{\sigma }^2}\right) \right. \nonumber \\{} & {} \times \left. (1-\bar{y}^{2})-\bar{\sigma }\bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}[(2m-1)\beta ]\right] , \end{aligned}$$

(59b)

$$\begin{aligned} \frac{d\bar{\sigma }}{d\bar{\tau }}= & {} 3\bar{\sigma }(1-\bar{\sigma }^{2}) \left[ \frac{(2m-3)m+1}{(4m-2)m} \left( \left( \frac{2m}{2m-1}\right) \bar{x}\right. \right. \nonumber \\{} & {} \times \left. \left. \sqrt{1-\bar{\sigma }^{2}}(1-\bar{y}^{2})-\bar{\sigma }\bar{y}^{2}\sqrt{1-\bar{x}^{2}}\right) \right] . \end{aligned}$$

(59c)

Table 15 Existence and the physical viability condition for fixed points at infinity for a noncanonical scalar field with kinetic term $F(X)=\beta X^m$ ($m\ne 1$) and $V(\phi )=V_{0}e^{-\lambda \phi /M_{Pl}}$, calculated from the system (59)

The system (59) presents six invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=0,1$ and $\bar{\sigma }=\pm 1$, same as in the power law case. Their stability is calculated in Appendix C. The fixed points for the system (59) are presented in Table 15. In the presence of pressureless dust, their stability conditions and corresponding cosmologies are presented in Table 16. The six isolated fixed points $B_{3\pm }^{a,b,c}$ exist only for $\bar{\sigma }=0$, for which they fall back on the lines of fixed points $B_{1\pm }$. The $\{x,y\}$ coordinates of the fixed points $B_{1\pm }$ and $B_{2\pm }^c$ are the same as the fixed points $B_{1\pm }$ and $B_{2\pm }$ of the 2D phase space for the case $m=1$ (see Table 7). Since for the case $m\ne 1$, $\sigma$ is a dynamical variable, we get the fixed points $B_{2\pm }^{a,b}$ and $B_{5\pm }$ at the boundary $\sigma \rightarrow \pm \infty$. We did not get any corresponding point for the case $m=1$ because in that case, $\sigma$ was a parameter and we considered it to be finite. In Fig. 7 we present the 3D phase portrait of the system (59) in the compact phase space for two cases, showing different types of physically viable nonsingular bouncing trajectories.

Table 16 Stability and the cosmological behavior of the physically viable fixed points defined in Table 15 in presence of pressureless dust ($\omega _m=0$). Stability of the fixed points (or the line of fixed points) $B_{1\pm }$, $B_{2\pm }^{a,b,c}$ can be determined by investigating the stability of the invariant submanifolds $\bar{x}=\pm 1$ and $\bar{y}=0,1$ (see Appendix C)

The phase trajectories in the right panel of Fig. 7 shows nonsingular bouncing trajectories for the case $F(X)=\beta X^3\,(\beta <0),\,V(\phi )=V_0 e^{-\lambda \phi /M_{pl}}$. In this case, two types of bouncing trajectories are possible, e.g. trajectories connecting a contracting De-Sitter phase to an expanding De-Sitter phase ($B_{1-}\rightarrow B_{1+}$), trajectories connecting a static and a De-Sitter phase ($B_{1-}\rightarrow B_{5-}$ and $B_{5+}\rightarrow B_{1+}$). However, one cannot say this is the generic behaviour of phase trajectories for $F(X)=\beta X^3\,(\beta <0),\,V(\phi )=V_0 e^{-\lambda \phi /M_{pl}}$, as $B_{1-}$ and $B_{1+}$ are not global repellers and attractors. For $m=3$, there are another attractor/repeller pair $B_{5\pm }$. Therefore, there can exist heteroclinic trajectories connecting them, which may or may not correspond to a nonsingular bouncing cosmology, depending on how the trajectories evolve.

The phase trajectories in the left panel of Fig. 7, which corresponds to $F(X)=\beta X^{2/3}\,(\beta <0),\,V(\phi )=V_0 e^{-\lambda \phi /M_{pl}}$, show nonsingular bouncing trajectories connecting two De-Sitter phases ($B_{1-} \rightarrow B_{1+}$). For $m=\frac{2}{3}$, $B_{1+}$ and $B_{1-}$ are the only attractor/repeller pair possible, i.e., they are global attractors and repellers. The trajectories connecting them must necessarily undergo a nonsingular bounce and we have shown explicitly in the figure five such characteristic trajectories. In this case, one can confidently say that a nonsingular bounce is a generic behaviour of the phase trajectories.

Fig. 7

Phase portrait of the system (59) with shaded region represents a non-physical region of the phase space. The plots are for $m=2/3$ (left panel) and $m=3$ (right panel) moreover ${{\,\textrm{sgn}\,}}((2m-1)\beta )=-1$

4 Summary

A careful analysis of the phase space structure of the models considered reveals some important physical aspects of the models. In this section, we summarize some of the main points of our phase space analysis.

• Matter dominated epoch: At the fixed points $B_{2\pm }$, $(x,y)\rightarrow (\pm \infty ,0)$ and one can write from Eqs. (11) and (19a)

  $$\begin{aligned} \frac{\Omega _m}{x^2} = \frac{\rho _m}{3M_{Pl}^{2}H^2} \rightarrow 1. \end{aligned}$$

  (60)

  Therefore these two points are matter-dominated fixed points having a power law evolution (see Tables 4, 8, 12, 16). In all the cases considered above $B_{2\pm }$ are saddles. The saddle fixed point $B_{2+}$ can be interpreted as the intermediate matter-dominated epoch during the expanding phase of the universe.

• Genericity of nonsingular bounces: By genericity of nonsingular bouncing solutions, what we mean is, no matter what initial state one may choose during the contracting phase (i.e. whatever phase space point one chooses in the region $\bar{x}<0$), one always ends up with a nonsingular bounce. If there is a global repeller in the region $\bar{x}<0$ and a global attractor in the region $\bar{x}>0$, then all the heteroclinic trajectories must necessarily represent a nonsingular bounce and one can say that nonsingular bounce is generic. In this case, one can also say that the bouncing solutions are stable in both past and future directions. No perturbation of initial conditions alters the past and future asymptotic of the evolution.

  Nonsingular bounces are only possible for phantom scalar fields within GR. When the kinetic term in the Lagrangian of the scalar field is $F(X)=\beta X$ ($\beta <0$), nonsingular bounce is generic for an exponential potential but not for a power law potential, as we have discussed in Sects. 3.1.1 and 3.1.2 and as is also clear from Figs. 2, 3, 4. This is because the fixed points $B_{1\pm }$ are global attractors/repellers for an exponential potential but not for a power law potential.

  When the kinetic term in the Lagrangian of the scalar field is $F(X)=\beta X^m$ with $\beta ,\,m$ belonging to the range specified in (51), nonsingular bounce is not generic for either power law or exponential potentials, as we have discussed in Sects. 3.2.1 and 3.2.2 and as is also clear from the examples presented in Figs. 6, 7. This is because the fixed points $B_{1\pm }$ are not global attractors/repellers for either case. It is possible for heteroclinic trajectories to exist which do not undergo any bounce. However, there exist specific ranges for the model parameters for which nonsingular bounce is generic. For $F(X)=\beta X^m$ ($m\ne 1$) and $V(\phi )=V_0 \phi ^n$, when $\{m,\,n\}$ is within the range

  $$\begin{aligned} \frac{1}{2}< m< 1, \quad n < \frac{2m}{m-1}, \end{aligned}$$

  (61)

  the points $B_{3\pm }$ and $B_{5\pm }$ are saddles and $B_{1\pm }$ becomes global attractors and repellers (see Table 12). For $F(X)=\beta X^m$ ($m\ne 1$) and $V(\phi )=V_0 e^{-\lambda \phi /M_Pl}$, when

  $$\begin{aligned} \frac{1}{2}< m < 1, \end{aligned}$$

  (62)

  the points $B_{5\pm }$ are saddles and $B_{1\pm }$ becomes global attractors and repellers (see Table 16). In these cases, nonsingular bounce becomes generic, as can be seen in the examples presented in the left panels of Figs. 6 and 7.

• Cosmic future when the bounce is generic: It is interesting to investigate from the phase space the cosmic future of the nonsingular bouncing cosmologies when they are generic.

  For $F(X)=\beta X$ ($\beta <0$) and an exponential potential, although mathematically the fixed points $A_{2\pm }$ are non-hyperbolic, from Fig. 4 one can see that they are actually global future and past attractors. The end state of a nonsingular bouncing cosmology is a future attractor which can be either a De-Sitter phase or a big-rip cosmology, depending on the choice of the model parameters. If $\sigma =\frac{\sqrt{{2}/{3}}\lambda }{\sqrt{-\beta }}\ge 0$, the future attractor is a De-Sitter phase given by the point $B_{1+}$. In this case, the cosmology can be matched with the $\Lambda$CDM at the asymptotic future. If, however, $\sigma =\frac{\sqrt{{2}/{3}}\lambda }{\sqrt{-\beta }}<0$, then the future attractor is a big-rip singularity given by the point $A_{2+}$.

  On the contrary, when $F(X)=\beta X^m$ for either a power law or an exponential potential, and in cases when the bounce is generic, the only end state that is possible is the De-Sitter phases represented by the point $B_{1+}$. In these cases, the cosmology can be matched with $\Lambda$CDM at the asymptotic future.

5 Conclusion

Nonsingular bouncing solutions are important candidates for early universe cosmology and it is necessary to investigate different aspects of them. This article deals with investigating them from the phase space point of view. For our analysis, we have taken nonsingular bouncing models in $F(X)-V(\phi )$ theory, considering two simple choices for the potential, namely power law and exponential potential. Our main motivation behind doing this exercise is to find how generic nonsingular bouncing solutions are. More precisely, even if for a theory a bouncing solution may exist, whether or not it arises from a large number of initial conditions. In the phase space picture, this question can be recast as to whether phase trajectories representing nonsingular bouncing solutions come from a large area of the phase space or only from some small patches.

For the purpose of a dynamical system analysis we have used the formulation of references [21, 22]. The formulation employs density-normalized dimensionless dynamical variables instead of the usual Hubble-normalized dynamical variables. This is because the Hubble normalized dynamical variables diverge at the bounce. We extend the phase space analysis of [21] by compactifying the phase space. Compactification of the phase space helps us visualize its global structure and answer questions regarding past and future asymptotic of a cosmological model. In our case, the compact phase space analysis helps us investigate the genericity of solutions as well as answer questions about their past and future asymptotic dynamics.

For both potentials, we prove the existence of intermediate matter-dominated epochs which arise as saddle fixed points in the phase space. We recover the result that for a nonsingular bounce to happen in an $F(X)-V(\phi )$ type scalar field Lagrangian in GR, the scalar field needs to be phantom.¹This result is also expected out of the fact that nonsingular bounce in spatially flat FLRW cosmology in GR requires violation of NEC [2, 3]. We showed that, in general, nonsingular bounce in these models is not generic due to the non-existence of global past or future attractors. However, we identify the parameter range for which nonsingular bounce can be generic. For the kinetic term $F(X)=\beta X^m$ ($m\ne 1/2,\,\beta <0$), these ranges are as follows:

• For a power law potential $V(\phi ) = V_0 \phi ^n$, the range is $\left\{ \frac{1}{2}<m<1,\,n<\frac{2m}{m-1}\right\}$.

• For an exponential potential $V(\phi ) = V_0 e^{-\lambda \phi /M_{Pl}}$, the range is $\left\{ \frac{1}{2}<m\le 1\right\}$.

When the model parameters lie outside this specific range, obtaining a bouncing solution is still feasible by carefully selecting a particular set of initial conditions and numerically evolving the system. However, in such cases, achieving a bounce may require fine-tuning the initial conditions. Even a slight arbitrary alteration in the initial conditions may lead to a phase trajectory that does not exhibit a bounce.

Conversely, when the model parameters fall within the mentioned range, any arbitrary change in the initial conditions will still result in a bounce. In this scenario, there is no need for precise fine-tuning of the initial conditions to achieve the bouncing behavior. For the special case when $F(X) = \beta X$ ($\beta <0$) and $V(\phi ) = V_0 e^{-\lambda \phi /M_{Pl}}$, the asymptotic future of the bouncing cosmology can be either De-Sitter or a big-rip. In all the other cases when the bounce is generic, the asymptotic future is definitely De-Sitter.

Despite the challenges arising from inhomogeneous cosmological perturbations in the $F(X)-V(\phi )$ model, we have chosen to work with this simpler model to highlight the significance of exploring the genericity of bouncing solutions, specifically their stability concerning initial condition perturbations. Our investigation employs the compact phase space formulation to address this question. A case that we left out of consideration is the very special case of Cuscuton ($F(X)\propto \pm \sqrt{X}$). Phase space analysis of Cuscuton cosmology should be done with care as the Cuscuton field is nondynamical and provides an additional constraint. Since the field is nondynamical, it has been suggested that a bounce with a phantom Cuscuton does not lead to ghost instabilities [26]. It is therefore interesting to do the same exercise for a Cuscuton bounce, something which we plan to address in a future work.

We assert that, within the context of the theory being considered, cases where the bounce is generic hold greater interest than those where it is not. By examining the genericity of bounces, we gain valuable insights into the stability and robustness of these solutions against variations in initial conditions. This knowledge is crucial in understanding bouncing cosmological models’ overall behaviour and viability. Furthermore, employing a compact phase space analysis enables us to explore the cosmic fate of nonsingular bouncing cosmologies when they are considered as generic solutions. The existing literature is replete with various nonsingular bouncing models, spanning both in GR and modified gravity theories (for a comprehensive overview, refer to the reviews [3, 6]). Recent works [10, 30, 31] have also contributed to this intriguing area of research. Given that a compact phase space analysis captures the system’s global dynamics, it holds tremendous potential for studying theories that feature nonsingular bouncing solutions. It would be particularly captivating to explore the compact phase space for such theories, identifying cases where the bounce is generic and comprehending the cosmic future within these scenarios. These fascinating ideas remain as promising avenues for future projects, enriching our understanding of the broader implications of nonsingular bouncing cosmologies in the cosmos.

Appendices

Appendix A: Center manifold dynamics for point $A_{1-}$ of model $F(X)=\beta X, V(\phi )=V_0 \phi ^n$

In this appendix, we apply the center manifold theory to study the dynamics of the system (25) near a point $A_{1-}(-1,0,0)$. For detailed mathematical background on center manifold theory, we refer the reader to Ref. [32].

Firstly, we translate the point $(-1,0,0)$ to (0, 0, 0) via a transformation $x\rightarrow x-1$, $y\rightarrow y$, $\sigma \rightarrow \sigma$ under which the system (25) becomes

$$\begin{aligned} \frac{dx}{d\tau }= & {} \frac{3}{2}(x-1)\left[ 2(x-1) - \sigma y^2 \right] \nonumber \\{} & {} - \frac{3}{2}\left[ (1-\omega _{m}) + (1+\omega _{m})((x-1)^{2}-y^2 )\right] , \end{aligned}$$

(A1a)

$$\begin{aligned} \frac{dy}{d\tau }= & {} \frac{3}{2}y\left[ -\sigma + 2(x-1) - \sigma y^2 \right] ,\end{aligned}$$

(A1b)

$$\begin{aligned} \frac{d\sigma }{d\tau }= & {} \frac{3}{n}\sigma ^{2}, \end{aligned}$$

(A1c)

which can be written as

$$\begin{aligned} \left( \begin{array}{c} \frac{dx}{d\tau } \\ \frac{dy}{d\tau }\\ \frac{d\sigma }{d\tau } \end{array} \right) =\left( \begin{array}{ccc} -3(1-\omega _{m}) &{} 0 &{} 0 \\ 0 &{} -3&{} 0\\ 0 &{} 0 &{} 0 \end{array} \right) \left( \begin{array}{c} x\\ y\\ \sigma \end{array} \right) +\left( \begin{array}{c} g_1\\ g_2\\ f \end{array} \right) , \end{aligned}$$

where

$$\begin{aligned} g_1= & {} -\frac{3}{2}\,\sigma \,x{y}^{2}+\frac{3}{2}\,{x}^{2}+\frac{3}{2}\,\sigma \,{y}^{2}-\frac{3}{2}{x}^{2}\omega _{{m}}\\{} & {} \quad +\frac{3}{2}\,{y}^{2}\omega _{{m}}+\frac{3}{2}\,{y}^{2},\\ g_2= & {} -\frac{3}{2}\,\sigma \,{y}^{3}-\frac{3}{2}\sigma \,y+3\,xy,\\ f= & {} \frac{3}{n}\sigma ^{2}. \end{aligned}$$

The local center manifold is given by

$$\begin{aligned} \lbrace \textbf{z}=\textbf{h}(\sigma ): \textbf{h}(0)=\textbf{0}, D\textbf{h}(\textbf{0})=\textbf{0}\rbrace , \end{aligned}$$

(A2)

where

$$\begin{aligned} \textbf{h}=\left( \begin{array}{c} h_1\\ h_2 \end{array} \right) ,\quad \textbf{z}=\left( \begin{array}{c} x\\ y \end{array} \right) . \end{aligned}$$

The function $\textbf{h}$ can be approximated in terms of a power series as

$$\begin{aligned} h_1(\sigma )= & {} a_2 \sigma ^2+a_3 \sigma ^3+\mathcal {O}(\sigma ^4), \end{aligned}$$

(A3)

$$\begin{aligned} h_2(\sigma )= & {} b_2 \sigma ^2+b_3 \sigma ^3+\mathcal {O}(\sigma ^4), \end{aligned}$$

(A4)

which is determined by the quasilinear partial differential equation [32]

$$\begin{aligned} D \textbf{h}(\sigma )\left[ A \sigma +\textbf{F}(\sigma ,\textbf{h}(\sigma ))\right] -B \textbf{h}(\sigma )-\textbf{g}(\sigma ,\textbf{h}(\sigma ))=\textbf{0} , \end{aligned}$$

(A5)

with

$$\begin{aligned} \textbf{g}=\left( \begin{array}{c} g_1\\ g_2 \end{array} \right) ,\quad \textbf{F}=f, \quad B= \left( \begin{array}{cc} -3 (1-\omega _{m}) &{} 0 \\ 0 &{}-3 \end{array} \right) ,\quad A=0. \end{aligned}$$

In order to solve the Eq. (A5), we substitute A, ${\textbf {h}}$, $\textbf{F}$, B, $\textbf{g}$ in to it. On comparing like powers of $\sigma$ from Eq. (A5) we obtain that constants $a_2=0$, $a_3=0$, $b_2=0$, $b_3=0$. Thus, the local center manifold is the $\sigma$-axis which coincides with the center subspace (a subspace generated by the eigenvectors corresponds to a vanishing eigenvalue of the corresponding Jacobian matrix).

Finally, the dynamics in a local center manifold is determined by the equation

$$\begin{aligned} \frac{d\sigma }{d\tau } =A\,\sigma +\textbf{F}(\sigma ,\textbf{h}(\sigma )), \end{aligned}$$

(A6)

which is simply

$$\begin{aligned} \frac{d\sigma }{d\tau } =\frac{3}{n}\sigma ^2. \end{aligned}$$

(A7)

Hence, point $A_{1-}$ is always a saddle as we obtain an even-parity order term.

Appendix B: Stability at invariant sub-manifold

It is important to investigate the stability of invariant submanifolds because that helps us infer the nature of stability of the fixed points that lie at the intersection of the invariant submanifolds. If a fixed point is located at the intersection of N invariant submanifolds in an N-dimensional phase space, then the point is

• Stable if all the invariant submanifolds are of attracting nature.

• Unstable if all the invariant submanifolds are of repelling nature.

• Saddle if some of the invariant submanifolds are attracting and some are repelling.

For model II considered in Sect. 3.1.2, the fixed points $A_{1\pm }$ are at the intersection of the invariant submanifolds $y=0,\,\Omega =0$. For both the models, the fixed points (or the lines of fixed points) $B_{1\pm }$ are at the intersection of the invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=1$ and the fixed points $B_{2\pm }$ are at the intersection of the invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=0$.

Below we investigate the stability of different invariant submanifolds.

1.1 B.1 $\bar{x}=\pm 1$

Consider $\bar{x}$ in the vicinity of $+1$ or $-1$, i.e. $0<\epsilon \equiv (1-\bar{x}^{2})\ll 1$. One can then rewrite Eq. (41a) as

$$\begin{aligned} \frac{d\bar{x}}{d\bar{\tau }}= & {} \frac{3}{2}(\omega _{k}-\omega _{m})\left( 1-\epsilon -\epsilon {{\,\textrm{sgn}\,}}(\beta )\right) \epsilon (1-\bar{y}^{2}) \nonumber \\{} & {} + \frac{3}{2}\left( (1+\omega _{m})\epsilon -\sigma \bar{x}{{\,\textrm{sgn}\,}}(\beta )\sqrt{\epsilon }\right) \epsilon \bar{y}^{2}. \end{aligned}$$

(B1)

Since $\omega _{k}=1$, to the lowest order of $\epsilon$ one can write

$$\begin{aligned} \frac{d\bar{x}}{d\bar{\tau }} \simeq \frac{3}{2}(1-\omega _{m})(1-\bar{y}^{2})\epsilon . \end{aligned}$$

(B2)

Since all quantities on the right-hand side are positive. The flow will always be in positive x-direction in the vicinity of $\bar{x}=\pm 1$. Therefore the invariant submanifold $\bar{x}=1$ is attracting while the invariant submanifold $\bar{x}=-1$ is repelling.

1.2 B.2. $\bar{y}=1$

Consider $\bar{y}$ in the vicinity of 1, i.e. $0<\epsilon \equiv (1-\bar{y}^{2})\ll 1$. We will consider the cases $\sigma =0$ and $\sigma \ne 0$ separately.

1.2.1 B.1.1 $\sigma =0$

One can rewrite Eq. (41b) as

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }} = \frac{3}{2}\bar{y}(1-\bar{y}^{2})^{2}(\omega _{k}+1)\bar{x} = \frac{3}{2}(\omega _{k}+1)\bar{x}\epsilon ^{2}\sqrt{1-\epsilon }, \end{aligned}$$

(B3)

Since $\omega _k=1$, to the lowest order of $\epsilon$ one can write

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }} \simeq 3\bar{x}\epsilon ^{2}. \end{aligned}$$

(B4)

From the above one can conclude that

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is attracting for $\bar{x}>0$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is repelling for $\bar{x}<0$.

1.2.2 B.1.2 $\sigma \ne 0$

One can rewrite Eq. (41b) as

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= & {} \frac{3}{2}\bar{y}(1-\bar{y}^{2}) \left[ (-\sigma \sqrt{1-\bar{x}^{2}} +(\omega _{k}+1)\bar{x})(1-\bar{y}^{2})\right. \nonumber \\{} & {} \left. - \sigma \bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )\right] \nonumber \\= & {} \frac{3}{2}\epsilon \sqrt{1-\epsilon }\left[ (-\sigma \sqrt{1-\bar{x}^{2}}+(\omega _{k}+1) \bar{x})\epsilon \right. \nonumber \\{} & {} \left. - \sigma (1-\epsilon )\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )\right] , \end{aligned}$$

(B5)

Since $\omega _k=1$, to the lowest order of $\epsilon$ one can write

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }} \simeq -\frac{3}{2}\epsilon \sigma \sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta ). \end{aligned}$$

(B6)

From the above one can conclude that

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is attracting if ${{\,\textrm{sgn}\,}}(\beta )\ne {{\,\textrm{sgn}\,}}(\sigma )$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is repelling if ${{\,\textrm{sgn}\,}}(\beta )={{\,\textrm{sgn}\,}}(\sigma )$.

1.3 B.3 $\bar{y}=0$

Consider $\bar{y}$ in the vicinity of $\bar{y}=0$, i.e. $0<\bar{y}\ll 1$. We will consider the cases $\sigma =0$, $\sigma >0$ and $\sigma <0$ separately.

1.3.1 B.3.1 $\sigma =0$

One can rewrite Eq. (41b) as

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }} = \frac{3}{2}\bar{y}(1-\bar{y}^{2})^{2}(\omega _{k}+1)\bar{x}. \end{aligned}$$

(B7)

Since $\omega _k=1$, to leading order of $\bar{y}$ one can write

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }} \simeq 3\bar{x}\bar{y}. \end{aligned}$$

(B8)

From the above one can conclude that

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is repelling for $\bar{x}>0$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is attracting for $\bar{x}<0$.

1.3.2 B.3.2 $\sigma \ne 0$

One can rewrite Eq. (41b) as

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= & {} \frac{3}{2}\bar{y}(1-\bar{y}^{2})\left[ (-\sigma \sqrt{1-\bar{x}^{2}} +(\omega _{k}+1)\bar{x})(1-\bar{y}^{2}) \right. \nonumber \\{} & {} \left. - \sigma \bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )\right] . \end{aligned}$$

(B9)

Since $\omega _k=1$, to the leading order of $\bar{y}$ one can write

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }} \simeq \frac{3}{2}\bar{y}\left[ -\sigma \sqrt{1-\bar{x}^{2}}+2\bar{x} \right] . \end{aligned}$$

(B10)

From the above one can conclude that

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is repelling for $\bar{x}>\frac{\sigma }{\sqrt{4+\sigma ^{2}}}$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is attracting for $\bar{x}<\frac{\sigma }{\sqrt{4+\sigma ^{2}}}$.

1.4 B.3.4 $\Omega _{m}=0$

Apart from $\bar{x}=\pm 1$ and $\bar{y}=0,1$, there is another invariant submanifold given by $\Omega _m=0$. The existence of this invariant submanifold is not apparent from the dynamical system (41). One could have guessed the existence of this submanifold from the physical argument that, if cosmology is initially a vacuum, it remains so as there is no mechanism for matter creation in classical physics. That $\Omega _m=0$ is an invariant submanifold can be explicitly shown if one tries to write a dynamical equation for $\Omega _m$ in terms of $\bar{\tau }$ using the dynamical equations (19) and the definitions (31), (40)

$$\begin{aligned} \frac{d\Omega _{m}}{d\bar{\tau }}=3\Omega _{m}[(1-\omega _{m})\bar{x}(1-\bar{y}^{2})-\sigma \bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )] . \end{aligned}$$

(B11)

From the above one can conclude that the invariant submanifold $\Omega _{m}=0$ is

• attracting for $(1-\omega _{m})\bar{x}(1-\bar{y}^{2})-\sigma \bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )<0$.

• repelling for $(1-\omega _{m})\bar{x}(1-\bar{y}^{2})-\sigma \bar{y}^{2}\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}(\beta )>0$.

Appendix C: Stability at invariant sub-manifold of $F(X)=\beta X^{m}$

As we have seen the stability through the invariant submanifold, the fixed points (or the lines of fixed points) $B_{1\pm }$ are at the intersection of the invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=1$, the fixed points $B_{2\pm }^{a},\, B_{2\pm }^{b},\, B_{2\pm }^{c}$ are at the intersection of the invariant submanifolds $\bar{x}=\pm 1,\,\bar{y}=0$ and the fixed points $B_{3\pm }$ are at the intersection of the invariant submanifolds $\bar{y}=0,\,\bar{\sigma }=\pm 1$.

1.1 C.1 $\bar{x}=\pm 1$

Consider $\bar{x}$ in the vicinity of $+1$ or $-1$, i.e. $0<\epsilon \equiv (1-\bar{x}^{2})\ll 1$.

$$\begin{aligned} \frac{d\bar{x}}{d\bar{\tau }}=\frac{3}{2}\epsilon \left( \frac{1}{2m-1}-\omega _{m}\right) (1-\bar{y}^{2})\sqrt{1-\bar{\sigma }^{2}} \end{aligned}$$

(C1)

Since all the quantities $\epsilon ,(1-\bar{y}^{2}), \sqrt{1-\bar{\sigma }^{2}}$ are positive in the right hand side. If $\omega _{m}=0$, the quantity $\frac{1}{2m-1}$ will decide the attracting or repelling behaviour of the $\bar{x}=\pm 1$.

• $\frac{d\bar{x}}{d\bar{\tau }}>0$ if $m>\frac{1}{2}$ the invariant submanifold $\bar{x}=1$ is attracting and the invariant submanifold $\bar{x}=-1$ is repelling.

• $\frac{d\bar{x}}{d\bar{\tau }}<0$ if $m<\frac{1}{2}$ the invariant submanifold $\bar{x}=1$ is repelling and the invariant submanifold $\bar{x}=-1$ is attracting.

1.2 C.2 $\bar{y}= 1$

1.2.1 C.2.1 $\bar{\sigma }\ne 0$

Consider $\bar{y}$ in the vicinity of 1, i.e. $0<\epsilon \equiv (1-\bar{y}^{2})\ll 1$.

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}=-\frac{3}{2}\epsilon \bar{\sigma }\sqrt{1-\bar{x}^{2}}{{\,\textrm{sgn}\,}}((2m-1)\beta ) \end{aligned}$$

(C2)

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is attracting if ${{\,\textrm{sgn}\,}}((2m-1)\beta )\ne {{\,\textrm{sgn}\,}}(\bar{\sigma })$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is repelling if ${{\,\textrm{sgn}\,}}((2m-1)\beta )= {{\,\textrm{sgn}\,}}(\bar{\sigma })$.

1.2.2 C.2.2 $\bar{\sigma }= 0$

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= \left( \frac{3m}{2m-1}\right) \bar{x}\epsilon ^{2} \end{aligned}$$

(C3)

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is attracting if $(\bar{x}>0)\wedge \left( (m<0)\vee (m>\frac{1}{2})\right)$ or $(\bar{x}<0)\wedge (0<m<\frac{1}{2})$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=1$ i.e. the invariant submanifold $\bar{y}=1$ is repelling if $(\bar{x}<0)\wedge \left( (m<0)\vee (m>\frac{1}{2})\right)$ or $(\bar{x}>0)\wedge (0<m<\frac{1}{2})$.

1.3 C.3 $\bar{y}=0$

1.3.1 C.3.1 $\bar{\sigma }\ne 0$

Consider $\bar{y}$ in the vicinity of 0, i.e. $0<\bar{y}\ll 1$.

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= \frac{3}{2}\bar{y}\left( -\bar{\sigma }\sqrt{1-\bar{x}^{2}}+\left( \frac{2m}{2m-1}\right) \bar{x}\sqrt{1-\bar{\sigma }^{2}}\right) \end{aligned}$$

(C4)

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is repelling if $\bar{\sigma }=-1$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is attracting if $\bar{\sigma }=1$.

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is repelling if $(\bar{x}=1)\wedge \left( (m<0)\vee (m>\frac{1}{2})\right)$ or $(\bar{x}=-1)\wedge (0<m<\frac{1}{2})$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is attracting if $(\bar{x}=-1)\wedge \left( (m<0)\vee (m>\frac{1}{2})\right)$ or $(\bar{x}=1)\wedge (0<m<\frac{1}{2})$.

1.3.2 C.3.2 $\bar{\sigma }= 0$

$$\begin{aligned} \frac{d\bar{y}}{d\bar{\tau }}= \left( \frac{3m}{2m-1}\right) \bar{x}\bar{y} \end{aligned}$$

(C5)

• $\frac{d\bar{y}}{d\bar{\tau }}>0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is repelling if $(\bar{x}>0)\wedge \left( (m<0)\vee (m>\frac{1}{2})\right)$ or $(\bar{x}<0)\wedge (0<m<\frac{1}{2})$.

• $\frac{d\bar{y}}{d\bar{\tau }}<0$ in the vicinity of $\bar{y}=0$ i.e. the invariant submanifold $\bar{y}=0$ is attracting if $(\bar{x}<0)\wedge \left( (m<0)\vee (m>\frac{1}{2})\right)$ or $(\bar{x}>0)\wedge (0<m<\frac{1}{2})$.

1.4 C.3.4 $\bar{\sigma }=\pm 1$

$$\begin{aligned}{} & {} \frac{d\bar{\sigma }}{d\bar{\tau }}=-3\epsilon \sqrt{1-\bar{x}^{2}}\left( (1-\bar{y}^{2})(\Gamma -1)+\frac{(2m -3)m+1}{(4m-2)m}\bar{y}^{2}\right) \nonumber \\{} & {} \quad \frac{d\bar{\sigma }}{d\bar{\tau }}={\left\{ \begin{array}{ll} -3\epsilon \sqrt{1-\bar{x}^{2}}\left( -\frac{(1-\bar{y}^{2})}{n}+\frac{(2m -3)m+1}{(4m-2)m}\bar{y}^{2}\right) , &{} \text {For power law potential}\\ -\frac{3((2m -3)m+1)}{(4m-2)m}\epsilon \sqrt{1-\bar{x}^{2}}\bar{y}^{2}, &{} \text {For exponential law potential} \end{array}\right. } \end{aligned}$$

(C6)

For the power law potential

1.4.1 C.3.1 $\bar{y}\ne 0$

• $\frac{d\bar{\sigma }}{d\bar{\tau }}>0$ if $-\frac{1}{n}+\left( \frac{(2m -3)m+1}{(4m-2)m}+\frac{1}{n}\right) \bar{y}^{2}<0$ the invariant submanifold $\bar{\sigma }=-1$ is repelling and $\bar{\sigma }=1$ is attracting.

• $\frac{d\bar{\sigma }}{d\bar{\tau }}<0$ if $-\frac{1}{n}+\left( \frac{(2m -3)m+1}{(4m-2)m}+\frac{1}{n}\right) \bar{y}^{2}>0$ the invariant submanifold $\bar{\sigma }=-1$ is attracting and $\bar{\sigma }=1$ is repelling.

1.4.2 C.3.2 $\bar{y}= 0$

• $\frac{d\bar{\sigma }}{d\bar{\tau }}>0$ if $\frac{1}{n}>0$ the invariant submanifold $\bar{\sigma }=-1$ is repelling and $\bar{\sigma }=1$ is attracting.

• $\frac{d\bar{\sigma }}{d\bar{\tau }}<0$ if $\frac{1}{n}<0$ the invariant submanifold $\bar{\sigma }=-1$ is attracting and $\bar{\sigma }=1$ is repelling.

For the exponential law potential

• $\frac{d\bar{\sigma }}{d\bar{\tau }}>0$ if $\frac{(2m -3)m+1}{(4m-2)m}<0$ i.e. $(0<m<1)$ the invariant submanifold $\bar{\sigma }=-1$ is repeling and $\bar{\sigma }=1$ is attracting.

• $\frac{d\bar{\sigma }}{d\bar{\tau }}<0$ if $\frac{(2m -3)m+1}{(4m-2)m}>0$ i.e. $(m<0)\vee (m>1)$ the invariant submanifold $\bar{\sigma }=-1$ is attracting and $\bar{\sigma }=1$ is repelling.