# A model-independent compact dynamical system formulation for exploring
bounce and cyclic cosmological evolutions in $f(R)$ gravity

A. S. Agrawal  0000-0003-4976-8769<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani, Hyderabad
Campus, Hyderabad-500078, India.  Charlotte Louw  0009-0009-7835-9631
<EMAIL_ADDRESS>Department of Mathematics and Applied Mathematics,
Cosmology and Gravity Group, University of Cape Town, Rondebosch, 7701, Cape
Town, South Africa  Saikat Chakraborty  0000-0002-5472-304X<EMAIL_ADDRESS>The Institute for Fundamental Study “The Tah Poe Academia Institute”, Naresuan
University, Phitsanulok 65000, Thailand  Center for Space Research, North-West
University, Mahikeng 2745, South Africa  Peter K.S. Dunsby
0000-0002-6271-9585<EMAIL_ADDRESS>Department of Mathematics and
Applied Mathematics, Cosmology and Gravity Group, University of Cape Town,
Rondebosch, 7701, Cape Town, South Africa  Center for Space Research, North-
West University, Mahikeng 2745, South Africa  South African Astronomical
Observatory, Observatory 7925, Cape Town, South Africa

###### Abstract

Using the dynamical systems approach together with the cosmographic
parameters, we present a _model-independent_ dynamical system formulation for
cosmology in $f(R)$ gravity. The formulation is model-independent in the sense
that one needs to specify not a particular functional form of $f(R)$ a-priori,
but rather a particular cosmological evolution, which fixes the cosmography.
In a sense, our approach is the way-around the reconstruction method. This is
shown using both non-compact and compact dynamical variables. The focus in
this paper is on the compact analysis, since we demonstrate the applicability
of this formulation using examples of bouncing and cyclic cosmology.

Keywords: Dynamical System, Bouncing Scenario, Scalar Field, $f(R)$ Gravity.

##  I Introduction

The standard cosmological model is built upon two key assumptions. Firstly, we
assume that General Relativity correctly describes the interaction between the
geometry and matter of the universe. Secondly, we assume that on large scales,
the universe is homogeneous and isotropic, which is known as the Cosmological
principle. While this may seem counter-intuitive, based on our living
experience of clumped matter, i.e., inhomogeneity on small scales, we choose
to work on a homogeneous and isotropic background and introduce perturbations
to account for the large-scale structure that we observe. In order to
describe, geometrically, this kind of background, we make use of the
Friedmann-Lemaître-Robertson-Walker (FLRW) metric. The primary evidence for
the Cosmological Principle is the observed uniformity of the temperature of
the Cosmic-Microwave Background (CMB)[1]. Thus if we adopt the homogeneity and
isotropy as an appropriate description of our universe, we then need to
develop our understanding of the gravitational interaction in order to
describe both the background expansion history and the evolution of large-
scale structure.

General Relativity (GR) has been extremely successful at correctly describing
the gravitational interaction on scales ranging from the solar system to
extragalactic scales. This is supported by a growing database of observations.
However, if one accepts the geometry presented by the FLRW metric, then
observational evidence seems to suggest that our universe is undergoing
accelerated expansion in the present epoch, thus contradicting the attractive
nature of gravity in the standard GR model. There seem to be two lines of
thought as to why this is the case. The first proposition involves the
existence of ”Dark Energy”, a mysterious matter component which currently
dominates the energy density of the universe. This component appears to affect
gravity in such a way that it becomes repulsive in its presence. While this
possibility keeps the GR description of gravity, it leads to the problem of
finding a physical description of Dark Energy. A popular candidate is that of
the cosmological constant $\Lambda$, which, when added to the Einstein-Hilbert
action, presents itself with a constant energy density and equation of state
parameter $\omega=-1$. The famous $\Lambda$CDM model of the universe is made
up of $\Lambda$ as well as another mysterious component, Cold Dark Matter
(CDM)[1]. While it has been shown that this model agrees mostly with
observations, it suffers from some problems, particularly the fine-tuning of
the energy density of $\Lambda$. The quantum vacuum energy is the only known
possibility for describing $\Lambda$, however, this value calculated in
quantum field theory differs from the observed value by more than $120$ orders
of magnitude.

Thus, we consider the second line of thought attempting to explain our
universe’s accelerated expansion. While the first approach involved altering
the right-hand side of Einstein’s equation $G_{\mu\nu}=\kappa T_{\mu\nu}$, one
may also consider altering the left side i.e., the acceleration of the
universe could be a sign of a break down of the description of the
gravitational action according to GR [2] on cosmological scales. Modified
theories of gravity ask the question of whether our understanding of how
gravity behaves at large scales is correct. The majority of these models
involve modifying the Einstein-Hilbert Lagrangian by introducing, in general,
two scalar degrees of freedom in the following form:

$R\to f(R,\phi).$

These degrees of freedom present themselves as an explicit scalar field $\phi$
as well as a hidden propagating degree of freedom given by
$F(R,\phi)=\frac{df}{dR}$ [3], while the standard matter Lagrangian remains
unchanged.

Of particular interest (both in general relativity and its various extensions)
is whether the Universe originated from a singularity at the beginning or
began expanding after a non-singular bounce-like phase. Since the inflationary
scenario can generate a power spectrum that is almost scale-invariant and in
line with the most recent observations, it is an appealing early Universe
scenario. The bouncing cosmology idea is an attractive option for inflation.
In the case of the bouncing scenario, the universe begins in a contracting
period, bounces and when the scale factor acquires minimum size, and then
begins to expand again. As a result, the bounce happens when $H=0$ and
$\dot{H}>0$. In some bouncing cosmology scenarios, in addition to producing an
observationally compatible power spectrum, having an expansion history free of
an initial singularity is another benefit of such a model.

The modified Raychaudhuri equation, the modified Friedmann equation, and the
energy conservation equation for the cosmic field make up the field equations
for $f(R)$ gravity. In general, they make up a complicated collection of
coupled non-linear partial differential equations. To retain the cosmological
principle, we primarily use homogeneous metrics in cosmology, which causes the
cosmological field equations to collapse into a collection of coupled
nonlinear ordinary differential equations. In order to extract quantitative as
well as qualitative information about the evolution of such nonlinear systems,
the dynamical system approach can be a useful strategy for studying these
equations. The cosmological field equations can be written as a closed set of
coupled first-order nonlinear differential equations, which constitutes the
dynamical system, by appropriately defining a set of dimensionless normalised
dynamical variables in the dynamical system formulation of cosmology.
Dynamical systems can be divided into two categories: continuous dynamical
systems, whose evolution is described by a collection of ordinary differential
equations and time-discrete dynamical systems, whose evolution is described by
a map or difference equations.

In this paper, we generalise a method first introduced in [3] in which the
reconstruction methods are completely avoided by making use of a dynamical
systems approach that does not require the knowledge of the functional form of
the $f(R)$ model. The dynamical systems approach has the advantage of
providing us with a qualitative understanding of cosmological models and have
seen much success in cosmology (see [4] for an application to the $\Lambda$
CDM model). Almost all past autonomous (not explicitly time-dependent)
dynamical system approaches to $f(R)$ models in literature have required the
functional form of $f(R)$ to be specified in order to close the system. One
can avoid this by introducing a set of dimensionless cosmographic parameters
which relate to the kinematical description of the universe [5]. These
parameters are obtained by taking a Taylor expansion of the scale factor
around the present redshift. As was shown in [3], if a kinematical description
of the universe is specified, then these parameters define a new constraint
which allows the dynamical system to be closed. Therefore, one can solve the
system, find the fixed points and study the phase space dynamics without
knowledge of the functional form of the $f(R)$ gravity and in fact, use this
method to explicitly determine the form of $f(R)$ from the cosmological
dynamics.

This cosmography approach is applied to several expansion histories that admit
cyclic and bouncing behaviour. By fixing the expansion history, one can easily
determine the relationship between the higher-order cosmographic parameters,
allowing one to close the system of dynamical equations in a way which is
independent of the form of $f(R)$.

This paper begins by deriving the field equations for a general $f(R)$ gravity
under the usual isotropic and homogeneous FLRW background with a perfect
fluid. In Section III, we present the model-independent dynamical systems
formulation and compactify the phasespace in Sections IV and V In Section VI,
we apply our formalism to three examples where the expansion history is fixed
by choosing the explicit form of the scale factor. This allows one to
determine the relation between the deceleration and jerk parameters and thus
close the system. Finally, in Section VII, we present our conclusions and
comment on possible future work.

##  II Field equations and the dynamical system for $f(R)$ gravity

$f(R)$ gravity can be derived from the following action [6, 7]:

$S=\frac{1}{2\kappa}\int\sqrt{-g}d^{4}x[f(R)+2\mathcal{L}_{m}]\,,$ (1)

where $\kappa=8\pi G$, $G$ being Newton’s gravitational constant, $g$ is the
metric determinant, $f(R)$ is a function of the Ricci scalar $R$,
$\mathcal{L}_{m}$ is the matter action. Varying the action with respect to the
metric yields the field equations

$G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\frac{\kappa
T_{\mu\nu}}{f^{\prime}(R)}+g_{\mu\nu}\frac{\left[f(R)-Rf^{\prime}(R)\right]}{2f^{\prime}(R)}+\frac{\nabla_{\mu}\nabla_{\nu}f^{\prime}(R)-g_{\mu\nu}\square
f^{\prime}(R)}{f^{\prime}(R)}\,,$ (2)

where the prime denotes a derivative with respect to $R$, $T_{\mu\nu}$ is the
energy momentum tensor. $f(R)$ theories of gravity propagate an extra scalar
degree of freedom in the gravity sector, sometimes dubbed as a ‘curvaton’, as
made clear from the trace field equation

$RF(R)-2f(R)+3\square F(R)=\kappa T\,,$ (3)

where $F(R)=f^{\prime}(R)$. $f(R)=R$ is the trivial case for which the action
reduces to GR, in which case the above equation becomes algebraic. Thus, there
is no extra propagating scalar degree of freedom. The absence of ghost and
tachyonic instability requires $f^{\prime}(R)=F(R)>0$ and
$f^{\prime\prime}(R)=F^{\prime}(R)\geq 0$ respectively.

For a Friedmann-Lemaître-Robertson Walker (FLRW) metric

$ds^{2}=-dt^{2}+a^{2}(t)\left[\frac{dr^{2}}{\left(1-kr^{2}\right)^{2}}+r^{2}d\theta^{2}+r^{2}sin^{2}\theta
d\phi^{2}\right]\,,$ (4)

and in the presence of a perfect fluid $T^{\mu}_{\nu}=(-\rho,P,P,P)$, the
modified Friedmann and Raychaudhuri equations can be neatly expressed as

$\displaystyle
3F\left(H^{2}+\frac{k}{a^{2}}\right)=\kappa\rho_{eff}=\kappa(\rho+\rho_{R})\,,$
(5a) $\displaystyle-F\left(2\dot{H}+3H^{2}+\frac{k}{a^{2}}\right)=\kappa
P_{eff}=\kappa(P+P_{R})\,,$ (5b)

where $\rho,\,P$ are the fluid energy density and pressure, and the curvaton
energy density and pressure are

$\displaystyle\kappa\rho_{R}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(RF-f)-3H\dot{F}$ (6a) $\displaystyle\kappa P_{R}$
$\displaystyle=$ $\displaystyle\ddot{F}+2H\dot{F}-\frac{1}{2}(RF-f)$ (6b)

Over dot represents the differentiation with respect to cosmic time throughout
the paper. The perfect fluid follows the standard continuity equation

$\dot{\rho}=-3H(1+\omega)\rho\,.$ (7)

The effective equation of state parameter and the curvaton equation of state
parameter are defined, respectively, as

$\displaystyle\omega_{eff}=\frac{P_{eff}}{\rho_{eff}}=-\frac{2\dot{H}+3H^{2}+\frac{k}{a^{2}}}{3H^{2}+\frac{k}{a^{2}}}\,,$
(8a)
$\displaystyle\omega_{R}=\frac{P_{R}}{\rho_{R}}=\frac{\ddot{F}+2H\dot{F}-\frac{1}{2}(RF-f)}{\frac{1}{2}(RF-f)-3H\dot{F}}\,.$
(8b)

The usual dynamical system formulation for cosmology in $f(R)$ gravity in
terms of expansion-normalized (or Hubble-normalized) dynamical variables was
introduced in [8, 9, 10] 111For a slightly different formulation, also in
terms of Hubble-normalized dynamical variables, see [11]. The formulation
starts with the introduction of dimensionless dynamical variables

$\displaystyle
x=\frac{\dot{F}}{HF},~{}~{}~{}y=\frac{R}{6H^{2}},~{}~{}~{}z=\frac{f}{6FH^{2}},~{}~{}~{}\Omega=\frac{\kappa\rho}{3FH^{2}},~{}~{}~{}K=\frac{k}{a^{2}H^{2}}\,.$
(9)

They are obtained by dividing the modified Friedmann equation (5a) by
$3H^{2}$, and satisfy the constraint

$-x+y-z+\Omega-K=1\,.$ (10)

Introducing the phase space-time variable $\tau=\epsilon\ln a$
($\epsilon=H/|H|=\pm 1$ for expanding/contracting universe) and the auxiliary
quantity

$\Gamma(R)\equiv\frac{d\ln R}{d\ln F}\equiv\frac{F}{RF^{\prime}}\,,$ (11)

the dynamical system can be written as

$\displaystyle\frac{dx}{\epsilon d\tau}$ $\displaystyle=$
$\displaystyle-4z-2x^{2}-(z+2)x+2y+\Omega(x+1-3\omega)\,,$ (12a)
$\displaystyle\frac{dy}{\epsilon d\tau}$ $\displaystyle=$ $\displaystyle
y[2\Omega-2(z-1)+x(\Gamma-2)]\,,$ (12b) $\displaystyle\frac{dz}{\epsilon
d\tau}$ $\displaystyle=$ $\displaystyle z(-2z+2\Omega-3x+2)+xy\Gamma\,,$ (12c)
$\displaystyle\frac{d\Omega}{\epsilon d\tau}$ $\displaystyle=$
$\displaystyle-\Omega(3x+2(z-\Omega)+(3\omega+1))\,,$ (12d)

where $K$ has been eliminated using the constraint (10). For a given
functional form of $f(R)$, one can write the relation

$\frac{y}{z}=\frac{RF}{f}\,.$ (13)

To make the dynamical system autonomous, it is necessary to invert the above
relation to find $R=R(y/z)$, provided, of course, that it is invertible. When
this is possible, one can close the dynamical system by writing
$\Gamma(R)=\Gamma(R(y/z))$. In practice, this is possible only for very
specific forms of $f(R)$.

The cosmological solutions with vanishing Ricci scalar form an invariant
submanifold, $y=0$, which divides the entire phase space into two disjoint
sectors containing solutions with positive and negative Ricci scalar. This
fact will be useful in Section IV, where we will present the compact dynamical
system formulation.

Note that the phase space-time variable is defined such that $d\tau=|H|dt$
rather than $d\tau=Hdt$. This is done to make the formulation valid for
applications in both expanding and contracting universes. The definition of
the phase space-time variable for an expanding and contracting universe is
different. In particular, $\tau$ is not continuous across a bounce or re-
collapse. This flaw is rectified in the compact dynamical system formulation
presented in Section IV.

##  III A model-independent dynamical system formulation

A _model-independent_ dynamical system formulation was proposed in Ref.[12],
which does not require one to specify a-priori an $f(R)$ model, i.e., a
functional form of $f(R)$, while still allowing one to obtain a closed
autonomous system. The formulation utilizes the cosmographic deceleration and
jerk parameters

$q\equiv-\frac{1}{aH^{2}}\frac{d^{2}a}{dt^{2}}\,,\qquad
j\equiv\frac{1}{aH^{3}}\frac{d^{3}a}{dt^{3}}\,,\qquad
s\equiv\frac{1}{aH^{4}}\frac{d^{4}a}{dt^{4}}\,,$ (14)

which are related to each other by

$j=2q^{2}+q-\frac{dq}{d\tau}\,,\qquad s=\frac{dj}{d\tau}-j(2+3q)\,.$ (15)

In terms of the same dynamical variables, as defined in Eq. (9), the Friedmann
equation and the definition of the Ricci scalar provide two constraints

$\displaystyle y$ $\displaystyle=$ $\displaystyle 1-q+K,$ (16a) $\displaystyle
z$ $\displaystyle=$ $\displaystyle\Omega-x-q$ (16b)

The resulting dynamical system for barotropic perfect fluid of equation of
state parameter $\omega$ is

$\displaystyle\frac{dx}{\epsilon d\tau}$ $\displaystyle=$ $\displaystyle
x(q-x)+2(x+K+q)-3\Omega(1+\omega)+2\,,$ (17a)
$\displaystyle\frac{d\Omega}{\epsilon d\tau}$ $\displaystyle=$
$\displaystyle-\Omega(x-2q+1+3\omega)\,,$ (17b)
$\displaystyle\frac{dK}{\epsilon d\tau}$ $\displaystyle=$ $\displaystyle
2qK\,,$ (17c) $\displaystyle\frac{dq}{\epsilon d\tau}$ $\displaystyle=$
$\displaystyle 2q^{2}+q-j\,,$ (17d) $\displaystyle\frac{dj}{\epsilon d\tau}$
$\displaystyle=$ $\displaystyle j(2+3q)+s\,,$ (17e)

where $y,\,z$ have been eliminated by the constraints (16). In principle, one
can obtain an infinite number of equations since there is an infinite
hierarchy of cosmographic parameters. The catch here is that this particular
formulation can be used to study the phase space of all possible $f(R)$
theories that can reproduce a given cosmic evolution. A given cosmic evolution
can always be specified by an algebraic relation between a finite number of
cosmographic parameters. For example, the cosmic evolution for the
$\Lambda$CDM model in the generic non-flat case is specified by the
cosmographic relation [13]

$s+2(q+j)+qj=1$ (18)

The dynamical system (17) reveals the existence of another invariant
submanifold $K=0$, containing spatially flat solutions. The existence of this
invariant submanifold was not apparent from the dynamical system (12) because
we chose to eliminate the variable $K$ using the Friedmann constraint. This
submanifold divides the entire phase space into two disjoint regions,
containing spatially positively curved and spatially negatively curved
solutions. Present-day observations marginally disfavor spatially negatively
curved solutions (see, for example, [14]), so we can safely assume that the
physically interesting cosmological dynamics take place in the $K\geq 0$
region.

For a cosmological solution to be physically realistic, it should be close to
the present-day observable universe at $z=0$ ($z$ being the redshift). At
$z=0$, $K\gtrsim 0,\,q\approx-0.55$, giving a positive value for $y$. Since
$y<0$ solutions can never cross into the $y>0$ region and vice versa, we can
safely assume that the physically viable cosmological solutions lie in the
region $y>0$. Even for non-singular bouncing cosmologies, which we will
investigate later, $q\rightarrow-\infty$ at the bounce gives $y>0$. From the
above discussions, we can safely say that the physically interesting
cosmological dynamics happen within the disjoint region $y\geq 0,\,K\geq 0$.
This fact will be useful in the formulation of the compact phase space in the
next section.

It can be noted that

$xy\Gamma=6(x+z-\Omega)+j+3q-2y$ (19)

Eliminating $y$ and $z$ on the right-hand side by using the constraint (16a)
and (16b) and using the definition of $\Gamma$ one can write

$\frac{1}{y\Gamma}=\frac{6F^{\prime}H^{2}}{F}=\frac{x}{j-q-2K-2}\;.$ (20)

We recall that the absence of ghost and tachyonic instability requires
$F>0,\,F^{\prime}\geq 0$. Assuming that the condition $F>0$ is met, demanding
$F^{\prime}\geq 0$ puts the following constraints on the phase space

$\frac{x}{j-q-2K-2}\geq 0\,.$ (21)

Therefore, not the entire phase space is physically viable. The physically
viable region of the phase space, where the theory is free from ghost and
tachyonic instability, must _necessarily_ satisfy the constraint (21). One
must be careful, though; this is _not_ a sufficient condition. It is certainly
possible that a region of the phase space satisfying the condition (21) is
plagued by both ghost and tachyonic instability ($F<0$ and $F^{\prime}<0$). To
be sure, one must explicitly calculate $F,\,F^{\prime}$ along a phase
trajectory of interest. In a region of the phase space, if the condition (21)
is satisfied and also the dynamical variable
$\Omega\equiv\frac{\kappa\rho}{3FH^{2}}\geq 0$, then one can definitely say
that the region is physically viable. We will see that this will be important
in example III of section VI.

##  IV Compact Dynamical System Formulation

In this section, following [15, 16, 17, 18, 19], we present a compact
dynamical systems formulation for $f(R)$ gravity. In this approach, all the
dimensionless dynamical variables are by definition compact, giving rise to a
compact phase space. Such a compact phase space formulation allows us to
obtain a global picture of the dynamics, allowing us to investigate fixed
points at infinity. Moreover, the formulation is made in such a way that it
can be used to investigate the phenomena of non-singular bounce and re-
collapse from the phase space point of view.

The compact phase space formulation presented here is valid for solutions with
non-negative global spatial curvature ($k\geq 0$) and non-negative Ricci
scalar ($R\geq 0$). This does not pose any limitation from the physical point
of view because, as we have already discussed in the previous sections, this
is a disjoint region in the phase space where all the physically interesting
cosmological dynamics take place. Two other crucial assumptions that are
required to make this formulation work are $f\geq 0$ and the ghost instability
never arises ($F>0$). In light of the above requirements, this formulation is
very suitable to investigate the phase space of theories like $R+\alpha R^{n}$
($\alpha>0$) [16] or the Hu-Sawicki model [18].

The compact phase space formulation starts with rewriting the Friedmann
equation as

$\left(3H+\frac{3}{2}\frac{\dot{F}}{F}\right)^{2}+\frac{3}{2}\left(\frac{f}{F}+\frac{6k}{a^{2}}\right)=\frac{3\rho}{F}+\frac{3R}{2}+\left(\frac{3}{2}\frac{\dot{F}}{F}\right)^{2}\,,$
(22)

allowing for the introduction of a new normalization factor $D$ defined as

$D^{2}=\left(3H+\frac{3}{2}\frac{\dot{F}}{F}\right)^{2}+\frac{3}{2}\left(\frac{f}{F}+\frac{6k}{a^{2}}\right)\,.$
(23)

Note that, for $k>0$, the quantity $D$ is positive definite $D>0$. A new set
of dynamical variables are then defined as,

$\bar{x}=\frac{3}{2}\frac{\dot{F}}{F}\frac{1}{D},~{}~{}~{}\bar{y}=\frac{3}{2}\frac{R}{D^{2}},~{}~{}~{}\bar{z}=\frac{3}{2}\frac{f}{F}\frac{1}{D^{2}},~{}~{}~{}\bar{Q}=\frac{3H}{D},~{}~{}~{}\bar{\Omega}=\frac{3\rho}{F}\frac{1}{D^{2}},~{}~{}~{}\bar{K}=\frac{9k}{a^{2}}\frac{1}{D^{2}}.$
(24)

We now have two independent constraints equations, one derived from the
Friedmann equation (22) and the other just from the definition of $D$ in (23),

$\displaystyle\bar{\Omega}+\bar{y}+\bar{x}^{2}$ $\displaystyle=$
$\displaystyle 1\,,$ (25a)
$\displaystyle(\bar{Q}+\bar{x})^{2}+\bar{z}+\bar{K}$ $\displaystyle=$
$\displaystyle 1\,.$ (25b)

The above constraints, together with the assumptions $k\geq 0,\,R\geq
0,\,f\geq 0,\,F>0$, compactify the phase space

$\displaystyle\bar{x}\in[-1,1],~{}~{}~{}~{}~{}~{}~{}~{}\bar{Q}\in[-2,2]~{}~{}~{}~{}~{}~{}~{}\bar{y},~{}\bar{z},~{}\bar{\Omega},~{}\bar{K}\in[0,1]$

A new phase space-time variable $\bar{\tau}$ is also defined as

$d\bar{\tau}=Ddt\,.$ (27)

The compact and the usual non-compact dynamical variables and time variables
are related as follows

$x=\frac{2\bar{x}}{\bar{Q}},~{}~{}y=\frac{\bar{y}}{\bar{Q}^{2}},~{}~{}z=\frac{\bar{z}}{\bar{Q}^{2}},~{}~{}K=\frac{\bar{K}}{\bar{Q}^{2}},~{}~{}\Omega=\frac{\bar{\Omega}}{\bar{Q}^{2}},~{}~{}\epsilon
d\tau=Hdt=\frac{\bar{Q}}{3}d\bar{\tau}.$ (28)

Using the constraint equations, we can eliminate $\bar{z},\,\bar{\Omega}$ to
arrive at the following autonomous system

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}=$
$\displaystyle\frac{1}{6}\Big{[}-2\bar{x}^{2}\bar{y}\Gamma+(1-3\omega)(1-\bar{x}^{2}-\bar{y})+2\bar{y}+4\left(\bar{x}^{2}-1\right)\left(1-\bar{Q}^{2}-\bar{x}\bar{Q}\right)+\bar{x}(\bar{Q}+\bar{x})\Big{]}$
(29a) $\displaystyle\frac{d\bar{y}}{d\bar{\tau}}=$
$\displaystyle-\frac{1}{3}\bar{y}\Big{[}(\bar{Q}+\bar{x})\Big{(}2\bar{y}-(1+3\omega)(1-\bar{x}^{2}-\bar{y})+4\bar{x}\bar{Q}\Big{)}-2\bar{Q}-4\bar{x}+2\bar{x}\Gamma(\bar{y}-1)+4\bar{x}\bar{K}$
(29b)
$\displaystyle\left((1+3\omega)(1-\bar{x}^{2}-\bar{y})-2\bar{y}\right)+4\bar{K}(1-\bar{x}^{2})\Big{]}\,,$
$\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}=$
$\displaystyle\frac{1}{6}\Big{[}-4\bar{x}\bar{Q}^{3}+\bar{x}\bar{Q}(5+3\omega)(1-\bar{x}\bar{Q})-\bar{Q}^{2}(1-3\omega)-\bar{Q}\bar{x}^{3}(1+3\omega)-3\bar{y}\bar{Q}(1+\omega)(\bar{Q}+\bar{x})$
(29c)
$\displaystyle+2\bar{y}(1-\Gamma\bar{Q}\bar{x})-2\bar{K}(1+2\bar{x}\bar{Q})\Big{]}\,,$
$\displaystyle\frac{d\bar{K}}{d\bar{\tau}}=$
$\displaystyle-\frac{1}{3}\bar{K}\Big{[}(\bar{Q}+\bar{x})\Big{(}-(1+3\omega)(1-\bar{x}^{2}-\bar{y})+4\bar{x}\bar{Q}+2\bar{y}\Big{)}+4\bar{x}(\bar{K}-1)+2\bar{x}\bar{y}\Gamma\Big{]}\,,$
(29d)

where $\Gamma\equiv\frac{F}{RF^{\prime}}$ is the same auxiliary quantity as
defined in Eq.(11). Just like the usual dynamical system formulation of
Section II, the applicability of this compact formulation is also limited by
the requirement that it should be possible to express $\Gamma$ in terms of
$\bar{y}$ and $\bar{z}$, thus closing the system (29).

##  V A model-independent compact dynamical system formulation

Inspired by the model-independent dynamical system formulation presented in
Section III, we now try to construct a similar model-independent formulation,
but in terms of the compact dynamical variables. Such a formulation will allow
us to investigate the global phase space structure of those $f(R)$ theories
that satisfy a particular cosmographic constraint without requiring us to
reconstruct their exact mathematical form, which might be too complicated for
any practical purpose. The only restriction on the mathematical form of the
$f(R)$ theories are the ones required for the compact formulation to work,
namely $f\geq 0,\,F>0$. We will again make use of the cosmographic
deceleration and jerk parameter.

In constructing the compact dynamical system (29), we chose to eliminate
$\bar{\Omega},\,\bar{z}$. If we instead choose to eliminate
$\bar{y},\,\bar{z}$, then we arrive at the following dynamical system:

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+2\bar{K}+\bar{x}\bar{Q}-2\bar{x}^{2}\bar{K}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)-\bar{\Omega}(1+3\omega)+3\bar{x}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})$
(30a)
$\displaystyle-2\bar{x}^{2}\left(1-\bar{\Omega}-\bar{x}^{2}\right)\Gamma\Big{]},$
$\displaystyle\frac{d\bar{\Omega}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{\Omega}}{3}\Big{[}3(\bar{x}\bar{\Omega}+\bar{Q}\bar{\Omega}-\bar{Q})(1+\omega)-2\bar{x}\left(1-\bar{\Omega}-\bar{x}^{2}\right)\Gamma-4\bar{x}\bar{Q}^{2}-4\bar{x}\bar{K}-2\bar{x}^{2}\bar{Q}+2\bar{x}^{3}\Big{]},$
(30b) $\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}3\bar{Q}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-2\bar{K}-4\bar{Q}^{2}+2\bar{x}\bar{Q}-2(\bar{Q}\bar{x}\Gamma-1)\left(1-\bar{\Omega}-\bar{x}^{2}\right)$
(30c)
$\displaystyle-4\bar{x}\bar{Q}\bar{K}-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}-4\bar{x}\bar{Q}^{3}\Big{]},$
$\displaystyle\frac{d\bar{K}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{K}}{6}\Big{[}4\bar{x}-4\bar{Q}-8\bar{x}\bar{K}-8\bar{x}\bar{Q}^{2}-4\bar{x}^{2}\bar{Q}+4\bar{x}^{3}+6\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-4\bar{x}\left(1-\bar{\Omega}-\bar{x}^{2}\right)\Gamma\Big{]}.$
(30d)

The auxiliary variable $\Gamma$, which encodes the information about the
explicit form of $f(R)$, enters the dynamical system only through the
combination

$\bar{x}(1-\bar{\Omega}-\bar{x}^{2})\Gamma=\bar{x}\bar{y}\Gamma\,.$ (31)

Using the mapping (28), the relations (19) and the constraints (25), it is
possible to express the above expression in terms of the cosmographic
deceleration and jerk parameter

$\bar{x}\left(1-\bar{\Omega}-\bar{x}^{2}\right)\Gamma=\frac{\bar{Q}^{3}}{2}(j-q-2)-\bar{K}\bar{Q}\,.$
(32)

Substituting this back into the system (30) gives the dynamical system

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+2\bar{K}+\bar{x}\bar{Q}-2\bar{x}^{2}\bar{K}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)-\bar{\Omega}(1+3\omega)+3\bar{x}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})$
(33a)
$\displaystyle-\bar{x}\bar{Q}^{3}(j-q-2)+2\bar{x}\bar{K}\bar{Q}\Big{]}\,,$
$\displaystyle\frac{d\bar{\Omega}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{\Omega}}{3}\Big{[}3(\bar{x}\bar{\Omega}+\bar{Q}\bar{\Omega}-\bar{Q})(1+\omega)-\bar{Q}^{3}(j-q-2)+2\bar{K}\bar{Q}-4\bar{x}\bar{Q}^{2}-4\bar{x}\bar{K}-2\bar{x}^{2}\bar{Q}+2\bar{x}^{3}\Big{]}\,,$
(33b) $\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}3\bar{Q}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-2\bar{K}-4\bar{Q}^{2}+2\bar{x}\bar{Q}-\bar{Q}^{4}(j-q-2)+2\bar{K}\bar{Q}^{2}+2\left(1-\bar{\Omega}-\bar{x}^{2}\right)$
(33c)
$\displaystyle-4\bar{x}\bar{Q}\bar{K}-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}-4\bar{x}\bar{Q}^{3}\Big{]}\,,$
$\displaystyle\frac{d\bar{K}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{K}}{6}\Big{[}4\bar{x}-4\bar{Q}-8\bar{x}\bar{K}-8\bar{x}\bar{Q}^{2}-4\bar{x}^{2}\bar{Q}+4\bar{x}^{3}+6\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-2\bar{Q}^{3}(j-q-2)+4\bar{K}\bar{Q}\Big{]}\,,$
(33d)

which contains no $\Gamma$-term, but instead contains the cosmographic
parameters $q,\,j$. The last piece of the trick is to use the definition of
the Ricci scalar in Eq.(16a), the mapping (28) and the Friedmann constraint
(25a) in such a way as to be able to write

$\displaystyle
q=1-\frac{\bar{y}-\bar{K}}{\bar{Q}^{2}}=1-\frac{1-\bar{\Omega}-\bar{x}^{2}-\bar{K}}{\bar{Q}^{2}}\,.$
(34)

If a cosmic evolution can be specified by an algebraic relation of the form
$j=j(q)$, then the dynamical system (33), together with the relation (34),
provides a model-independent dynamical system, in the sense that it can be
used to investigate the phase space of all such $f(R)$ theories that support
the cosmology $j=j(q)$.

We notice the existence of two invariant submanifolds corresponding to the
spatially flat case ($\bar{K}=0$) and the vacuum case ($\bar{\Omega}=0$). At
the intersection of these two invariant submanifolds, which is an invariant
submanifold itself corresponding to the spatially flat vacuum case, the flow
is 2-dimensional.

The physically viable region of the phase space, where the theory is free from
ghost and tachyonic instability, is constrained by the conditions
$F>0,\,F^{\prime}\geq 0$. Noting that

$\frac{2}{3}\frac{F^{\prime}D^{2}}{F}=\frac{1}{\left(1-\bar{\Omega}-\bar{x}^{2}\right)\Gamma}=\frac{2\bar{x}}{\bar{Q}^{3}(j-q-2)-2\bar{K}\bar{Q}}\,,$
(35)

the physically viable region of the phase space must necessarily satisfy the
condition

$\frac{2\bar{x}}{\bar{Q}^{3}(j-q-2)-2\bar{K}\bar{Q}}\geq 0\,.$ (36)

The above condition is equivalent to the condition (21). We recall that the
condition (21) was only a necessary condition but not a sufficient condition
for the physical viability of a region of the phase space. This is because the
condition (21) is also satisfied when both ghost and tachyonic instability are
present simultaneously. However, we recall that $F>0$ was actually a key
assumption to construct the compact phase space. Without this assumption, the
dynamical variables defined in (24) would not be bounded. Therefore, the
condition for the absence of ghost instability, namely $F>0$, is already
implicitly present in the compact phase space analysis. Thus, for a compact
phase space analysis, we can safely take the condition (36) to single out the
physically viable regions of phase space.

##  VI Applications

To illustrate the applicability of the model-independent compact dynamical
system formulation constructed in the previous section, we now apply the
formulation to three cosmological scenarios: two nonsingular bounce scenarios
and a cyclic cosmology. This is achieved by finding a cosmographic constraint
$j=j(q)$ corresponding to the specific cosmology. To do this, we proceed as
follows: Firstly, we specify a cosmological scenario as a given function of
time, $a(t)$ or $H(t)$. From there, it is straightforward to calculate the
deceleration parameter

$q(t)=-\frac{1}{aH^{2}}\frac{d^{2}a}{dt^{2}}=-1-\frac{\dot{H}}{H^{2}}\,,$ (37)

which we will invert to find $t=t(q)$. Next, we calculate the jerk parameter

$j(t)=\frac{1}{aH^{3}}\frac{d^{3}a}{dt^{3}}=-2-3q+\frac{\ddot{H}}{H^{3}}\,,$
(38)

and substitute $t=t(q)$ in that expression. As a result, we get $j=j(q)$. We
can now use the expression (34) for $q$ to close the dynamical system. In
practice, we note that $q$ and $j$ appear in the system only through the
combination $Q^{3}(j-q-2)$, so for each of the cosmological scenarios, we will
present the expression of this quantity explicitly.

The cosmic evolution corresponding to the fixed points can be obtained by
calculating the deceleration parameter $q=-1-\frac{\dot{H}}{H^{2}}$ using the
expression (34), except for those points for which the limiting value of the
deceleration parameter takes the form $q=1\pm\frac{0}{0}$. This occurs when

$Q=0\,\,,\,\,\bar{\Omega}+\bar{x}^{2}+\bar{K}=1\,.$ (39)

For such points, we cannot determine the cosmic evolution by this method. Let
us elaborate on this point from a more geometric point of view. Eq.(34) can be
written as

$\displaystyle\bar{\Omega}+\bar{x}^{2}+\bar{K}+(1-q)\bar{Q}^{2}=1\,,$ (40)

which shows that $q=constant$ corresponds in general to a family of
3-dimensional hypersurfaces within the 4-dimensional phase space. These
$q=constant$ hypersurfaces intersect at the curve given by Eq.(39). If a fixed
point lies on a specific $q=constant$ hypersurface, then its cosmology can be
determined unambiguously. However, if a fixed point lies at the intersection
of $q=constant$ hypersurfaces, i.e., on the curve (39), then one cannot
uniquely identify the cosmology corresponding to such a fixed point.222One can
think of this situation as being similar to what one encounters is plane polar
coordinates $(r,\theta)$. The origin does not have a unique $\theta$
coordinate as all $\theta=constant$ lines intersect at the origin. This is an
artefact of the particular choice of the compact variables (24).333One may not
encounter this issue if one uses, say, the Poincaré compactification
prescription [10]. The particular compactification prescription used here is
used because of its advantage in the clear depiction of bouncing trajectories.
As we will see below, all three examples we consider have fixed points like
this, so we need to be careful while interpreting the cosmology corresponding
to such points.

###  VI.1 Example-I

As our first example, let us take the ansatz

$H(t)=\frac{\beta t}{\left(\alpha^{2}+t^{2}\right)}\,,$ (41)

Where $-\infty<t<\infty$, $\alpha$ being the characteristic time-scale, i.e.,
$[\alpha]=[t]$ and $\beta>0$ is the dimensionless (stretching) constant. For
$|t|\gg\alpha$, $H(t)\approx\beta/t$ and for $|t|\ll\alpha$, $H(t)\approx\beta
t/\alpha^{2}$. This ansatz therefore represents a cosmic evolution that
connects a power law contraction phase $a(t)\sim(-t)^{\beta}$ to a power law
expanding phase $a(t)\sim t^{\beta}$ through a nonsingular bounce of the form
$a(t)\sim e^{\beta t^{2}/(2\alpha^{2})}$. This ansatz was considered in the
Ref.[20] where the authors were concerned about the behaviour of small shear
during such an evolution. Following [20], we consider the numerical value of
$\beta$ to be the same as $\alpha$ for the sake of simplicity.

The deceleration parameters can be calculated from Eq.(37)

$q=-1-\frac{\alpha^{2}-t^{2}}{\alpha t^{2}}\,,$ (42)

inverting, which we get

$t^{2}=\frac{\alpha^{2}}{1-\alpha(1+q)}\,.$ (43)

The jerk parameter can be calculated from Eq.(38). Using the above expression
for $t^{2}$ in terms of $q$ we get

$j=\frac{\alpha-2}{\alpha^{2}}[2-\alpha(2+3q)]\,.$ (44)

Finally, we obtain

$\bar{Q}^{3}(j-q-2)=\frac{2(2\alpha-3)}{\alpha}\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})-\frac{4\bar{Q}^{3}}{\alpha^{2}}(2\alpha^{2}-3\alpha+1)\,,$
(45)

which gives us the following dynamical system

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+2\bar{K}+\bar{x}\bar{Q}-2\bar{x}^{2}\bar{K}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)-\bar{\Omega}(1+3\omega)+3\bar{x}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})$
(46a)
$\displaystyle-\frac{2(2\alpha-3)}{\alpha}\bar{x}\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})+\frac{4\bar{x}\bar{Q}^{3}}{\alpha^{2}}\left(2\alpha^{2}-3\alpha+1\right)+2\bar{x}\bar{K}\bar{Q}\Big{]}\,,$
$\displaystyle\frac{d\bar{\Omega}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{\Omega}}{3}\Big{[}3(\bar{x}\bar{\Omega}+\bar{Q}\bar{\Omega}-\bar{Q})(1+\omega)-\frac{2(2\alpha-3)}{\alpha}\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})+\frac{4\bar{Q}^{3}}{\alpha^{2}}\left(2\alpha^{2}-3\alpha+1\right)+2\bar{K}\bar{Q}-4\bar{x}\bar{Q}^{2}$
(46b) $\displaystyle-4\bar{x}\bar{K}-2\bar{x}^{2}\bar{Q}+2\bar{x}^{3}\Big{]},$
$\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}3\bar{Q}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-2\bar{K}-4\bar{Q}^{2}+2\bar{x}\bar{Q}-\frac{2(2\alpha-3)}{\alpha}\bar{Q}^{2}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})+\frac{4\bar{Q}^{4}}{\alpha^{2}}\left(2\alpha^{2}-3\alpha+1\right)$
(46c)
$\displaystyle+2\bar{K}\bar{Q}^{2}+2\left(1-\bar{\Omega}-\bar{x}^{2}\right)-4\bar{x}\bar{Q}\bar{K}-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}-4\bar{x}\bar{Q}^{3}\Big{]}\,,$
$\displaystyle\frac{d\bar{K}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{K}}{6}\Big{[}4\bar{x}-4\bar{Q}-8\bar{x}\bar{K}-8\bar{x}\bar{Q}^{2}-4\bar{x}^{2}\bar{Q}+4\bar{x}^{3}+6\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-\frac{4(2\alpha-3)}{\alpha}\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})$
(46d)
$\displaystyle+\frac{8\bar{Q}^{3}}{\alpha^{2}}\left(2\alpha^{2}-3\alpha+1\right)+4\bar{K}\bar{Q}\Big{]}\,.$

The fixed points of this dynamical system for the case of dust ($\omega=0$)
and assuming the parameter value $\alpha=2/3$ are listed in Table 1.

| Fixed
---
points
Coordinates $(\bar{x},\bar{\Omega},\bar{Q},\bar{K})$ | Jacobian eigenvalues | Stability | Cosmology
$B_{1\pm}$ | $(\pm 1,0,0,0)$ | $\big{\\{}\pm\frac{4}{3},\pm\frac{4}{3},\pm\frac{2}{3},\pm\frac{2}{3}\big{\\}}$ |  |  $(-)$ attractor
---
$(+)$ repeller
|  $|t|\rightarrow\infty$ limit of
---
$a(t)\sim|t|^{\frac{2}{3}}$
$B_{2\pm}$ | $(\pm 0.53,0,\mp 1.20,0)$ | $\big{\\{}\mp 1.70,\pm 1.20,\mp 0.40,\mp 0.35\big{\\}}$ | Saddle | $a(t)\sim|t|^{\frac{2}{3}}$
$B_{3\pm}$ | $(\pm 0.923,0,\pm 0.545,0)$ | $\big{\\{}\mp 0.775,\mp 0.615,\mp 0.544,\pm 0.181\big{\\}}$ | Saddle | $a(t)\sim|t|^{\frac{2}{3}}$
$B_{4\pm}$ | $\left(0,\frac{2}{3},\pm\sqrt{\frac{2}{3}},0\right)$ | $\left\\{\pm\frac{\sqrt{73}+5}{6\sqrt{6}},\mp\sqrt{\frac{2}{3}},\pm\frac{1}{3}\sqrt{\frac{2}{3}},\mp\frac{\sqrt{73}-5}{6\sqrt{6}}\right\\}$ | Saddle | $a(t)\sim|t|^{\frac{2}{3}}$
$B_{5\pm}$ | $\left(\pm\frac{3}{\sqrt{17}},\frac{12}{17},\pm\frac{2}{\sqrt{17}},0\right)$ | $\left\\{\mp\frac{2(\sqrt{17}+\sqrt{170})}{51},\pm\frac{8}{3\sqrt{17}},\pm\frac{2}{\sqrt{17}},\pm\frac{2(\sqrt{170}-\sqrt{17})}{51}\right\\}$ | Saddle | $a(t)\sim|t|^{\frac{2}{3}}$
$B_{6\pm}$ | $\left(\pm\frac{1}{\sqrt{3}},0,0,\frac{2}{3}\right)$ | $\left\\{\mp\frac{4}{3\sqrt{3}},\mp\frac{2}{3\sqrt{3}},0,0\right\\}$ |  |  $(+)$ attractor
---
or saddle
$(-)$ repeller
or saddle
Indeterminate
$B_{7}$ | $\left(0,\frac{2}{5},0,\frac{3}{5}\right)$ | $\left\\{-\sqrt{\frac{2}{15}},\sqrt{\frac{2}{15}},0,0\right\\}$ | Saddle | Indeterminate
Table 1: Fixed points, Jacobian eigenvalues, stability and the respective
cosmologies for the cosmic evolution corresponding to the ansatz (41) with
$\omega=0$ and $\alpha=2/3$. Cosmic evolution corresponding to the fixed
points $B_{1\pm},\,B_{6\pm},\,B_{7}$ cannot be determined as, according to the
expression (34), $q$ is indeterminate.

In general, the phase space is 4-dimensional. Figure 1 shows different
projections of a phase trajectory that corresponds to a non-singular bounce.
The last panel in the figure proves that throughout its entire evolution
history, the bouncing cosmology encounters neither ghost nor tachyonic
instabilities.

One can note that the cosmology corresponding to the fixed points
$B_{1\pm},\,B_{6\pm},\,B_{7}$ cannot be determined from the deceleration
parameter, but has to be found in other ways. We note that, for all these
fixed points, $Q=0$, i.e. $H\equiv\frac{\dot{a}}{a}=0$. This can imply either
a static solution $a(t)=constant$, or the $|t|\rightarrow\infty$ limit of the
asymptotic form of the evolution ($a(t)\sim|t|^{2/3}$) where $\dot{a}(t)$
diverges faster than $a(t)$, essentially making $H(t)\rightarrow 0$. We
identify the cosmology corresponding to the fixed points $B_{1\pm}$ with the
latter option, which is supported by the first panel of Figure 1. This gives a
picture of a bouncing cosmology consistent with the ansatz (41).
Unfortunately, for $B_{6\pm},\,B_{7}$, we still cannot exactly identify the
corresponding cosmology.

Figure 1: Parametric plots $\bar{x}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$
(upper left), $\bar{\Omega}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$ (upper
right) and $\bar{K}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$ (lower left) for a
trajectory passing through the point
$(\bar{x},\bar{\Omega},\bar{Q},\bar{K})=(0,0.14,0,0.15)$ for the ansatz (41)
with $\omega=0,\,\alpha=2/3$. The lower right panel shows the quantity
$\frac{2\bar{x}}{\bar{Q}^{3}(j-q-2)-2\bar{K}\bar{Q}}$ remains positive
throughout the entire evolution history of the bouncing cosmology, signalling
a stable bouncing solution.

For the special case when the bounce occurs in a spatially flat and vacuum
scenario, $\bar{\Omega}=0$ and $\bar{K}=0$, the dynamical system reduces to

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+\bar{x}\bar{Q}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)-\frac{2(2\alpha-3)}{\alpha}\bar{x}\bar{Q}(1-\bar{x}^{2})+\frac{4\bar{x}\bar{Q}^{3}}{\alpha^{2}}\left(2\alpha^{2}-3\alpha+1\right)\Big{]}\,,$
(47a) $\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\bar{x}\bar{Q}-4\bar{Q}^{2}+\Big{(}2-\frac{2(2\alpha-3)\bar{Q}^{2}}{\alpha}\Big{)}(1-\bar{x}^{2})+\frac{4\left(2\alpha^{2}-3\alpha+1\right)\bar{Q}^{4}}{\alpha^{2}}-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}-4\bar{x}\bar{Q}^{3}\Big{]}\,.$
(47b)

In this case, the phase space is 2-dimensional. The fixed points are presented
in Table 2.

Fixed points | Coordinates $(\bar{x},\bar{Q})$ | Jacobian eigenvalues | Stability | Cosmology
---|---|---|---|---
$A_{1\pm}$ | $(\pm 1,0)$ | $\left\\{\pm\frac{4}{3},\pm\frac{2}{3}\right\\}$ |  |  $(-)$ attractor
---
$(+)$ repeller
Indeterminate
$A_{2\pm}$ | $(\pm 0.53,\mp 1.20)$ | $\big{\\{}\mp 1.71,\pm 1.12\big{\\}}$ | Saddle | $a(t)\sim|t|^{\frac{2}{3}}$
$A_{3\pm}$ | $(\pm 0.923,\pm 0.545)$ | $\\{\mp 0.78,\mp 0.54\\}$ |  |  $(+)$ attractor
---
$(-)$ repeller
$a(t)\sim|t|^{\frac{2}{3}}$
Table 2: Fixed points, Jacobian eigenvalues, stability and the respective
cosmologies corresponding to the ansatz (41) with $\alpha=\frac{2}{3}$ for the
spatially flat vacuum case.

The phase portrait is presented in Figure 2, where we also highlight the
physically viable region according to the condition (21). One can note that,
in general, all the non-singular bouncing trajectories do go out of the
physically viable region for some time during their course of evolution,
except for the single trajectory passing through the origin. All other
bouncing trajectories encounter either ghost or tachyonic instability, at
least in the vicinity of the bounce. Therefore, one can conclude that even
though globally stable non-singular bouncing solutions of the form (41) are in
principle possible in $f(R)$ gravity, they are very fine-tuned.

Figure 2: Phase portrait for the spatially flat vacuum case for the ansatz
(41) with $\alpha=2/3$. The shaded region is the physically viable region
where the $f(R)$ theory is free from ghost and tachyonic instability
($F>0,\,F^{\prime}>0$). The green ellipse corresponds to the invariant
submanifold $q=\frac{1-\alpha}{\alpha}=\frac{1}{2}$.

The time evolution of the dynamical variables for the trajectory passing
through the origin is shown in Figure 3.

Figure 3:  The left and right panels depict the time evolution of the
variables $\bar{x}$ and $\bar{Q}$ with respect to $\bar{\tau}$ for the
trajectory passing through the origin. $\alpha=2/3$ is considered.

The $f(R)$ phase space corresponding to the ansatz (41) has two invariant
submanifolds corresponding to
$q=\frac{1-\alpha}{\alpha},\,\frac{2-\alpha}{\alpha}$. This can be seen
clearly by inserting the expression of $j=j(q)$ from Eq.(44) into the first
equation of (15)

$\frac{dq}{d\tau}=2\left(q-\frac{1-\alpha}{\alpha}\right)\left(q-\frac{2-\alpha}{\alpha}\right)\,.$
(48)

They represent 3-dimensional surfaces in the 4-dimensional phase space spanned
by the variables $\\{\bar{x},\bar{\Omega},\bar{Q},\bar{K}\\}$, whose equation
can be found using Eq.(34). The intersection of the submanifold
$q=\frac{1-\alpha}{\alpha}$ with the spatially flat vacuum invariant
submanifold for $\alpha=\frac{2}{3}$ is shown in Figure 2 as the green curve.
As can be seen from the figure, the bouncing trajectories corresponding to
Eq.(41) are contained within the region bounded by this submanifold. Moreover,
the submanifold is asymptotic to all the bouncing trajectories, which is
consistent with the fact that the past and future asymptotic state of the
cosmology (41) have a power law evolution $a(t)\sim|t|^{\alpha}$ (recall that
we have taken the same numerical value for $\beta$ and $\alpha$ for
simplicity). The trajectories that lie outside the region bounded by the
invariant submanifold $q=\frac{1-\alpha}{\alpha}$ are totally different
cosmological solutions that still satisfy the same cosmographic relation (44).
The other invariant submanifold $q=\frac{2-\alpha}{\alpha}$ lies outside the
compact region and is not of any physical interest.

###  VI.2 Example-II

The second example that we consider has a scale factor that evolves according
to

$a(t)=a_{0}\cosh{\lambda t}\,,$ (49)

where again $-\infty<t<\infty$ and $\lambda$ is a positive constant with
$[\lambda]=[t]^{-1}$. In asymptotic past and asymptotic future, the behaviour
of the scale factor is $e^{-\lambda t}$ and $e^{\lambda t}$ respectively.
Therefore, this ansatz represents a cosmic evolution that connects a
contracting De-Sitter (deflationary) phase to an expanding de-Sitter
(inflationary) phase through a non-singular bounce. This ansatz was introduced
in Ref.[21] as a model of inflation preceded by a bounce. The authors in
Ref.[21] try to reconstruct the $f(R)$, which comes out in a complicated form.
Our approach presents a way around the cumbersome exercise of reconstruction.

Hubble and deceleration parameters for the above-mentioned scale factor are

$\displaystyle H$ $\displaystyle=$ $\displaystyle\lambda\tanh{{\lambda t}}\,,$
(50a) $\displaystyle q$ $\displaystyle=$
$\displaystyle-1-\operatorname{csch}^{2}{\lambda t}\,.$ (50b)

Eq.(50a) can be inverted to get

$\operatorname{csch}^{2}{\lambda t}=-1-q\,.$ (51)

The jerk parameter can be calculated from Eq.(38). Using the above expression
for $\operatorname{csch}^{2}{\lambda t}$, we obtain

$j=-q\,.$ (52)

Finally, we obtain

$\bar{Q}^{3}(j-q-2)=-2\bar{Q}^{3}(q+1)=2\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})-4\bar{Q}^{3}\,,$
(53)

which gives us the dynamical system

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+2\bar{K}+\bar{x}\bar{Q}-2\bar{x}^{2}\bar{K}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)-\bar{\Omega}(1+3\omega)+3\bar{x}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})$
(54a)
$\displaystyle-2\bar{x}\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})+4\bar{x}\bar{Q}^{3}+2\bar{x}\bar{K}\bar{Q}\Big{]}\,,$
$\displaystyle\frac{d\bar{\Omega}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{\Omega}}{3}\Big{[}3(\bar{x}\bar{\Omega}+\bar{Q}\bar{\Omega}-\bar{Q})(1+\omega)-2\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})+4\bar{Q}^{3}+2\bar{K}\bar{Q}-4\bar{x}\bar{Q}^{2}-4\bar{x}\bar{K}-2\bar{x}^{2}\bar{Q}+2\bar{x}^{3}\Big{]}\,,$
(54b) $\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}3\bar{Q}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-2\bar{K}-4\bar{Q}^{2}+2\bar{x}\bar{Q}-2\bar{Q}^{2}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})+4\bar{Q}^{4}+2\bar{K}\bar{Q}^{2}+2\left(1-\bar{\Omega}-\bar{x}^{2}\right)$
(54c)
$\displaystyle-4\bar{x}\bar{Q}\bar{K}-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}-4\bar{x}\bar{Q}^{3}\Big{]}\,,$
$\displaystyle\frac{d\bar{K}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{K}}{6}\Big{[}4\bar{x}-4\bar{Q}-8\bar{x}\bar{K}-8\bar{x}\bar{Q}^{2}-4\bar{x}^{2}\bar{Q}+4\bar{x}^{3}+6\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-4\bar{Q}(1-\bar{\Omega}-\bar{x}^{2}-\bar{K})+8\bar{Q}^{3}+4\bar{K}\bar{Q}\Big{]}\,.$

The fixed points of this dynamical system are listed in Table 3 and their
stability and respective cosmologies are listed in Table 4.

| Fixed
---
points
Coordinates $(\bar{x},\bar{\Omega},\bar{Q},\bar{K})$ | Existence | Jacobian eigenvalues
$B_{1\pm}$ | $(\pm 1,0,0,0)$ | Always | $\left\\{\pm\frac{4}{3},\pm\frac{4}{3},\pm\frac{2}{3},\pm\frac{2}{3}\right\\}$
$B_{2\pm}$ | $\left(0,0,\pm\frac{1}{\sqrt{2}},0\right)$ | Always | $\left\\{\mp\frac{1}{\sqrt{2}},\mp\frac{\sqrt{2}}{3},\mp\frac{\sqrt{2}}{3},\pm\frac{1}{3\sqrt{2}}\right\\}$
$B_{3\pm}$ | $\left(\pm\frac{1}{3},0,\pm\frac{2}{3},0\right)$ | Always | $\left\\{\mp\frac{8}{9},\mp\frac{4}{9},\mp\frac{4}{9},\mp\frac{2}{9}\right\\}$
$B_{4+}$ | $\left(\frac{\bar{Q}+\sqrt{3+\bar{Q}^{2}}}{3},0,\bar{Q},\frac{6-11\bar{Q}^{2}-2\bar{Q}\sqrt{3+\bar{Q}^{2}}}{9}\right)$ | $-0.939\leq\bar{Q}\leq 0.59$ | $\left\\{0,\frac{2\bar{Q}}{3},\frac{2\left(\bar{Q}-2\sqrt{\bar{Q}^{2}+3}\right)}{9},\frac{\left(-2\sqrt{\bar{Q}^{2}+3}-5\bar{Q}\right)}{9}\right\\}$
$B_{4-}$ | $\left(\frac{\bar{Q}-\sqrt{3+\bar{Q}^{2}}}{3},0,\bar{Q},\frac{6-11\bar{Q}^{2}+2\bar{Q}\sqrt{3+\bar{Q}^{2}}}{9}\right)$ | $-0.59\leq\bar{Q}\leq 0.939$ | $\left\\{0,\frac{2\bar{Q}}{3},\frac{2\left(\bar{Q}+2\sqrt{\bar{Q}^{2}+3}\right)}{9},\frac{\left(2\sqrt{\bar{Q}^{2}+3}-5\bar{Q}\right)}{9}\right\\}$
$B_{5}$ | $\left(-\frac{\bar{Q}}{2},\frac{4-7\bar{Q}^{2}}{10},\bar{Q},\frac{12-11\bar{Q}^{2}}{20}\right)$ | $-\frac{2}{\sqrt{7}}\leq\bar{Q}\leq\frac{2}{\sqrt{7}}$ | $\left\\{0,\frac{2\bar{Q}}{3},\frac{2\bar{Q}}{3}-\frac{1}{3}\sqrt{\frac{19\bar{Q}^{2}}{10}+\frac{6}{5}},\frac{\left(\sqrt{190\bar{Q}^{2}+120}+20\bar{Q}\right)}{30}\right\\}$
Table 3: Fixed points, their existence conditions and the Jacobian eigenvalues for the cosmic evolution corresponding to the ansatz (49) with $\omega=0$. Fixed points | Coordinates $(\bar{x},\bar{\Omega},\bar{Q},\bar{K})$ | Stability | Cosmology
---|---|---|---
$B_{1\pm}$ | $(\pm 1,0,0,0)$ |  |  $(-)$ attractor
---
$(+)$ repeller
Indeterminate
$B_{2\pm}$ | $\left(0,0,\pm\frac{1}{\sqrt{2}},0\right)$ | Saddle |  | $a(t)=a_{0}e^{ct}$
---
(De-sitter)
$B_{3\pm}$ | $\left(\pm\frac{1}{3},0,\pm\frac{2}{3},0\right)$ |  |  $(+)$ attractor
---
$(-)$ repeller
| $a(t)=a_{0}e^{ct}$
---
(De-sitter)
$B_{4+}$ | $\left(\frac{\bar{Q}\pm\sqrt{3+\bar{Q}^{2}}}{3},0,\bar{Q},\frac{6-11\bar{Q}^{2}\mp 2\bar{Q}\sqrt{3+\bar{Q}^{2}}}{9}\right)$ |  | Stable for $-\frac{2}{\sqrt{7}}<\bar{Q}<0$
---
saddle otherwise
$a(t)=a_{0}+Ct$
$B_{4-}$ | $\left(\frac{\bar{Q}\pm\sqrt{3+\bar{Q}^{2}}}{3},0,\bar{Q},\frac{6-11\bar{Q}^{2}\mp 2\bar{Q}\sqrt{3+\bar{Q}^{2}}}{9}\right)$ |  | Unstable for $0<\bar{Q}<\frac{2}{\sqrt{7}}$
---
saddle otherwise
$a(t)=a_{0}+Ct$
$B_{5}$ | $\left(-\frac{\bar{Q}}{2},\frac{4-7\bar{Q}^{2}}{10},\bar{Q},\frac{12-11\bar{Q}^{2}}{20}\right)$ | Saddle | $a(t)=a_{0}+Ct$
Table 4: Fixed points, stability and the respective cosmologies for the cosmic
evolution corresponding to the ansatz (49) with $\omega=0$. For those fixed
points which have $q=1\pm\frac{0}{0}$ form, cosmology is indeterminate.

For the points $B_{1\pm}$, we cannot use the same argument that we used for
model-I (41), as for the model (49), both $a(t)$ and $\dot{a}(t)$ diverges at
the same rate as $|t|\rightarrow\infty$. So, unfortunately, we are unable to
determine their cosmology unambiguously. We recognize this as a limitation of
our method.

Figures 4 show different projections of two-phase trajectories that correspond
to a non-singular bounce. The last panel in the figure proves that throughout
its entire evolution history, bouncing cosmology encounters neither ghost nor
tachyonic instability.

Figure 4: Parametric plots $\bar{x}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$
(upper left), $\bar{\Omega}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$ (upper
right) and $\bar{K}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$ (lower left) for a
trajectory passing through the point
$(\bar{x},\bar{\Omega},\bar{Q},\bar{K})=(0,0.14,0,0.45)$ for the ansatz (49)
with $\omega=0$. The lower right panel shows the quantity
$\frac{2\bar{x}}{\bar{Q}^{3}(j-q-2)-2\bar{K}\bar{Q}}$ remains positive
throughout the entire evolution history of the bouncing cosmology, signaling a
stable bouncing solution.

For the spatially flat vacuum case, the dynamical system reduces to

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+\bar{x}\bar{Q}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)-2\bar{x}\bar{Q}(1-\bar{x}^{2})+4\bar{x}\bar{Q}^{3}\Big{]}\,,$
(55a) $\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}-4\bar{Q}^{2}+2\bar{x}\bar{Q}+\left(2-2\bar{Q}^{2}\right)(1-\bar{x}^{2})+4\bar{Q}^{4}-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}-4\bar{x}\bar{Q}^{3}\Big{]}\,.$
(55b)

whose fixed points are listed in Table 5.

Fixed points | Coordinates $(\bar{x},\bar{Q})$ | Jacobian eigenvalues | Stability | Cosmology
---|---|---|---|---
$A_{1\pm}$ | $(\pm 1,0)$ | $\left\\{\pm\frac{4}{3},\pm\frac{2}{3}\right\\}$ |  |  $(-)$ attractor
---
$(+)$ repeller
Indeterminate
$A_{2\pm}$ | $\left(0,\pm\frac{1}{\sqrt{2}}\right)$ | $\left\\{\mp\frac{\sqrt{2}}{3},\pm\frac{1}{3\sqrt{2}}\right\\}$ | Saddle |  | $a(t)=a_{0}e^{ct}$
---
(De-sitter)
$A_{3\pm}$ | $\left(\pm\frac{1}{3},\pm\frac{2}{3}\right)$ | $\left\\{\mp\frac{4}{9},\mp\frac{2}{9}\right\\}$ |  |  $(+)$ attractor
---
$(-)$ repeller
| $a(t)=a_{0}e^{ct}$
---
(De-sitter)
$A_{4\pm}$ | $(\approx\pm 0.81,\approx\pm 0.59)$ | $\\{\mp 0.68,\pm 0.39\\}$ | Saddle | $a(t)=a_{0}+Ct$
$A_{5\pm}$ | $(\approx\pm 0.34,\approx\mp 0.94)$ | $\\{\mp 1.08,\mp 0.625\\}$ |  |  $(+)$ attractor
---
$(-)$ repeller
$a(t)=a_{0}+Ct$
Table 5: Fixed points, Jacobian eigenvalues, stabilities and the respective
cosmologies for the ansatz (49) for the spatially flat vacuum case.

The phase portrait is presented in Figure 5, where the physically viable
region is highlighted. Amazingly, one can notice that there are no bouncing
trajectories that remain completely within the physically viable region for
the entire course of evolution. If one considers a spatially flat evolution of
the form (49) in a vacuum scenario, within the framework of $f(R)$ gravity, it
is bound to encounter either ghost or tachyonic instability at least in the
vicinity of the bounce. This makes the model considered in Ref.[21] not
physically viable.

Figure 5: Phase portrait for the spatially flat vacuum case for the ansatz
(49). The shaded region is the physically viable region where the $f(R)$
theory is free from ghost and tachyonic instability ($F>0,\,F^{\prime}>0$).
The green ellipse corresponds to the invariant submanifold $q=-1$. The brown
ellipse corresponds to the invariant submanifold $q=0$.

The time evolution of the dynamical variables for the trajectory passing
through the origin is shown in Figure 6.

Figure 6: The left and right panels, respectively, depict the time evolution
of the variables $\bar{x}$ and $\bar{Q}$ with respect to $\bar{\tau}$ for the
trajectory passing through the origin.

The $f(R)$ phase space corresponding to the ansatz (49) has two invariant
submanifolds corresponding to $q=0,\,-1$. This can be seen clearly by
inserting the expression of $j=j(q)$ from Eq.(44) into the first equation of
(15)

$\frac{dq}{d\tau}=2q(q+1)\,.$ (56)

The equations of the corresponding 3-dimensional surfaces can be found using
Eq.(34). The intersections of these invariant submanifolds with the spatially
flat vacuum invariant submanifold are shown in Figure 5 as the green and brown
curves. As can be seen from the figure, the bouncing trajectories
corresponding to Eq.(49) are contained within the region bounded by the
submanifold $q=-1$. This submanifold is asymptotic to all the bouncing
trajectories, which is consistent with the fact that the past and future
asymptotic states of the cosmology (49) are De-Sitter. The trajectories that
lie outside the region bounded by the invariant submanifold $q=-1$ are
different cosmological solutions that still satisfy the same cosmographic
relation (52).

###  VI.3 Example III

The third and the last example that we consider is a non-singular cyclic
cosmology

$a(t)=\frac{1}{2\Lambda}\left[2-\cos(2\zeta t)\right]\,,$ (57)

where $\Lambda$ is a dimensionless constant and $\zeta$ is a constant with
dimension $[\zeta]=[t]^{-1}$. It has been argued that a cyclic cosmological
model is compatible with current observations if the cosmological constant
happens to be decaying [22]. The simple example we use here is taken from
Ref.[23], where the authors investigated whether such a cyclic cosmology can
be supported by a third-order polynomial form of $f(R)$. The authors in
Ref.[23] did not investigate the stability of such a model. As we will see
below, a simple cyclic ansatz like (57) is totally unstable in $f(R)$ gravity,
irrespective of its form.

The deceleration and the jerk parameter is calculated from Eq.(37) and Eq.(38)
respectively

$\displaystyle q$ $\displaystyle=$ $\displaystyle\frac{\cos^{2}(2\zeta
t)-2\cos(2\zeta t)}{\sin^{2}(2\zeta t)}\,,$ (58a) $\displaystyle j$
$\displaystyle=$ $\displaystyle-[\cot(2\zeta t)-2\csc(2\zeta
t)]^{2}=\frac{\cos^{2}(2\zeta t)-4}{1-\cos^{2}(2\zeta t)}-2q\,.$ (58b)

From the expression of the deceleration parameter above one can solve

$\cos(2\zeta t)=\frac{1-\sqrt{1+q+q^{2}}}{1+q}\,,$ (59)

using which one arrives at the expression

$j=-2\left(\sqrt{q^{2}+q+1}+1\right)-q\,.$ (60)

Finally, one obtains

$\bar{Q}^{3}(j-q-2)=-2\bar{Q}\left(\sqrt{{\left(\bar{K}+\bar{x}^{2}+\bar{\Omega}-1\right)\left(\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)}+3\bar{Q}^{4}}+\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)\,.$
(61)

We are now in a position to write the dynamical system

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+2\bar{K}+\bar{x}\bar{Q}-2\bar{x}^{2}\bar{K}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)-\bar{\Omega}(1+3\omega)+3\bar{x}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})$
(62a)
$\displaystyle+2\bar{x}\bar{K}\bar{Q}+2\bar{x}\bar{Q}\left(\sqrt{{\left(\bar{K}+\bar{x}^{2}+\bar{\Omega}-1\right)\left(\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)}+3\bar{Q}^{4}}+\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)\Big{]},$
$\displaystyle\frac{d\bar{\Omega}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{\Omega}}{3}\Big{[}3(\bar{x}\bar{\Omega}+\bar{Q}\bar{\Omega}-\bar{Q})(1+\omega)+2\bar{K}\bar{Q}-4\bar{x}\bar{Q}^{2}-4\bar{x}\bar{K}-2\bar{x}^{2}\bar{Q}+2\bar{x}^{3}$
(62b)
$\displaystyle+2\bar{Q}\left(\sqrt{{\left(\bar{K}+\bar{x}^{2}+\bar{\Omega}-1\right)\left(\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)}+3\bar{Q}^{4}}+\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)\Big{]},$
$\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}3\bar{Q}\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})-2\bar{K}-4\bar{Q}^{2}+2\bar{x}\bar{Q}+2\bar{K}\bar{Q}^{2}+2\left(1-\bar{\Omega}-\bar{x}^{2}\right)-4\bar{x}\bar{Q}\bar{K}-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}$
(62c)
$\displaystyle-4\bar{x}\bar{Q}^{3}+2\bar{Q}^{2}\left(\sqrt{{\left(\bar{K}+\bar{x}^{2}+\bar{\Omega}-1\right)\left(\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)}+3\bar{Q}^{4}}+\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)\Big{]},$
$\displaystyle\frac{d\bar{K}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{\bar{K}}{6}\Big{[}4\bar{x}-4\bar{Q}-8\bar{x}\bar{K}-8\bar{x}\bar{Q}^{2}-4\bar{x}^{2}\bar{Q}+4\bar{x}^{3}+6\bar{\Omega}(1+\omega)(\bar{Q}+\bar{x})+4\bar{K}\bar{Q}$
(62d)
$\displaystyle+4\bar{Q}\left(\sqrt{{\left(\bar{K}+\bar{x}^{2}+\bar{\Omega}-1\right)\left(\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)}+3\bar{Q}^{4}}+\bar{K}+3\bar{Q}^{2}+\bar{x}^{2}+\bar{\Omega}-1\right)\Big{]}.$

An important point is worth mentioning here. The method we devised in this
paper depends crucially on the deceleration and the jerk parameter, both of
which are periodic for this model. One has to be careful while dealing with
periodic functions. The fixed point structure that will be presented below
strictly holds within one period
$t\in\left[-\frac{\pi}{\zeta},\frac{\pi}{\zeta}\right]$.

The fixed points are listed in Table 6.

| Fixed
---
points
Coordinates $(\bar{x},\bar{\Omega},\bar{Q},\bar{K})$ | Jacobian eigenvalues | Stability | Cosmology
$B_{1+}$ | $(1,0,0,0)$ | $\left\\{\frac{4}{3},\frac{4}{3},\frac{2}{3},\frac{2}{3}\right\\}$ |  | always
---
repeller
Indeterminate
$B_{1-}$ | $(-1,0,0,0)$ | $\left\\{-\frac{4}{3},-\frac{4}{3},-\frac{2}{3},-\frac{2}{3}\right\\}$ |  | always
---
attractor
Indeterminate
$B_{2+}$ | $\left(\frac{1}{\sqrt{3}},0,0,\frac{2}{3}\right)$ | $\\{-0.561285,\,0.0458372,\,0.116841\,\pm 0.4318i\\}$ | Saddle | Indeterminate
$B_{2-}$ | $\left(-\frac{1}{\sqrt{3}},0,0,\frac{2}{3}\right)$ | $\\{-0.296296,\,-0.133618,\,0.458974\,\pm 0.50229i\\}$ | Saddle | Indeterminate
$B_{3}$ | $\left(0,\frac{2}{5},0,\frac{3}{5}\right)$ | $\\{0.66674,0.431396,0.158996,-0.017242\\}$ | Saddle | Indeterminate
Table 6: Fixed points, Jacobian eigenvalues and stabilities corresponding to
the ansatz (57) with $\omega=0$. Interestingly, all the fixed points have a
deceleration parameter of the form $q=1\pm\frac{0}{0}$. Therefore, their
cosmologies are indeterminate.

We note that, interestingly, all the fixed points have the deceleration
parameter of the form $q=1\pm\frac{0}{0}$. We cannot use the argument that we
used for model-I (41), as the cosmology is cyclic. So, unfortunately, we are
unable to determine their cosmology unambiguously. We recognize this as a
limitation of our method.

Parametric plots, however, clearly show the cyclic nature of the cosmology. In
Figure 7, we show the parametric plots $\bar{x}(\bar{\tau})$ v/s
$\bar{Q}(\bar{\tau})$, $\bar{\Omega}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$
and $\bar{K}(\bar{\tau})$ v/s $\bar{Q}(\bar{\tau})$ for a trajectory passing
through the point $(\bar{x},\bar{\Omega},\bar{Q},\bar{K})=(0,0.08,-0.4,0.68)$.
The last panel of Figure 7 investigates the physical viability of the cyclic
cosmology (57) in $f(R)$ gravity. What the plot reveals is that the expansion
phases of the cyclic cosmology are always unstable. This makes a simplistic
cyclic model like (57) totally non-viable in $f(R)$ gravity.

Figure 7: Parametric plots for the ansatz (57) of $\bar{x}(\bar{\tau})$ v/s
$\bar{Q}(\bar{\tau})$ (upper left), $\bar{\Omega}(\bar{\tau})$ v/s
$\bar{Q}(\bar{\tau})$ (upper right) and $\bar{K}(\bar{\tau})$ v/s
$\bar{Q}(\bar{\tau})$ (lower left) for a trajectory passing through the point
$(\bar{x},\bar{\Omega},\bar{Q},\bar{K})=(0,0.08,-0.4,0.68)$. The lower right
panel shows the time evolution of $\bar{Q}$ and the quantity
$\frac{2\bar{x}}{\bar{Q}^{3}(j-q-2)-2\bar{K}\bar{Q}}$ with respect to time
variable $\bar{\tau}$.

Considering the spatially flat vacuum case, one can write the dynamical system
as

$\displaystyle\frac{d\bar{x}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}2\left(-1+\bar{x}\bar{Q}+2\bar{Q}^{2}-2\bar{x}^{2}\bar{Q}^{2}-\bar{x}^{3}\bar{Q}+\bar{x}^{4}\right)+2\bar{x}\bar{Q}(3\bar{Q}^{2}+\bar{x}^{2}-1)$
(63a)
$\displaystyle+2\bar{x}\bar{Q}\sqrt{\left(\bar{x}^{2}-1\right)\left(3\bar{Q}^{2}+\bar{x}^{2}-1\right)+3\bar{Q}^{4}}\Big{]},$
$\displaystyle\frac{d\bar{Q}}{d\bar{\tau}}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\Big{[}-4\bar{Q}^{2}+2\bar{x}\bar{Q}+2\left(1-\bar{x}^{2}\right)-2\bar{x}^{2}\bar{Q}^{2}+2\bar{x}^{3}\bar{Q}-4\bar{x}\bar{Q}^{3}+2\bar{Q}^{2}(3\bar{Q}^{2}+\bar{x}^{2}-1)$
(63b)
$\displaystyle+2\bar{Q}^{2}\sqrt{\left(\bar{x}^{2}-1\right)\left(3\bar{Q}^{2}+\bar{x}^{2}-1\right)+3\bar{Q}^{4}}\Big{]}\,,$

whose fixed points are listed in Table 7.

| Fixed
---
points
| Coordinates
---
$(\bar{x},\bar{Q})$
Stability | Eigenvalues | Cosmology
$A_{1\pm}$ | $(\pm 1,0)$ |  |  $(-)$ attractor
---
$(+)$ repeller
$\left\\{\pm\frac{4}{3},\pm\frac{2}{3}\right\\}$ | Indeterminate
Table 7: Stability and cosmology of critical points for example III.

The corresponding phase portrait is shown in Figure 8.

Figure 8: Phase portrait for the spatially flat vacuum case for the ansatz
(57). The shaded region represents absence of ghost and tachyonic instability.

The time evolution of the dynamical variables are shown in Figure 9.

Figure 9: The left and right panels depict the temporal evolution of the
variables $\bar{x}$ and $\bar{Q}$ for the variable $\bar{\tau}$.

We note that the cyclic nature of the evolution is not apparent from these
plots, as these plots are generated using the numerical solutions of the
dynamical system, which gives the fixed point structure strictly within one
period $t\in\left[-\frac{\pi}{\zeta},\frac{\pi}{\zeta}\right]$.

##  VII Conclusion

In this paper, a new model independent dynamical systems approach based on
cosmographic parameters [12] was developed further by compactifying the phase
space for a general $f(R)$ theory of gravity and giving the corresponding
physical viability conditions on the phase space.

In order to showcase the applicability of this approach, three examples
displaying cyclic and/or bounce behaviour were studied. Fixed points as well
as the associated stability for all systems were found. For the examples
involving bouncing cosmology, it was found that the past and future asymptotic
states of the bounce manifest themselves as invariant submanifolds in the
phase space. The phase trajectories corresponding to bounces are confined
within the region bounded by these submanifolds. An interesting observation
made was that cyclic and bounce behaviour was only seen clearly when studying
the parametric plots with positive spatial curvature initial conditions. This
is consistent with a recent paper [18] in which it was found through the use
of the explicit functional form of the $f(R)$ as well as the constraints on
the compact phase space variables that the Hu-Sawicki model with $n=1=C_{1}$
displays no physical viable bounces in a spatially flat universe. The question
remains whether this is a general statement, i.e., whether $f(R)$ models
associated to a cyclic or bounce scale factor in FLRW cosmology can only
display physically viable cyclic or bounce behaviour when the curvature is
positive.

All three examples showcased fixed points in which $q=1\pm\frac{0}{0}$ (34)
and thus the cosmology could not be determined at these points. This showcases
the disadvantage of compactifying the phase space.

It would be interesting to extend this work by writing the system of equations
in terms of redshift (or scale factor). This idea involves evolving the
resulting equations and using the definition of the variables to plot the
corresponding $f(R)$ as a function of the Ricci scalar. In other words, this
new reconstruction method could be used to find the functional forms of $f(R)$
models associated with different cosmography as well as comparing them to
$\Lambda$CDM and observational data. The present work provides a basis for
reconstructing backgrounds that may be used in studies of large-scale
structures, independent of gravitational theory.

## Acknowledgement

ASA acknowledges the financial support provided by the University Grants
Commission (UGC) through the Senior Research Fellowship (File No. 16-9 (June
2017)/2018 (NET/CSIR)) to carry out the research work. SC acknowledges funding
support from the NSRF via the Program Management Unit for Human Resources and
Institutional Development, Research and Innovation (Thailand) [grant number
B01F650006].

## References

## References

  * J. P. Ostriker and P. J. Steinhardt [1995] J. P. Ostriker and P. J. Steinhardt, Cosmic concordance (1995), arXiv:astro-ph/9505066 [astro-ph] .
  * Capozziello  _et al._ [2008] S. Capozziello, V. F. Cardone, and V. Salzano, Cosmography of $f({R})$ gravity, Physical Review D 78, 10.1103/physrevd.78.063504 (2008).
  * Chakraborty  _et al._ [2021a] S. Chakraborty, K. MacDevette, and P. Dunsby, Model independent approach to the study of $f({R})$ cosmologies with expansion histories close to $\Lambda$CDM, Physical Review D 103, 10.1103/physrevd.103.124040 (2021a).
  * Goliath and Ellis [1999] M. Goliath and G. F. R. Ellis, Homogeneous cosmologies with a cosmological constant, Physical Review D 60, 10.1103/physrevd.60.023502 (1999).
  * Dunsby and Luongo [2016] P. K. S. Dunsby and O. Luongo, On the theory and applications of modern cosmography, Int. J. Geom. Meth. Mod. Phys. 13, 1630002 (2016), arXiv:1511.06532 [gr-qc] .
  * Sotiriou and Faraoni [2010] T. P. Sotiriou and V. Faraoni, $f(R)$ theories of gravity, Rev. Mod. Phys. 82, 451 (2010), arXiv:0805.1726 [gr-qc] .
  * De Felice and Tsujikawa [2010] A. De Felice and S. Tsujikawa, $f(R)$ theories, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Carloni  _et al._ [2005] S. Carloni, P. K. S. Dunsby, S. Capozziello, and A. Troisi, Cosmological dynamics of $R^{n}$ gravity, CQG 22, 4839 (2005).
  * Carloni  _et al._ [2009] S. Carloni, A. Troisi, and P. K. S. Dunsby, Some remarks on the dynamical systems approach to fourth order gravity, Gen. Rel. Grav. 41, 1757 (2009), arXiv:0706.0452 [gr-qc] .
  * Abdelwahab  _et al._ [2008] M. Abdelwahab, S. Carloni, and P. K. S. Dunsby, Cosmological dynamics of exponential gravity, Class. Quant. Grav. 25, 135002 (2008), arXiv:0706.1375 [gr-qc] .
  * Amendola  _et al._ [2007] L. Amendola, R. Gannouji, D. Polarski, and S. Tsujikawa, Conditions for the cosmological viability of $f(R)$ dark energy models, Phys. Rev. D 75, 083504 (2007), arXiv:gr-qc/0612180 .
  * Chakraborty  _et al._ [2021b] S. Chakraborty, K. MacDevette, and P. Dunsby, A model independent approach to the study of $f(R)$ cosmologies with expansion histories close to $\Lambda$CDM, Phys. Rev. D 103, 124040 (2021b), arXiv:2103.02274 [gr-qc] .
  * Dunajski and Gibbons [2008] M. Dunajski and G. Gibbons, Cosmic Jerk, Snap and Beyond, Class. Quant. Grav. 25, 235012 (2008), arXiv:0807.0207 [gr-qc] .
  * Di Valentino  _et al._ [2019] E. Di Valentino, A. Melchiorri, and J. Silk, Planck evidence for a closed Universe and a possible crisis for cosmology, Nature Astron. 4, 196 (2019), arXiv:1911.02087 [astro-ph.CO] .
  * Goheer  _et al._ [2008] N. Goheer, J. A. Leach, and P. K. S. Dunsby, Compactifying the state space for alternative theories of gravity, Class. Quant. Grav. 25, 035013 (2008), arXiv:0710.0819 [gr-qc] .
  * Abdelwahab  _et al._ [2012] M. Abdelwahab, R. Goswami, and P. K. S. Dunsby, Cosmological dynamics of fourth order gravity: A compact view, Phys. Rev. D 85, 083511 (2012), arXiv:1111.0171 [gr-qc] .
  * Kandhai and Dunsby [2015] S. Kandhai and P. K. S. Dunsby, Cosmological dynamics of viable $f(R)$ theories of gravity (2015), arXiv:1511.00101 [gr-qc] .
  * MacDevette  _et al._ [2022] K. MacDevette, P. Dunsby, and S. Chakraborty, Comprehensive analysis of the compact phase space for Hu-Sawicki $f(R)$ dark energy models including spatial curvature, Phys. Rev. D 106, 103533 (2022), arXiv:2208.15002 [gr-qc] .
  * Chakraborty  _et al._ [2022] S. Chakraborty, P. K. S. Dunsby, and K. Macdevette, A note on the dynamical system formulations in $f(R)$ gravity, IJGMMP 19, 10.1142/s0219887822300033 (2022).
  * Arora  _et al._ [2022] S. Arora, S. Mandal, S. Chakraborty, G. Leon, and P. K. Sahoo, Can $f(R)$ gravity isotropise a pre-bounce contracting universe?, JCAP 2022 (09), 042, arXiv:2207.08479 [gr-qc] .
  * Bouhmadi-Lopez  _et al._ [2013] M. Bouhmadi-Lopez, J. Morais, and A. B. Henriques, Smoking guns of a bounce in modified theories of gravity through the spectrum of gravitational waves, Phys. Rev. D 87, 103528 (2013), arXiv:1210.1761 [astro-ph.CO] .
  * Ellis  _et al._ [2016] G. F. R. Ellis, E. Platts, D. Sloan, and A. Weltman, Current observations with a decaying cosmological constant allow for chaotic cyclic cosmology, JCAP 2016 (04), 026, arXiv:1511.03076 [astro-ph.CO] .
  * Pavlović and Sossich [2021] P. Pavlović and M. Sossich, Dynamic properties of cyclic cosmologies, Phys. Rev. D 103, 023529 (2021), arXiv:2009.03625 [gr-qc] .