# A comprehensive analysis of the compact phase space for Hu-Sawicki $f(R)$
dark energy models including spatial curvature

Kelly MacDevette<EMAIL_ADDRESS>Department of Mathematics and Applied
Mathematics, University of Cape Town, Rondebosch 7700, Cape Town, South Africa
Peter Dunsby<EMAIL_ADDRESS>Department of Mathematics and Applied
Mathematics, University of Cape Town, Rondebosch 7700, Cape Town, South Africa
Saikat Chakraborty<EMAIL_ADDRESS>Center for Space Research, North-
West University, Mahikeng 2745, South Africa

(August 27, 2024)

###### Abstract

We present a comprehensive dynamical systems analysis of homogeneous and
isotropic Friedmann-Laîmatre-Robertson-Walker cosmologies in the Hu-Sawicki
$f(R)$ dark energy model for the parameter choice $\\{n,C_{1}\\}=\\{1,1\\}$.
For a generic $f(R)$ theory, we outline the procedures of compactification of
the phase space, which in general is 4-dimensional. We also outline how, given
an $f(R)$ model, one can determine the coordinate of the phase space point
that corresponds to the present day universe and the equation of a surface in
the phase space that represents the $\Lambda$CDM evolution history. Next, we
apply these procedures to the Hu-Sawicki model under consideration. We
identify some novel features of the phase space of the model such as the
existence of invariant submanifolds and 2-dimensional sheets of fixed points.
We determine the physically viable region of the phase space, the fixed point
corresponding to possible matter dominated epochs and discuss the possibility
of a non-singular bounce, re-collapse and cyclic evolution. We also provide a
numerical analysis comparing the $\Lambda$CDM evolution and the Hu-Sawicki
evolution.

##  I Introduction

A hundred years after it was formulated, Einstein’s theory of General
Relativity (GR) remains the best description of the gravitational action,
surviving observational tests ranging from the solar system to extragalactic
scales. Despite this success, high precision astrophysical and cosmological
observations conducted over the past two decades (supernovae Ia, cosmic
microwave background anisotropies, large scale structure formation, baryon
oscillations and weak lensing measurements) [1,
*SupernovaCosmologyProject:1998vns, *SupernovaSearchTeam:1998fmf,
*SupernovaSearchTeam:2003cyd, *SupernovaCosmologyProject:2003dcn,
*SupernovaSearchTeam:2004lze, *SNLS:2005qlf, *WMAP:2003elm, *WMAP:2006bqn,
*SDSS:2003eyi, *SDSS:2004kqt, *2dFGRS:2005yhx, *SDSS:2005xqv, *Blake:2005jd,
*Jain:2003tba] seems to suggest that, at least on cosmological scales, General
Relativity might not represent a complete description of gravity. In
particular, the late-time acceleration of the Universe cannot be described
within the framework of General Relativity without introducing additional
exotic matter degrees of freedom - known as Dark Energy. Although inferred by
fitting observational data to our preferred geometry of the universe (the
Robertson-Walker metric), Dark Energy has not been directly observed.
Currently the Concordance (or $\Lambda$CDM) model is regarded as the
theoretical framework that best fits all available data. Here a Cosmological
Constant dominates the present day universe, driving late-time acceleration,
while ordinary matter is dominated by Cold Dark Matter (CDM), which form the
potential wells for large-scale structure formation. Despite its success, the
$\Lambda$CDM model is plagued by significant fine-tuning problems related to
the vacuum energy scale, so it is important to consider other theoretical
alternatives, which are able to describe the main features of the universe,
without the introduction of a cosmological constant or dynamical Dark Energy.

One of the most popular alternatives to the Concordance Model is based on
modifications of standard Einstein gravity [16, 17]. Such models became
popular in the 1980’s because it was shown that they naturally admit a phase
of accelerated expansion which could be associated with inflation in the early
universe. This property has led to the idea that late-time acceleration could
also have a geometrical origin and that there is a connection between Dark
Energy and a non-standard behaviour of gravitation on cosmological scales.

Among one of the most widely studied of these geometrical dark energy models
is the so called $f(R)$ theories of gravity [18]. These models are based on
gravitational actions which are non-linear in the Ricci curvature R. Here the
higher order curvature corrections are written as an energy-momentum tensor of
geometrical origin describing an “effective” source term on the right hand
side of the standard Einstein field equations. This leads naturally to a
geometrical origin of the late-time acceleration of the Universe. $f(R)$
gravity is characterised by the existence of a propagating scalar degree of
freedom (called scalaron or curvaton) $\varphi=f^{\prime}(R)$ as apparent from
the trace field equation

$Rf^{\prime}(R)-2f(R)+3\Box f^{\prime}(R)=T,$ (1)

$T$ being the trace of the energy momentum tensor [16, 17]. GR is the trivial
case of $f(R)$ for which $f^{\prime}(R)=1$ and the scalar degree of freedom
vanishes. One of the well-known late-time $f(R)$ models is the so-called Hu-
Sawicki model [19], which, in its most generic form, is a 3-parameter $f(R)$
model. We will be particularly focusing on this important model in this paper.

Because of the fourth order nature of these theories, the cosmological
dynamics is considerably richer (and more complicated) than their GR
counterparts, making the theory of dynamical systems (DS) ideally suited for
studying them. In fact, describing these cosmologies using the dynamical
systems approach provides a relatively simple method for finding exact
solutions (via their corresponding fixed points) and obtaining a (qualitative)
description of the global dynamics of these models [20, 21, 22, 23, 24, 25,
26]. This approach has been very successful in shedding light on the dynamics
of a range of homogeneous cosmologies including the Friedmann-Laîmatre-
Robertson-Walker (FLRW), Bianchi I & 5 [27, 25] and Kantowski-Sachs models
[28]. It has also proved useful in identifying those orbits in phase-space
(expansion histories) which are close to the $\Lambda$CDM model, which can
then be used as suitable background cosmologies for the study of the growth of
large scale structure [29, 30, 31, 32].

The paper is organised as follows: In Sec.II we consider a generic $f(R)$
theory and review the DS approach to studying these models. In particular, we
describe how to compactify the phase-space and identify surfaces containing
orbits which represent $\Lambda$CDM evolution, as well as the point in phase-
space corresponding to the present universe (based on the current values of
the Hubble and deceleration parameters). In III we apply this framework to the
Hu-Sawicki $f(R)$ dark energy model for the parameter choice
$\\{n,C_{1}\\}=\\{1,1\\}$, obtaining a complete description of the dynamics of
these models including spatial curvature. We determine the physically viable
region of the phase space, the fixed point corresponding to possible matter
dominated epochs and discuss the possibility of a non-singular bounce, re-
collapse and cyclic evolution. We also provide a numerical analysis comparing
the $\Lambda$CDM evolution and the Hu-Sawicki evolution, setting initial
conditions so that the HS model is dynamically indistinguishable from the
Concordance model. Throughout the paper we use the metric signature
$(-,+,+,+)$ and the convention $\kappa=8\pi G=1$.

##  II Dynamical systems approach for generic $f(R)$ cosmology

We start by writing the field equation for $f(R)$ gravity in the presence of a
perfect fluid in a homogeneous and isotropic FLRW background [16, 17]

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

where the prime denotes a derivative with respect to $R$ and we have defined
the curvaton energy density and pressure as

$\displaystyle\rho_{R}\equiv\frac{1}{2}(Rf^{\prime}-f)-3H\dot{f^{\prime}},$
(3a) $\displaystyle
P_{R}\equiv\ddot{f^{\prime}}+2H\dot{f^{\prime}}-\frac{1}{2}(Rf^{\prime}-f).$
(3b)

The curvaton equation of state parameter is

$w_{R}\equiv\frac{P_{R}}{\rho_{R}}=\frac{\ddot{f^{\prime}}+2H\dot{f^{\prime}}-\frac{1}{2}(Rf^{\prime}-f)}{\frac{1}{2}(Rf^{\prime}-f)-3H\dot{f^{\prime}}},$
(4)

and the effective equation of state parameter of the universe is

$w_{\rm eff}\equiv\frac{P_{\rm eff}}{\rho_{\rm
eff}}=\frac{P+P_{R}}{\rho+\rho_{R}}=-\frac{2\dot{H}+3H^{2}+k/a^{2}}{3\left(H^{2}+k/a^{2}\right)}.$
(5)

If the perfect fluid is a barotropic one with an equation of state parameter
$\frac{P}{\rho}=w$, then $w,\,w_{R}$ and $w_{\rm eff}$ are related as

$w_{\rm eff}=w\frac{\rho}{\rho_{\rm eff}}+w_{R}\frac{\rho_{R}}{\rho_{\rm
eff}}.$ (6)

Since matter is minimally coupled to curvature, the perfect fluid and the
curvaton evolve independently according to their own continuity equations [33]

$\displaystyle\dot{\rho}+3H(1+w)\rho=0,$ (7a)
$\displaystyle\dot{\rho}_{R}+\left(3H(1+w_{R})+\frac{\dot{R}f^{\prime\prime}}{f^{\prime
2}}\right)\rho_{R}=-\left(2H-3\frac{\dot{R}f^{\prime\prime}}{f^{\prime}}\right)\frac{k}{a^{2}}+3\frac{H^{2}\dot{R}f^{\prime\prime}}{f^{\prime}}.$
(7b)

The matter continuity equation gives the usual scaling

$\displaystyle\rho=3H_{0}^{2}\Omega_{0}a^{-3(1+w)},$ (8)

where we have used the definition of the observable quantity matter density
parameter $\Omega_{0}=\frac{\rho_{0}}{3H_{0}^{2}}$, with the subscript ‘$0$’
being used to denote the present day values.

It is worth mentioning here two important conditions for physical viability of
any $f(R)$ gravity model:

  * •

$f^{\prime}(R)<0$ makes the scalar degree of freedom appearing in the theory a
ghost. To eradicate the possibility of a ghost degree of freedom, one must
require $f^{\prime}(R)>0$ for all $R$.

  * •

$f^{\prime\prime}(R)<0$ is related to unstable growth of curvature
perturbation in the weak gravity limit (this is also known as the Dolgov-
Kawasaki instability, as this was first suggested by Dolgov and Kawasaki in
[34]). Therefore one requires $f^{\prime\prime}(R)>0$ at least during the
early epoch of matter domination.

###  II.1 Phase space analysis for $f(R)$ gravity

In this section we present a dynamical system formulation for FLRW cosmologies
in $f(R)$ gravity following the Refs.[21, 22]. We start by writing the
Friedmann equation (2a) as follows

$\displaystyle
H^{2}=\frac{1}{3}\frac{\rho}{f^{\prime}}+\frac{1}{6f^{\prime}}(Rf^{\prime}-f)-H\frac{\dot{f}^{\prime}}{f^{\prime}}-\frac{k}{a^{2}}.$
(9)

Next, we divide both side by $3H^{2}$ and define a set of Hubble-normalized
dimensionless dynamical variables

$\tilde{x}=\frac{\dot{f}^{\prime}}{f^{\prime}}\frac{1}{H},\qquad\tilde{v}=\frac{1}{6}\frac{R}{H^{2}},\qquad\tilde{y}=\frac{1}{6}\frac{f}{f^{\prime}}\frac{1}{H^{2}},\qquad\tilde{\Omega}=\frac{1}{3}\frac{\rho}{f^{\prime}}\frac{1}{H^{2}},\qquad\tilde{K}=\frac{k}{a^{2}}\frac{1}{H^{2}},$
(10)

that are related by the constraint equation

$\displaystyle\tilde{\Omega}+\tilde{v}-\tilde{y}-\tilde{x}-\tilde{K}=1.$ (11)

Differentiating the dynamical variables with respect to the e-folding number
$N\equiv\ln a$, one obtains the following autonomous system of equations

$\displaystyle\frac{d\tilde{x}}{dN}=-2\tilde{v}+4+3\tilde{x}-3(1+w)\tilde{\Omega}+4\tilde{K}-\tilde{x}^{2}-\tilde{x}\tilde{v}+\tilde{x}\tilde{K},$
(12a)
$\displaystyle\frac{d\tilde{v}}{dN}=\tilde{v}\left(\Gamma\tilde{x}-2\tilde{v}+2\tilde{K}+4\right),$
(12b)
$\displaystyle\frac{d\tilde{K}}{dN}=-2\tilde{K}\left(\tilde{v}-\tilde{K}-1\right),$
(12c)
$\displaystyle\frac{d\tilde{\Omega}}{dN}=-\tilde{\Omega}\left(-1+3w+\tilde{x}+2\tilde{v}-2\tilde{K}\right),$
(12d)
$\displaystyle\frac{d\tilde{y}}{dN}=\Gamma\tilde{x}\tilde{v}+\tilde{y}(2\tilde{K}-2\tilde{v}-\tilde{x}+4),$
(12e)

where we had defined an auxiliary quantity

$\Gamma\equiv\frac{f^{\prime}}{Rf^{\prime\prime}},$ (13)

that is explicitly dependent on the functional form of $f(R)$. The system of
dynamical equations (12), along with the Friedmann constraint (11),
constitutes the dynamical system for FLRW cosmologies in $f(R)$ gravity.
Because of the existence of the constraint, one of the dynamical variables is
redundant, and the resulting phase space is actually 4-dimensional.

In order to close the system, $\Gamma$ must be expressed as a function of the
dynamical variables. This can be done by noting that the relation

$\frac{\tilde{v}}{\tilde{y}}=\frac{Rf^{\prime}}{f}$ (14)

can, in principle, be inverted to obtain
$R=R\left(\frac{\tilde{v}}{\tilde{y}}\right)$. Using this one can obtain
$\Gamma=\Gamma\left(\frac{\tilde{v}}{\tilde{y}}\right)$. As mentioned in [22],
the success of this particular dynamical system formulation crucially depends
on the invertibility of the relation (14). This is a limitation of this
formulation. For example, the most generic case of the Hu-Sawicki $f(R)$ model
[19] cannot be treated with this approach and a different formulation is
necessary. A possible solution to this issue was devised in [23]. In this
paper we will focus on the only particular case of the Hu-Sawicki $f(R)$ that
can in fact be treated with this formulation [35] 111For a historical account
of different dynamical system formulations for FLRW cosmology in $f(R)$
gravity, see [36]..

The Hubble-normalized dynamical variables defined in Eq.(10) are non-compact,
i.e., their domain is infinite or semi-infinite. In particular, by definition,
all of them diverge at a cosmological bounce or recollapse, that necessarily
requires the condition $H=0$. Therefore, there can be interesting physical
properties at the infinity of the phase space that one misses using this
formulation. Also, one can justifiably take the e-folding $N\equiv\ln a$ as a
phase space time variable as long as one confines the attention to an
expanding universe only, as $\dot{N}=H$ is always positive. For a contracting
universe one should alter the definition as $N\equiv-\ln a$ and in general one
should define the phase space time variable as $N\equiv\frac{|H|}{H}\ln a$.
Clearly this definition is discontinuous at a cosmological bounce or
recollapse. As we will see in the next subsection, this shortcoming is not
present in the compact phase space formulation.

###  II.2 Compactifying the phase space

In this section we outline a prescription for the compactification of the
phase space. We follow the approach of Refs.[24, 35], but generalize it to the
case of non-negative global spatial curvature 222For a different
compactification approach via the Poincare compactification, see [20]. We
start by writing the Friedmann equation (2a) as follows

$D^{2}\equiv\left(3H+\frac{3}{2}\frac{\dot{f^{\prime}}}{f^{\prime}}\right)^{2}+\frac{3}{2}\left(\frac{f}{f^{\prime}}+\frac{6k}{a^{2}}\right)=\frac{3\rho}{f^{\prime}}+\frac{3}{2}R+\frac{9}{4}\left(\frac{\dot{f^{\prime}}}{f^{\prime}}\right)^{2}.$
(15)

Next, we divide both side by $D^{2}$ and define a set of _compact_ dynamical
variables

$x=\frac{3}{2}\frac{\dot{f}^{\prime}}{f^{\prime}}\frac{1}{D},\qquad
v=\frac{3}{2}\frac{R}{D^{2}},\qquad
y=\frac{3}{2}\frac{f}{f^{\prime}}\frac{1}{D^{2}},\qquad\Omega=\frac{3\rho}{f^{\prime}}\frac{1}{D^{2}},\qquad
Q=\frac{3H}{D},\qquad K=\frac{9k}{a^{2}}\frac{1}{D^{2}},$ (16)

which are now related by two constrained equations

$\displaystyle(Q+x)^{2}+K+y=1,$ (17a) $\displaystyle\Omega+x^{2}+v=1,$ (17b)

that comes directly from Eq.(15). The dynamical variables defined in Eq.(16)
are compact because, requiring $f^{\prime}(R)>0$, restricting ourselves to
non-negative global spatial curvature and keeping in mind the two constraints
(17a),(17b), one gets the following finite domains for the dynamical variables

$\displaystyle-1\leq x\leq 1,\quad 0\leq\Omega\leq 1,\quad-2\leq Q\leq 2,$
(18) $\displaystyle 0\leq v\leq 1,\quad 0\leq y\leq 1,\quad 0\leq K\leq 1.$

It is straightforward to find a relationship between the compact dynamical
variables and their non-compact counterparts

$\tilde{x}=2\frac{x}{Q},\qquad\tilde{v}=\frac{v}{Q^{2}},\qquad\tilde{y}=\frac{y}{Q^{2}},\qquad\tilde{K}=\frac{K}{Q^{2}},\qquad\tilde{\Omega}=\frac{\Omega}{Q^{2}}.$
(19)

To write an autonomous dynamical system, we need to define a new phase space
time variable $\tau$ such that

$d\tau\equiv Ddt,$ (20)

Differentiating the compact dynamical variables with respect to $\tau$ and
using the field equations results in a system of six first order autonomous
differential equations. Because of the existence of two constraint equations
(17a) and (17b), two of the dynamical variables are redundant. Therefore, the
compact phase space, just like it’s non-compact counterpart, is actually
4-dimensional. We choose to eliminate $y$ and $\Omega$ and write the dynamical
system as follows

$\displaystyle\frac{dv}{d\tau}=-\frac{1}{3}v\left((Q+x)\left(2v-(1+3w)(1-x^{2}-v)+4xQ\right)-2Q-4x+2x\Gamma(v-1)+4xK\right),$
(21a)
$\displaystyle\frac{dx}{d\tau}=\frac{1}{6}\Bigg{[}-2x^{2}v\Gamma+(1-3w)(1-x^{2}-v)+2v+4\left(x^{2}-1\right)\left(1-Q^{2}-xQ\right)+x(Q+x)\left((1+3w)(1-x^{2}-v)-2v\right)$
$\displaystyle\qquad+4K(1-x^{2})\Bigg{]},$ (21b)
$\displaystyle\frac{dQ}{d\tau}=\frac{1}{6}\Bigg{[}-4xQ^{3}+xQ(5+3w)(1-xQ)-Q^{2}(1-3w)-Qx^{3}(1+3w)-3vQ(1+w)(Q+x)$
$\displaystyle\qquad+2v(1-Qx\Gamma)-2K(1+2xQ)\Bigg{]},$ (21c)
$\displaystyle\frac{dK}{d\tau}=-\frac{1}{3}K\left((Q+x)(-(1+3w)(1-x^{2}-v)+4xQ+2v)+4x(K-1)+2xv\Gamma\right).$
(21d)

The system can be closed by inverting the relation

$\frac{v}{y}=\frac{Rf^{\prime}}{f},$ (22)

to find $R=R\left(\frac{v}{y}\right)$ and using it to write
$\Gamma=\Gamma\left(\frac{v}{y}\right)$, provided, of course, that the
relation (22) is invertible. Note that, unlike in the non-compact phase space
analysis, the definition of the phase space time variable $\tau$ is _not_
discontinuous at a cosmological bounce or recollapse, as $\dot{\tau}=D$ is
always positive definite by definition.

A compact phase space analysis helps analyze the global features of the
solution space of a theory, some of which might not be apparent from a non-
compact phase space analysis. From this point onwards, we will use only the
compact phase space analysis for all of our subsequent analysis.

###  II.3 $\Lambda$CDM point

In this section we will try to determine which point in the 4-dimensional
phase space $v$-$x$-$Q$-$K$ corresponds to the present day universe. In other
words, we will outline how to calculate the values of the compact dynamical
variables $v,\,x,\,Q,\,K$ based on the observed values of different
cosmological parameters. To this aim, let us first introduce two sets of
dimensionless parameters. The first set is the usual matter density parameter
($\Omega_{m}$) and spatial curvature density parameter ($\Omega_{k}$)

$\Omega_{m}\equiv\frac{\rho}{3H^{2}},\qquad\Omega_{k}\equiv\frac{k}{a^{2}H^{2}},$
(23)

where we have deliberately used the subscript ‘$m$’ in $\Omega_{m}$ to make it
apparent that it is different from the compact dynamical variable $\Omega$
defined in Eq.(16) ($\Omega\rightarrow\Omega_{m}$ only at the GR limit
$f^{\prime}\rightarrow 1$). The second set is the dimensionless cosmographic
parameters, namely the dimensionless Hubble parameter ($h$), deceleration
parameter ($q$), jerk parameter ($j$) and snap parameter ($s$) parameter [37,
38]

$h\equiv\frac{H}{H_{0}},\qquad q\equiv-\frac{\ddot{a}}{aH^{2}},\qquad
j\equiv\frac{\dddot{a}}{aH^{3}},\qquad s\equiv\frac{a^{(4)}}{aH^{4}},$ (24)

where $H_{0}$ is the present day value of the Hubble parameter. Both the
matter density parameters and the cosmographic parameters are sets of
observable quantities.

Assuming that the dominant hydrodynamic fluid content in the universe is given
by the cold dark matter (CDM) that is modelled by a dust fluid, the Hubble
parameter and its first derivative are given by

$\displaystyle
H^{2}=\frac{\rho}{3f^{\prime}}+\frac{1}{6f^{\prime}}(Rf^{\prime}-f)-H\frac{\dot{f}^{\prime}}{f^{\prime}}-\frac{k}{a^{2}},$
(25a)
$\displaystyle\dot{H}=-\frac{1}{2f^{\prime}}\left[\frac{1}{3}\rho+\ddot{R}f^{\prime\prime}+\dot{R}^{2}f^{\prime\prime\prime}-Hf^{\prime\prime}\dot{R}\right]+\frac{3}{2}\frac{k}{a^{2}},$
(25b)

It can be verified that Eqs.((25a),(25b)) are completely equivalent to
Eqs.((2a),(2b)). The Ricci scalar and its first two derivatives are as follows

$\displaystyle R=6\left(\dot{H}+2H^{2}+\frac{k}{a^{2}}\right),$ (26a)
$\displaystyle\dot{R}=6\left(\ddot{H}+4H\dot{H}-2H\frac{k}{a^{2}}\right),$
(26b)
$\displaystyle\ddot{R}=6\left[\dddot{H}+4H\ddot{H}+4\dot{H}^{2}+\left(4H^{2}-2\dot{H}\right)\frac{k}{a^{2}}\right].$
(26c)

Using the matter abundance parameters and the cosmographic parameters, the
Ricci scalar and its first two derivatives can be expressed as purely
algebraic combinations of them [33]

$\displaystyle\frac{R}{H_{0}^{2}}=6h^{2}(1-q+\Omega_{k}),$ (27a)
$\displaystyle\frac{\dot{R}}{H_{0}^{3}}=6h^{3}(j-q-2-2\Omega_{k}),$ (27b)
$\displaystyle\frac{\ddot{R}}{H_{0}^{4}}=6h^{4}\left(s+q^{2}+8q+6+2(3+q)\Omega_{k}\right).$
(27c)

Eqs.((25a),(25b)) can also be expressed using the cosmographic parameters as
purely algebraic equations i.e. not explicitly containing any time derivatives

$\displaystyle\mathcal{A}f^{\prime}+\frac{f}{H^{2}}+\mathcal{B}\left(H^{2}f^{\prime\prime}\right)-6\Omega_{m}=0,$
(28a)
$\displaystyle\mathcal{C}f^{\prime}-\mathcal{D}\left(H^{2}f^{\prime\prime}\right)-\mathcal{G}\left(H^{4}f^{\prime\prime\prime}\right)-\Omega_{m}=0,$
(28b)

with the coefficients
$\mathcal{A},\,\mathcal{B},\,\mathcal{C},\,\mathcal{D},\,\mathcal{G}$ being
different algebraic combinations of the cosmographic parameters and the
spatial curvature density parameter

$\displaystyle\mathcal{A}=6q,$ (29a)
$\displaystyle\mathcal{B}=36(j-q-2-2\Omega_{k}),$ (29b)
$\displaystyle\mathcal{C}=3\Omega_{k}+2+2q,$ (29c)
$\displaystyle\mathcal{D}=6(8-j+9q+q^{2}+2\Omega_{k}(4+q)+s),$ (29d)
$\displaystyle\mathcal{G}=36(j-q-2-2\Omega_{k})^{2}.$ (29e)

Henceforth we use the subscript ‘$0$’ to denote the parameter values in the
current epoch. Since the value of the scale factor itself bears no dynamical
significance, it is conventional to take $a_{0}=1$. Also, by definition, the
present day value of the dimensionless Hubble parameter is $h_{0}=1$. Other
cosmographic parameters are quantities that can be calculated from the
observational data once a parametrization for the dark energy equation of
state is fixed. Using the Chevallier-Polarski-Linder (CPL) parameterization
[39, 40] for the dark energy equation of state

$w_{\rm DE}=w_{0}+w_{a}(1-a)$ (30)

and fixing the free parameters to those corresponding to the $\Lambda$CDM
model $(w_{0},w_{a})=(-1,0)$ results in the following expressions for
$(q_{0},j_{0},s_{0})$

$\displaystyle q_{0}=\frac{1}{2}-\frac{3}{2}\Omega_{\Lambda 0},$ (31a)
$\displaystyle j_{0}=\Omega_{k}+1,$ (31b) $\displaystyle
s_{0}=1-\frac{9}{2}\Omega_{m0},$ (31c)

where $\Omega_{\Lambda}$ is used to denote the dark energy density parameter.
The values of different density parameters, taken from Planck 2018 data [41],
are

$(\Omega_{\Lambda 0},\Omega_{k0},\Omega_{m0})=(0.6843,0.0007,0.315).$ (32)

This gives the present day values for the cosmographic parameters within the
$\Lambda$CDM model as

$(q_{0},j_{0},s_{0})=(-0.52645,1.0007,-0.4175).$ (33)

The values in Eqs.((32),(33)) allow us to calculate the present day values of
the Ricci scalar and its derivatives using Eq.(27)

$\left(\frac{R_{0}}{H_{0}^{2}},\frac{\dot{R}_{0}}{H_{0}^{3}},\frac{\ddot{R}_{0}}{H_{0}^{4}}\right)=(9.1629,-2.8455,9.9091),$
(34)

as well as the present day values of the coefficients in (29)

$(\mathcal{A}_{0},\mathcal{B}_{0},\mathcal{C}_{0},\mathcal{D}_{0},\mathcal{G}_{0})=(-3.1587,-17.0730,0.9492,12.7546,8.0969).$
(35)

With the known values of these coefficients, one can see now that Eq.(28)
becomes a pair of algebraic equations with constant coefficients for
$\frac{f_{0}}{H_{0}^{2}},f^{\prime}_{0},H_{0}^{2}f^{\prime\prime}_{0},H_{0}^{4}f^{\prime\prime\prime}_{0}$.
If one now specifies an explicit functional form of $f(R)$, the two equations
(28a) and (28b) can be used to determine the value of up to two a-priori
unspecified model parameter of the theory (modulo a factor of $H_{0}$ that is
used to make the model parameters dimensionless).

Note that if there are more than two unspecified model parameters, one
requires more equations additional to (28a),(28b), which are obtained by
taking time derivatives of Eq.(25b) and considering cosmographic parameters of
order higher than the snap parameter. On the other hand, if there is only one
unspecified model parameter, then the two equations (28a),(28b) are redundant
to find its value. In this situation, one should consider the reduced set of
cosmographic parameters $\\{q,j\\}$ and use their present values to determine
the value of the model parameter. Eq.(28a) now suffices for this purpose.
Eq.(28b) in this case does not provide an additional equation to determine the
model parameter, as it contains $s$ whose present day value will not be
required anymore.

Once the value(s) of the a-priori unspecified model parameter(s) is determined
(modulo the factor of $H_{0}$), one can calculate the values of the compact
dynamical variables corresponding to the present day observable universe using
the values given in Eqs.((33),(34)). Henceforth we will call this point in the
phase space the $\Lambda$CDM point. In Sec.III.8 we will calculate the phase
space coordinates of the $\Lambda$CDM point for the specific case of the Hu-
Sawicki $f(R)$ model that we consider in this paper.

###  II.4 $\Lambda$CDM surface

The cosmological evolution in the $\Lambda$CDM model can be written as

$\displaystyle h(a)$ $\displaystyle=$
$\displaystyle\sqrt{\Omega_{m0}a^{-3}+\Omega_{\Lambda 0}-\Omega_{k0}a^{-2}}$
(36) $\displaystyle=$ $\displaystyle\sqrt{\Omega_{m0}a^{-3}+\Omega_{\Lambda
0}-(\Omega_{m0}+\Omega_{\Lambda 0}-1)a^{-2}},$

where $\Omega_{m0},\,\Omega_{\Lambda 0},\,\Omega_{k0}$ has the values
specified in (32) and the second line comes from the first line by writing the
equation at the present epoch. The $\Lambda$CDM evolution history can be
specified by the necessary and sufficient condition333By necessary and
sufficient condition, we imply that Eq.(36), differentiated twice, gives
Eq.(37) and Eq.(37), integrated twice, gives back Eq.(36). [38]

$j=\Omega_{k}+1.$ (37)

This is a purely kinematic condition specifying _how_ the universe evolves,
with no information on what is the inherent dynamical model for the evolution.
There can be models other than the $\Lambda$CDM that obey the same
cosmographic requirement (37) and therefore produces a cosmological evolution
identical to (36). Such models are kinematically degenerate with $\Lambda$CDM
at the background level. For an $f(R)$ model to exactly mimic the $\Lambda$CDM
evolution history at the background level, it must obey the condition
(37)444For reconstruction of $f(R)$ models based on the $\Lambda$CDM evolution
history (36), see [42]. For a generic dynamical system analysis of
$\Lambda$CDM mimicking $f(R)$ models, see [43].

The deceleration parameter $q$ can be expressed in terms of the dynamical
variables as

$q=1+\tilde{K}-\tilde{v}=1+\frac{K-v}{Q^{2}}.$ (38)

The jerk parameter $j$ can be expressed in terms of the deceleration parameter
and it’s time derivative.

$j=2q^{2}+q-\frac{1}{H}\frac{dq}{dt}=2q^{2}+q-\frac{3}{Q}\frac{dq}{d\tau}.$
(39)

Using the dynamical equations (21), $j$ can be expressed as a function of the
dynamical variables

$j=j(K,x,v,Q),$ (40)

whose explicit functional form we do not write here as the relevant
calculations can be tackled numerically. Eq.(37) now provides an algebraic
constraint which specifies a 3-dimensional hypersurface over the entire
4-dimensional phase space $K$-$x$-$v$-$Q$. Henceforth we will call it the
$\Lambda$CDM surface. If an $f(R)$ theory is to produce a cosmic evolution
history that is asymptotically $\Lambda$CDM, then there must be a De-Sitter
saddle or stable fixed point lying on this hypersurface and a phase trajectory
that asymptotically approaches this point. If an $f(R)$ theory is to produce a
cosmic evolution history that is indistinguishable from the $\Lambda$CDM
evolution (36) at all times, then there must be an unstable or a saddle fixed
point corresponding to a matter dominated phase on this surface, a De-Sitter
saddle or stable fixed point on this surface, and a phase trajectory
connecting these two fixed points it must produce a phase trajectory lying
entirely on this surface connecting the two fixed points.

##  III Dynamical systems analysis applied to the Hu-Sawicki $f(R)$ Model

A theoretically viable $f(R)$ alternative to dark energy must satisfy the
following criteria

  * •

$f(R)\to R$ for $R\to 0$.

  * •

$f(R)\to R-2\Lambda$ ($\Lambda>0$) for $R\to\infty$.

  * •

$f^{\prime}(R)>0$ for all $R$ so that the curvature degree of freedom is
always physical.

  * •

$f^{\prime\prime}(R)>0$ for all $R$ so that at no point of time there is an
unbounded growth of curvature perturbation.

The Hu-Sawicki $f(R)$ model [19] was specifically designed so as to satisfy
all of these conditions. In its most generic form, the Hu-Sawicki (henceforth
HS) model is a 3-parameter theory which has the following form

$f(R)=R-\frac{C_{1}R^{n}}{C_{2}R^{n}+1}.$ (41)

Taking the limit $R\to\infty$, one can read off $\Lambda=\frac{C_{1}}{C_{2}}$,
so that $\frac{C_{1}}{C_{2}}$ must be positive. As was pointed out in
Ref.[35], only the particular case $n=1=C_{1}$ can be studied with the present
dynamical system formulation since only for this special case the relation
(22) is invertible. Henceforth we will confine our study to this special case
only. Choosing to work with the particular case $n=1$ does not, in any way,
make our subsequent analysis any less relevant, because it was shown in
Ref.[44] that $n$ remains completely unconstrained by current cosmological
data. Since we live in a universe with a positive value of the cosmological
constant, choosing $C_{1}=1$ means now that $C_{2}$ must have a positive
value.

###  III.1 Compact dynamical system

For $n=1=C_{1}$, the relation (22) is invertible and one can find

$\displaystyle\Gamma=\frac{1}{2}\frac{vy}{(v-y)^{2}}.$ (42)

With the above expression for $\Gamma$, the dynamical system (21) becomes
singular with a pole of second order on the 3-dimensional hypersurface given
by

$\displaystyle z\equiv(v-y)=v+(Q+x)^{2}+K-1=0,$ (43)

which means one needs to be careful to construct a dynamical system here so as
to be able to apply the dynamical system technique. To this goal, we redefine
the phase space time variable as

$\displaystyle d\eta=\frac{d\tau}{\big{(}v+(Q+x)^{2}+K-1\big{)}^{2}},$ (44)

so that with respect to this redefined time variable, the new dynamical system
becomes regular everywhere. The explicit form of the dynamical system in the
compact phase space in presence of a generic perfect fluid and for the HS
model with $n=1=C_{1}$ is as follows

$\displaystyle\frac{dv}{d\eta}=-\frac{1}{3}vz^{2}\Big{(}(Q+x)\big{(}2v-(1+3w)(1-x^{2}-v)+4xQ\big{)}-2Q-4x+4xK\Big{)}-\frac{1}{3}v^{2}x(v-z)(v-1),$
(45a)
$\displaystyle\frac{dx}{d\eta}=\frac{1}{6}z^{2}\Bigg{[}(1-3w)(1-x^{2}-v)+2v+4(x^{2}-1)(1-Q^{2}-xQ)+x(Q+x)\big{(}(1+3w)(1-x^{2}-v)-2v\big{)}$
$\displaystyle+4K(1-x^{2})\Bigg{]}-\frac{1}{6}x^{2}v^{2}(v-z),$ (45b)
$\displaystyle\frac{dQ}{d\eta}=\frac{1}{6}z^{2}\Bigg{[}-4xQ^{3}+xQ(5+3w)(1-xQ)-Q^{2}(1-3w)-Qx^{3}(1+3w)-3vQ(1+w)(Q+x)$
$\displaystyle+2v-2K(1+2xQ)\Bigg{]}-\frac{1}{6}v^{2}(v-z)Qx,$ (45c)
$\displaystyle\frac{dK}{d\eta}=-\frac{1}{3}z^{2}K((Q+x)(-(1+3w)(1-x^{2}-v)+4xQ+2v)+4x(K-1))-\frac{1}{3}Kxv^{2}(v-z).$
(45d)

The variable $z$ from Eq.(43) follows the dynamical equation

$\displaystyle\frac{dz}{d\eta}=$
$\displaystyle\frac{1}{3}z\Bigg{[}v^{2}x\Big{(}k+(Q+x)^{2}+v-1\Big{)}-z\bigg{(}x\left(4k^{2}+k\left(8Q^{2}+7v-7\right)+4Q^{4}+Q^{2}(v-1)(9w+13)+3(v-1)^{2}(w+1)\right)$
(46)
$\displaystyle+Q\left(3k(v-1)+3Q^{2}(v-1)(w+1)+3v^{2}+3(v-1)^{2}w-4v+3\right)+x^{3}\Big{(}5k+3Q^{2}(3w+5)+2(v-1)(3w+2)\Big{)}$
$\displaystyle+Qx^{2}\left(13k+Q^{2}(3w+13)+2(v-1)(6w+7)\right)+Q(9w+7)x^{4}+(3w+1)x^{5}\bigg{)}+x\Big{(}k+(Q+x)^{2}-z-1\Big{)}^{3}\Bigg{]}.$

Interestingly, the singular submanifold $z=0$ now becomes an invariant
submanifold in the phase space, by virtue of the time redefinition (44).

The fixed points of the system for the case when the perfect fluid is CDM
($w=0$) are listed in Tab.1, along with their nature of stability.

Point | $(K,x,v,Q,\Omega,y)$ | Stability | Deceleration parameter | $w_{\rm eff}$ | Cosmology
---|---|---|---|---|---
$\mathcal{S}_{1\pm}$ | $\left(K,0,v,\pm\sqrt{1-K-v},1-v,v\right)$ | Saddle | $\frac{1-2v}{1-K-v}$ | $\frac{1}{3}\frac{(1-3v)}{(1-v)}$ | Depends on the point
$\mathcal{S}_{2\pm}$ | $\left(K,x,0,-x\pm\sqrt{1-K},1-x^{2},0\right)$ | $\begin{cases}&\text{Attractor for $3Q<x<0$},\\\ &\text{Repeller for $0<x<3Q$},\\\ &\text{Saddle otherwise}.\end{cases}$ | $\frac{x^{2}\mp 2x\sqrt{1-K}+1}{(-x\pm\sqrt{1-K})^{2}}$ | $\frac{1}{3}$ | Depends on the point
$\mathcal{E}_{\pm}$ | $\left(0,0,1,\pm\frac{1}{\sqrt{2}},0,\frac{1}{2}\right)$ | NH | $-1$ | $-1$ | De-Sitter
$\mathcal{F}$ | $\left(\frac{3}{5},0,\frac{3}{5},0,\frac{2}{5},\frac{2}{5}\right)$ | Spiral | undefined | $-\frac{1}{3}$ | Einstein static solution
$\mathcal{G}_{\pm}$ | $\left(0,\pm\frac{1}{3},\frac{8}{9},\pm\frac{2}{3},0,0\right)$ | Attractor/Repeller | $-1$ | $-1$ | De-Sitter
$\mathcal{H}_{\pm}$ | $\left(\frac{2}{3},\pm\frac{1}{\sqrt{3}},\frac{2}{3},0,0,0\right)$ | Saddle | undefined | $-\frac{1}{3}$ | Einstein static solution
$\mathcal{I}_{\pm}$ | $\left(\frac{3}{5},\mp\frac{1}{\sqrt{5}},\frac{4}{5},\pm\frac{1}{\sqrt{5}},0,\frac{2}{5}\right)$ | Saddle | $0$ | $-\frac{1}{3}$ | $a\sim t$
Table 1: Fixed points for the HS model with the parameter choice $n=C_{1}=1$
in presence of CDM ($w=0$). Note that $\mathcal{S}1_{\pm}$ and
$\mathcal{S}2_{\pm}$ are actually two pairs of sheets of fixed points.
Deceleration parameter $q=-1-\frac{\dot{H}}{H^{2}}$ and the effective equation
of state parameter $w_{\rm
eff}=-\frac{2\dot{H}+3H^{2}+\frac{k}{a^{2}}}{3\left(H^{2}+\frac{k}{a^{2}}\right)}$
can be determined in terms of the compact dynamical variables as
$q=1+\frac{K-v}{Q^{2}}$ and $w_{\rm eff}=\frac{Q^{2}+K-2v}{3(Q^{2}+K)}$.

As can be seen from the table, there are nine isolated fixed points
$\mathcal{E}_{\pm},\,\mathcal{F},\,\mathcal{G}_{\pm},\,\mathcal{H}_{\pm},\,\mathcal{I}_{\pm}$
and two 2-dimensional “sheets” of fixed points that can be expressed in a
compact way as

$\displaystyle\mathcal{S}_{1}\equiv(K,x,v,Q)=\left((1-Q^{2})\cos^{2}\theta,0,(1-Q^{2})\sin^{2}\theta,Q\right),$
(47)
$\displaystyle\mathcal{S}_{2}\equiv(K,x,v,Q)=\left(1-(Q+x)^{2},x,0,Q\right).$
(48)

Both $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ have one part corresponding to an
expanding phase ($Q>0$) and another part corresponding to a contracting phase
($Q<0$). In Tab.1 we explicitly write them as $\mathcal{S}_{1\pm}$ and
$\mathcal{S}_{2\pm}$. Stability of the isolated fixed points is found in the
usual way using the Hartman-Grobman theorem. Stability of the sheets of fixed
points need some further discussion which we provide in Appendix B. The entire
phase space is 4-dimensional.

The phase space of the HS model is very rich in structure and many of its
features deserve an entire subsection dedicated to its detailed discussion.

###  III.2 Sheets of fixed points

As mentioned previously, there are two pairs of sheets of fixed points
$\mathcal{S}1_{\pm}$ and $\mathcal{S}2_{\pm}$ in the 4-dimensional compact
phase space of the Hu-Sawicki model. Since all the points on these 2-D sheets
are fixed points, there is no phase flow either on the sheet or across the
sheet. Therefore each of these sheets of fixed points are also invariant
submanifolds. Each of the sheets $S1_{\pm}$ divides the 3-dimensional
submanifold $x=0$ into two disjoint volumes. Likewise, each of the sheets
$S2_{\pm}$ divides the 3-dimensional submanifold $v=0$ into two disjoint
volumes. It is to be kept in mind, however, that these sheets do _not_ divide
the entire 4-dimensional phase space into two disjoint 4-dimensional volumes,
as in that case one would require to have not a _sheet_ of fixed points but a
_volume_ of fixed points (as a lower dimensional analogy one could think of a
1-parameter family of fixed points in a complete 3-dimensional phase space.).

The Jacobian of the dynamical system calculated on all four sheets yields all
the four eigenvalues to be null. Therefore the Hartman-Grobman theorem cannot
be used to determine their stability. We expand on this in Appendix B.

###  III.3 Physically viable regions

Not the entire 4-dimensional phase space is physically viable. For the HS
model with $n=C_{1}=1$, one can obtain from the definitions of the phase space
variables $v$ and $y$ that

$C_{2}R+1=\frac{y}{v-y},$ (49)

and therefore

$\displaystyle
f^{\prime}(R)=1-\frac{1}{(C_{2}R+1)^{2}}=\frac{v(2y-v)}{y^{2}},$ (50a)
$\displaystyle
f^{\prime\prime}(R)=\frac{2C_{2}}{(C_{2}R+1)^{3}}=\frac{2C_{2}(v-y)^{3}}{y^{3}}.$
(50b)

Absence of ghost and tachyonic instability requires
$f^{\prime}>0,\,f^{\prime\prime}>0$. Keeping in mind that $C_{2}$ is a
positive constant, this constrains the physically viable region of the phase
space as

$\frac{1}{2}v<y<v,$ (51)

or, using the constraint equation (17a),

$\frac{1}{2}v<1-K-(Q+x)^{2}<v.$ (52)

Furthermore, one requires $0\leq y,\,\Omega\leq 1$, which, using the
constraint equations (17a) and (17b), translates into the conditions

$\displaystyle 0\leq(Q+x)^{2}+K\leq 1,$ (53) $\displaystyle 0\leq v+x^{2}\leq
1.$ (54)

Constraints (52),(53),(54) must be imposed additionally to single out the
physically viable region of the phase space.

###  III.4 Matter dominated epoch

An effectively matter dominated epoch is given by a fixed point which has
$w_{\rm eff}=0$. The only candidate for this is the 1-parameter family of
fixed points given by the curve

$\mathcal{M}\equiv(K,x,v,Q,\Omega,y)=\left(\frac{2}{3}-Q^{2},0,\frac{1}{3},Q,\frac{2}{3},\frac{1}{3}\right)\subset\mathcal{S}_{1}.$
(55)

If one confines the attention to the invariant submanifold $K=0$, i.e.
consider only the spatially flat cosmologies, then one gets two pairs of
isolated fixed points

$\mathcal{P}_{m\pm}\equiv(x,v,Q,\Omega,y)=\left(0,\frac{1}{3},\pm\sqrt{\frac{2}{3}},\frac{2}{3},\frac{1}{3}\right).$
(56)

which correspond to an effectively dark matter dominated epoch during the
expansion and contraction phase respectively. From Tab.1 one can note that
matter dominated fixed points are always saddle fixed points i.e. an
intermediate epoch of a cosmology.

###  III.5 Fixed points on the $\Lambda$CDM surface

It can be explicitly checked that the fixed points $\mathcal{E}_{\pm}$ and
$\mathcal{G}_{\pm}$ satisfy the cosmographic condition (37), i.e. they fall on
the $\Lambda$CDM surface. In an expanding universe, $\mathcal{E}_{+}$ is a
saddle (i.e. an intermediate phase) whereas $\mathcal{G}_{+}$ is a stable
fixed point (i.e. a future attractor). For the HS model under consideration to
asymptotically mimic the $\Lambda$CDM evolution at late times, the phase
trajectories near the matter dominated fixed point $\mathcal{M}$ must end at
the De-Sitter future attractor $\mathcal{G}_{+}$ via the intermediate De-
Sitter epoch $\mathcal{E}_{+}$.

Figure 1: A characteristic phase trajectory for the spatially flat case $K=0$
passing through the point $(x,v,Q)=\left(0,0.334,\sqrt{\frac{2}{3}}\right)$
(indicated by the red dot) in the vicinity of the matter dominated fixed point
$\mathcal{P}_{m+}\equiv(x,v,Q)=\left(0,\frac{1}{3},\sqrt{\frac{2}{3}}\right)$.
The trajectory passes near the intermediate De-Sitter phase $\mathcal{E}_{+}$
and ends at the De-Sitter future attractor $\mathcal{G}_{+}$. Both the De-
Sitter points lie on the $\Lambda$CDM surface.

Fig.1 shows such a trajectory for the spatially flat case. The trajectory is
plotted by numerically solving the system in eq.(45), setting
$\\{K,w\\}=\\{0,0\\}$ and using the initial conditions
$(x,v,Q)_{\eta=0}=\left(0,0.334,\sqrt{\frac{2}{3}}\right)$, which falls in the
vicinity of the matter dominated fixed point
$\mathcal{P}_{m+}\equiv(x,v,Q)=\left(0,\frac{1}{3},\sqrt{\frac{2}{3}}\right)$.

###  III.6 Comparison with $\Lambda$CDM evolution

In this section we present the comparison of the evolution of the
dimensionless Hubble parameter ($h$), deceleration parameter ($q$), effective
equation of state parameter ($w_{\rm eff}$) and the jerk parameter ($j$)
between the spatially flat $\Lambda$CDM cosmology and a HS cosmology given by
the phase trajectory of Fig.1, setting the same initial conditions for both.
To do the numerical analysis it proves to be easier to work with the redshift
$z=\frac{1-a(t)}{a(t)}$. For the spatially flat $\Lambda$CDM model the
Friedmann equation and the Raychaudhuri equation are

$\displaystyle 3H^{2}$ $\displaystyle=$ $\displaystyle\rho_{m}+\Lambda,$ (57a)
$\displaystyle 2\dot{H}+3H^{2}$ $\displaystyle=$ $\displaystyle\Lambda,$ (57b)

from which one can write the evolution of the dimensionless Hubble parameter
and the deceleration parameter as555For the spatially flat $\Lambda$CDM model,
the jerk parameter $j$ is identically unity [38]

$\displaystyle h(z)=\sqrt{\Omega_{m0}(1+z)^{3}+(1-\Omega_{m0})},$ (58a)
$\displaystyle q(z)=\frac{3\Omega_{m0}(1+z)^{3}}{2h^{2}(z)}-1,$ (58b)
$\displaystyle w_{\rm eff}(z)=-\frac{1}{3}(1-2q(z)),$ (58c) $\displaystyle
j(z)=1.$ (58d)

Noting that $dN=-\frac{dz}{1+z}$, one can express the non-compact dynamical
system in Eq.(12) with respect to $z$:

$\displaystyle\frac{d\tilde{x}}{dz}=$
$\displaystyle-\frac{1}{(z+1)}\Big{[}-2\tilde{v}+4+3\tilde{x}-3(1+w)\tilde{\Omega}+4\tilde{K}-\tilde{x}^{2}-\tilde{x}\tilde{v}+\tilde{x}\tilde{K}\Big{]},$
(59a) $\displaystyle\frac{d\tilde{v}}{dz}=$
$\displaystyle-\frac{\tilde{v}}{(z+1)}\left(\Gamma\tilde{x}-2\tilde{v}+2\tilde{K}+4\right),$
(59b) $\displaystyle\frac{d\tilde{K}}{dz}=$
$\displaystyle\frac{2\tilde{K}}{(z+1)}\left(\tilde{v}-\tilde{K}-1\right),$
(59c) $\displaystyle\frac{d\tilde{\Omega}}{dz}=$
$\displaystyle\frac{\tilde{\Omega}}{(z+1)}\left(-1+3w+\tilde{x}+2\tilde{v}-2\tilde{K}\right),$
(59d) $\displaystyle\frac{d\tilde{y}}{dz}=$
$\displaystyle-\frac{1}{(z+1)}\left(\Gamma\tilde{x}\tilde{v}+\tilde{y}(2\tilde{K}-2\tilde{v}-\tilde{x}+4)\right).$
(59e)

For the spatially flat case $\tilde{K}=0$, the constraint equation (11) gives

$\tilde{\Omega}=1+\tilde{x}+\tilde{y}-\tilde{v},$ (60)

and the non-compact dynamical system (59) reduces to

$\displaystyle\frac{d\tilde{x}}{dz}=$
$\displaystyle-\frac{1}{(z+1)}\Big{[}-2\tilde{v}+4+3\tilde{x}-3(1+w)(1+\tilde{x}+\tilde{y}-\tilde{v})-\tilde{x}^{2}-\tilde{x}\tilde{v}\Big{]},$
(61a) $\displaystyle\frac{d\tilde{v}}{dz}=$
$\displaystyle-\frac{\tilde{v}}{(z+1)}\left(\Gamma\tilde{x}-2\tilde{v}+4\right),$
(61b) $\displaystyle\frac{d\tilde{y}}{dz}=$
$\displaystyle-\frac{1}{(z+1)}\left(\Gamma\tilde{x}\tilde{v}+\tilde{y}(-2\tilde{v}-\tilde{x}+4)\right).$
(61c)

To solve the system in Eq.(61), one needs to set an initial condition
$(\tilde{x},\tilde{v},\tilde{y})_{z_{i}}$ at some $z_{i}$. To compare the
evolution of the spatially flat $\Lambda$CDM model and the spatially flat HS
model with $n=1=C_{1}$, one needs to set the same initial condition for the
two models at some time in the past. Following the idea adopted in Ref.[35],
we set the initial condition at $z=20$. In the remainder of this section we
explain how to set the same initial condition for both the cosmologies at
$z=20$.

Taking the parameter choice $(\Omega_{m0},\Omega_{\Lambda
0})=(0.315,0.685)$666Actual Planck values are $(\Omega_{\Lambda
0},\Omega_{k0},\Omega_{m0})=(0.6843\pm 0.0073,0.0007\pm 0.0019,0.315\pm
0.007)$ [41]. However, since we are comparing between spatially flat
cosmologies, it is a good enough approximation to take
$(\Omega_{m0},\Omega_{\Lambda 0})=(0.315,0.685)$. in Eq.(58), for the
$\Lambda$CDM evolution one gets

$(h,q)_{z=20}=(h_{\Lambda\rm CDM},q_{\Lambda\rm CDM})_{z=20}=(54.0176,0.499).$
(62)

Setting $w=0$ (CDM) and
$\Gamma=\frac{1}{2}\frac{\tilde{v}\tilde{y}}{(\tilde{v}-\tilde{y})^{2}}$ (for
HS model with $n=1=C_{1}$) and using the initial conditions, one can
numerically solve the system (59) to find the evolution of various dynamical
quantities along the trajectory. Evolution of different cosmographic
parameters and the effective equation of state parameters are then obtained as

$\displaystyle q(z)$ $\displaystyle=$
$\displaystyle-1-\frac{\dot{H}}{H^{2}}=1-\tilde{v}(z),$ (63a) $\displaystyle
w_{\textnormal{eff}}(z)$ $\displaystyle=$
$\displaystyle-\frac{2\dot{H}+3H^{2}}{3H^{2}}=\frac{1}{3}(1-2\tilde{v}(z)),$
(63b) $\displaystyle j(z)$ $\displaystyle=$ $\displaystyle
2q^{2}(z)+q(z)-\frac{dq(z)}{dN}=2q^{2}(z)+q(z)+(1+z)\frac{dq(z)}{dz},$ (63c)
$\displaystyle h(z)$ $\displaystyle=$
$\displaystyle\exp\left[\int_{0}^{z}\frac{1+q(z)}{1+z}dz\right].$ (63d)

One now gets from (63a) that

$\tilde{v}(z=20)=1-q(z=20)=0.501.$ (64)

Since we are setting the same initial conditions for the $\Lambda$CDM and HS
cosmology at $z=20$, we assume that gravity at $z=20$ is described well enough
by GR, so that

$\tilde{x}(z=20)=0.$ (65)

One still needs to determine the initial condition $\tilde{y}(z=20)$. To this
goal, we adopt the following procedure. The past asymptotic of the
$\Lambda$CDM model is the matter dominated phase;
$\Omega_{m}\gg\Omega_{\Lambda}$ for $z\gg 1$. We consider the $f(R)$ cosmology
to be at the matter dominated phase at $z=20$, so that

$Q(z=20)=Q|_{\mathcal{P}_{m+}}=\sqrt{\frac{2}{3}},$ (66)

$\mathcal{P}_{m+}$ being the matter dominated fixed point. Setting $K=0$
(since we are considering spatially flat cosmologies here) and $x(z=20)=0$
(since $\tilde{x}(z=20)=0$), $y(z=20)$ can be found from the constraint
equation (17a) to be

$y(z=20)=\frac{1}{3}.$ (67)

One now gets from Eq.(19) that

$\tilde{y}(z=20)=\frac{y(z=20)}{Q^{2}(z=20)}=\frac{1}{2}.$ (68)

Eqs.(64),(65),(68) provides the initial conditions necessary to solve the
system in eq.(59).

Once $Q(z=20)$ is known, using Eq.(19) one can obtain the numerical values of
the compact dynamical variables at $z=20$ from the numerical values of the
corresponding noncompact dynamical variables at $z=20$:

$\left(x,v,Q\right)_{z=20}=\left(Q^{2}\tilde{x},Q^{2}\tilde{v},Q\right)_{z=20}\left(0,0.334,\sqrt{\frac{2}{3}}\right),$
(69)

In the compact phase space this point lies in the vicinity of the matter
dominated fixed point $\mathcal{P}_{m+}$. The phase trajectory shown in Fig.1
passes through this point. Solving the system in Eq.(59) with the initial
conditions given by Eqs.(64),(65),(68) we find the evolution of various
quantities along this particular trajectory.

In Fig.2, we show the evolution of the dimensionless Hubble parameter,
deceleration parameter and the effective equation of state for the trajectory
of Fig.1 with respect to the redshift $z$, and compare it to the $\Lambda$CDM
evolution.

Figure 2: Evolution with respect to the redshift ($z=\frac{1-a}{a}$) of the
dimensionless Hubble parameter ($h$), deceleration parameter ($q$) and
effective equation of state parameter ($w_{\rm eff}$) along the phase
trajectory of Fig.1 (red curve) and the corresponding evolutions for the
$\Lambda$CDM model (blue curve). The initial conditions for the non-compact
variables, corresponding to those used in Fig.1, are
$(\tilde{x},\tilde{v},\tilde{y})=(0,0.501,0.5)$ and are set at $z=20$
following [35].

In Fig.3 we show the evolution of the jerk parameter with respect to the phase
space time variable $\eta$ defined in Eq.(44) and compare with that of the
$\Lambda$CDM model (i.e. $j=1$).

Figure 3: Evolution with respect to the redshift ($z=\frac{1-a}{a}$) of the
jerk parameter ($j$) along the phase trajectory of Fig.1 (red curve) and the
corresponding evolutions for the $\Lambda$CDM model (blue curve). The initial
conditions for the non-compact variables, corresponding to those used in
Fig.1, are $(\tilde{x},\tilde{v},\tilde{y})=(0,0.501,0.5)$ and are set at
$z=20$ following [35]. The plot reaffirms the fact that the phase trajectory
of Fig.1 passes through the De-Sitter saddle $\mathcal{E}_{+}$ and settles
ultimately on the De-Sitter attractor $\mathcal{G}_{+}$, both of which fall on
the surface $j=1$ (see subsection III.5).

###  III.7 Nonsingular bounce, recollapse and cyclic cosmology

A nonsingular bounce and a recollapse can be specified respectively by the
conditions

$\dot{H}|_{H=0}>0,\qquad\dot{H}|_{H=0}<0.$ (70)

Since $Q=\frac{3H}{D}=0$ is zero for a nonsingular bounce or a recollapse, a
phase trajectory corresponding to a bounce or a recollapse must cross the
$Q=0$ hypersurface. For a nonsingular bounce ($\dot{H}>0$) the crossing is
from $Q<0$ side to $Q>0$ side and for a recollapse ($\dot{H}<0$) the crossing
is from $Q>0$ side to $Q<0$ side. Starting from the definition of the Ricci
scalar, one can express $\dot{H}$ as

$\displaystyle\dot{H}$ $\displaystyle=\frac{D^{2}}{9}\left(v-2Q^{2}-K\right).$
(71)

Using the above expression, one can distinguish the regions on the $Q=0$
hypersurface in which the trajectories cross from left to right (bounce) and
from right to left (recollapse) as well as the region in which the
trajectories don’t cross at all. They are respectively

$\displaystyle v>K\qquad\qquad(\text{Bounce})$ (72) $\displaystyle
v<K\qquad\qquad(\text{Recollapse})$ $\displaystyle v=K\qquad\qquad(\text{No
crossing}).$

On the $Q=0$ hypersurface the physical viability conditions (52),(53),(54)
become

$\frac{1}{2}v<1-K-x^{2}<v,\qquad 0\leq K+x^{2}\leq 1,\qquad 0\leq v+x^{2}\leq
1.$ (73)

Clearly, for the spatially flat case $K=0$, the first and the third inequality
are incompatible. This implies a breakdown of the condition
$f^{\prime\prime}>0$. One therefore arrives at an important conclusion:

  * •

No physically viable nonsingular bounce or recollapse is possible in spatially
flat FLRW cosmology for the Hu-Sawicki model with $n=1=C_{1}$.

Note that this conclusion holds only for $K=0$. For $K\neq 0$ perfectly viable
bounces or recollapses can occur. Increasing the value of $K$ from $0$ to $1$
allows for more space for the first and third inequality to be simultaneously
satisfied, and consequently increases the possibility of a perfectly viable
bounce (see 2-dimensional slices of the phase space in Fig.4).

(a)

(b)

(c)

Figure 4: 2-D slices of the phase space for various values of $K\neq 0$ and
$x=0$. The sheet of fixed points $\mathcal{S}_{1}$ intersects these 2-D slices
at curves that are represented by the black curves in the plots. The sheet of
fixed points $\mathcal{S}_{2}$ intersects these 2-D slices only at the two
points $(Q,v)=(\pm\sqrt{1-K},0)$, which, coincidentally, also fall on
$\mathcal{S}_{1}$. Increasing the value of $K$ allows more trajectories to
pass through $Q=0$ as either bounces or recollapses. The green shaded region
corresponds to the subset of the three phase space viability constraints
(52),(53),(54), as projected onto the 2-D $Q$-$v$ slices.

In Fig.5, 6 and 7 we show explicitly 2-D parametric plots of bouncing and
recollapsing trajectories taking the initial conditions from the 2-dimensional
slices in Fig. 4.

(a) (b) (c)

Figure 5: Parametric plots corresponding to the initial conditions
$(K(0),x(0),v(0),Q(0))=(0.2,0,0.9,0)$ (coordinates labelled). The
corresponding point in the phase space is marked with a blue dot in Fig.4(a).
The trajectory evolves away from the saddle fixed point $(K,x,v,Q)=(0,0,1,0)$
on $\mathcal{S}_{1}$ in the past towards the saddle fixed point
$(K,x,v,Q)=(0,-1,0,2)$ on $\mathcal{S}_{2+}$. It then undergoes a recollapse
in the past and a bounce at the initial condition point
$(K,x,v,Q)=(0.2,0,0.9,0)$ (since $v>K$ at this point; see Eq.(72)) then
recollapses again and evolves towards the corresponding collapsing saddle
fixed point $(K,x,v,Q)=(0,1,0,-2)$ on $\mathcal{S}_{2-}$ in the future, before
returning to the point of origin on $\mathcal{S}_{1}$. This entire cycle is
repeated over and over again and this cyclic trajectory corresponds to a
cyclic cosmology. We note here that the recollapse points look to be identical
from Fig.5(b) and 5(c), however the $x$ coordinate switches sign through the
bounce between them, as seen in Fig.5(a).

Figure 6: Parametric plots corresponding to the initial conditions
$(K(0),x(0),v(0),Q(0))=(0.6,0,0.86,0)$ (coordinates labelled). The
corresponding point in the phase space is marked with a blue dot in Fig.4(b).
The trajectory evolves asymptotically from the past attractor
$\mathcal{G}_{-}$ to the future attractor $\mathcal{G}_{+}$, undergoing a
bounce at $(K,x,v,Q)=(0.6,0,0.86,0)$ (since $v>K$ at this point; see Eq.(72)).

Figure 7: Parametric plots corresponding to the initial conditions
$(K(0),x(0),v(0),Q(0))=(0.9,0,0.95,0)$ (coordinates labelled). The
corresponding point in the phase space is marked with a blue dot in Fig.4(c).
The trajectory evolves asymptotically from the past attractor
$\mathcal{G}_{-}$ to the future attractor $\mathcal{G}_{+}$, undergoing a
bounce at $(K,x,v,Q)=(0.9,0,0.95,0)$ (since $v>K$ at this point; see Eq.(72)).

Figure 8: Parametric plots corresponding to the initial conditions
$(K(0),x(0),v(0),Q(0))=(0.9,0,0.2,0)$ (coordinates labelled). The
corresponding point is marked with a blue dot in Fig.4(c). The trajectory
evolves asymptotically from the saddle fixed point $(K,x,v,Q)=(0,-1,0,2)$ on
$\mathcal{S}2_{+}$ in the past to the saddle fixed point
$(K,x,v,Q)=(0,1,0,-2)$ on $\mathcal{S}2_{-}$ in the future, undergoing a
recollapse at $(K,x,v,Q)=(0.9,0,0.2,0)$ (since $v<K$ at this point; see
Eq.(72)).

Here we find it worth mentioning that cyclic cosmologies are possible in the
HS model with $n=1=C_{1}$. This is the case, for example, for the evolution
represented by the parametric plots in Fig.5; it periodically undergoes a
distinct cycle. Cyclic universes have been studied in the context of
$f(R)$-gravity [45, 46, 47], however most results produce numerical solutions
or complex analytical functions for the form of $f(R)$ for which such a
solution is possible. As far as we are aware, explicit cyclic solutions have
not been investigated specifically for physically viable $f(R)$ models like
Hu-Sawicki.

###  III.8 $\Lambda$CDM point for the Hu-Sawicki $f(R)$ model

The particular case of the Hu-Sawicki model we are considering ((41)) has
$n=C_{1}=1$ and only one unspecified model parameter $C_{2}$. Following the
discussion of the last paragraph of Sec.II.3, one can use Eq.(28a) to
determine $C_{2}$, which gives a third order equation for $c=C_{2}H_{0}^{2}$.
However, only one of the three solutions is positive. Since we restrict this
study to $C_{2}>0$, the only viable solution is777We give the value here to
such a degree of accuracy since the dynamical system is very sensitive. It is
possible to keep only up to the third decimal place, however this slightly
alters some of the values found in subsequent calculations.

$c=C_{2}H_{0}^{2}=0.11841.$ (74)

One can then calculate values for $f(R_{0})$ and its derivatives (up to second
order since this is the highest order derivative required for calculation of
the dynamical variable values):

$\displaystyle\frac{f(R_{0})}{H_{0}^{2}}=4.76818,$ (75) $\displaystyle
f^{\prime}(R_{0})=0.76996,$
$\displaystyle\frac{f^{\prime\prime}(R_{0})}{(H_{0}^{2})^{-1}}=0.02613.$

Using these values one can calculate

$\displaystyle D_{0}=$
$\displaystyle\sqrt{\left(3H_{0}+\frac{3}{2}\frac{f^{\prime\prime}(R_{0})\dot{R}_{0}}{f^{\prime}(R_{0})}\right)^{2}+\frac{3}{2}\left(\frac{f(R_{0})}{f^{\prime}(R_{0})}+\frac{6k}{a_{0}^{2}}\right)}$
(76) $\displaystyle=$ $\displaystyle 4.17549H_{0},$

and, therefore, the present day values of the compact dynamical variables

$\displaystyle(x_{0},v_{0},y_{0},Q_{0},K_{0},\Omega_{0})$
$\displaystyle=\left(\frac{3}{2}\frac{f^{\prime\prime}(R_{0})\dot{R_{0}}}{f^{\prime}(R_{0})}\frac{1}{D_{0}},\frac{3}{2}\frac{R_{0}}{D_{0}^{2}},\frac{3}{2}\frac{f(R_{0})}{f^{\prime}(R_{0})}\frac{1}{D_{0}^{2}},\frac{3H_{0}}{D_{0}},1-(Q_{0}+x_{0})^{2}-y_{0},1-x_{0}^{2}-v_{0}\right)$
(77) $\displaystyle=(-0.03469,0.78833,0.53207,0.71848,0.00036,0.21047).$

This is the $\Lambda$CDM point for HS model with $n=C_{1}=1$. In Fig. 9 we
show parametric plots taking the initial conditions at this point.

Figure 9: Parametric plots corresponding to the initial conditions
$(K(0),x(0),v(0),Q(0))=(0.00036,-0.03469,0.78833,0.71848)$, which is the
$\Lambda$CDM point for the HS model with $n=C_{1}=1$ (Eq.(77)). The trajectory
evolves asymptotically from the point $(K,x,v,Q)=(0,0,1,0)$ on
$\mathcal{S}_{1+}$ in the past, passing through the $\Lambda$CDM point and the
point $\mathcal{E}_{+}$ and evolving asymptotically to $\mathcal{G}_{+}$ in
the future.

##  IV Conclusion

In this paper we have presented a comprehensive analysis of the phase space of
$\\{n,C_{1}\\}=\\{1,1\\}$ Hu-Sawicki $f(R)$-gravity model including spatial
curvature. We began by extending the compactified dynamical system formalism
for a general $f(R)$, presented in [35], to include non-spatially flat
cosmologies. We then introduced a method to determine a point and surface in
the phase space corresponding to the present day universe and the $\Lambda$CDM
evolution history respectively (referred to as the $\Lambda$CDM point and
surface). Applying these procedures to the Hu-Sawicki model, we identified
some novel features of the phase space; including invariant submanifolds, a
physically viable region and the existence and stability of two 2-dimensional
sheets of fixed points. We presented a trajectory mimicking the $\Lambda$CDM
evolution history and compared the evolution of their cosmographic parameters.
Following this, we investigated the existence of bouncing and recollapsing
cosmologies in the model and arrived at a _no-go_ conclusion stating that this
kind of trajectory is not possible in the spatially flat model for our
specific parameter choice. We also presented the case of a cyclic cosmology
which, to the best of our knowledge, has not been identified before in this
particular model.

A comment here is in order regarding Ref.[48], in which an analytical method
is developed for identifying bouncing solutions in $f(R)$ theory for spatially
flat FLRW cosmology. This study finds that the model used (a similar late-time
model) does permit bounces for the spatially flat case. The method developed
here provides a way to analytically test the _no-go_ conclusion for the Hu-
Sawicki model with $n=1=C_{1}$. Applying the analytic method to this model, we
find that it is possible to find bounces in the spatially flat case, however
this requires that the matter Equation of State parameter $w_{m}$ cross the
phantom barrier ($w_{m}<-1$). On the surface, it appears that bounces are
possible with $-1\leq w_{m}\leq 1$ but upon further investigation this results
in ghost instabilities since $f^{\prime}<0$. Therefore, if one requires that
$w_{m}$ is non-phantom and that $f^{\prime}>0$ throughout the evolution
history, the _no-go_ conclusion holds for the $n=1=C_{1}$ Hu-Sawicki model.

## Acknowledgements

KM thanks the University of Cape Town for financial support. This work is
based on the research supported in part by the National Research Foundation of
South Africa (Grant Number: 123055). S.C. acknowledges the financial
assistance provided by the North-West University, South Africa through the
postdoctoral grant NWU PDRF Fund NW.1G01487 as well as the accommodation and
financial assistance provided kindly by the Department of Mathematics, BITS-
Pilani, Hyderabad Campus. PD acknowledges financial support from the First
Rand Bank (SA).

##  Appendix A Stability of the invariant submanifolds

Suppose that in the 4-dimensional phase space $K$-$x$-$v$-$Q$, a 3-dimensional
invariant submanifold is given by the condition $\mathbb{S}(K,x,v,Q)=0$, i.e.

$\frac{d\mathbb{S}}{d\tau}\bigg{|}_{\mathbb{S}=0}=0.$ (78)

The stability of this submanifold can be determined as follows. Consider the
family of hypersurfaces $\mathbb{S}(K,x,v,Q)=\mathcal{C}$, which contains the
invariant submanifold as $\mathcal{C}=0$. Let us consider an arbitrary
submanifold with $\mathcal{C}=\delta$ which is infinitesimally close to the
invariant submanifold $\mathcal{C}=0$. Keeping in mind the condition for
$\mathbb{S}=0$ being an invariant submanifold, Taylor expanding
$d\mathbb{S}/d\tau$ around it, one gets

$\displaystyle\frac{d\mathbb{S}}{d\tau}\bigg{|}_{\mathbb{S}=\delta}$
$\displaystyle=$
$\displaystyle\frac{d(d\mathbb{S}/d\tau)}{d\mathbb{S}}\bigg{|}_{\mathbb{S}=0}\delta+\frac{1}{2}\frac{d^{2}(d\mathbb{S}/d\tau)}{d\mathbb{S}^{2}}\bigg{|}_{\mathbb{S}=0}\delta^{2}+\frac{1}{6}\frac{d^{3}(d\mathbb{S}/d\tau)}{d\mathbb{S}^{3}}\bigg{|}_{\mathbb{S}=0}\delta^{3}+...$
(79)

The first nonvanishing term in the above expression gives the directionality
(towards or away from) of the flow component normal to the invariant
submanifold $\mathbb{S}=0$ for the phase flow in its vicinity, which
determines its nature of stability. The invariant submanifold is of attracting
(repelling) nature if the phase flow in its vicinity is towards (away from)
it.

In the HS model under consideration there are three invariant submanifolds,
namely

$K=0,\qquad v=0,\qquad z\equiv v+(Q+x)^{2}+K-1=0.$

###  A.1 Stability of $K=0$

Using the dynamical equation of $K$, the Taylor expansion of $dK/d\tau$ around
$K=0$ up to the first nonvanishing term is

$\frac{dK}{d\tau}\bigg{|}_{K=\delta}=\frac{1}{3}\left[v^{2}x\left((Q+x)^{2}-1\right)-\left((Q+x)^{2}-1+v\right)^{2}\left(-4x+(Q+x)(x^{2}+4xQ+3v-1)\right)\right]\delta+\mathcal{O}(\delta^{2})$
(80)

The invariant submanifold is attracting (repelling) when the coefficient of
$\delta$ is negative (positive).

###  A.2 Stability of $v=0$

Using the dynamical equation of $v$, the Taylor expansion of $dv/d\tau$ around
$v=0$ up to the first nonvanishing term is

$\frac{dv}{d\tau}\bigg{|}_{v=\delta}=\begin{cases}\frac{1}{3}z^{2}\left[Q\left(3-4xQ-5x^{2}\right)+x\left(5-4k-x^{2}\right)\right]\delta+\mathcal{O}(\delta^{2})\hskip
40.00006pt(\text{For $z>0$}),\\\
\frac{1}{3}x\delta^{3}+\mathcal{O}(\delta^{4})\hskip 210.00032pt\,(\text{For
$z=0$}).\end{cases}$ (81)

The leading nonvanishing term in the Taylor expansion is $\mathcal{O}(\delta)$
for $z>0$ and $\mathcal{O}(\delta^{3})$ for $z=0$. The invariant submanifold
is attracting (repelling) when the leading coefficient is negative (positive).

###  A.3 Stability of $z=0$

Using Eq.(46), the Taylor expansion of $dz/d\tau$ around $z=0$ up to the first
nonvanishing term is

$\frac{dz}{d\tau}\bigg{|}_{z=\delta}=\begin{cases}-\frac{1}{3}v^{3}x\delta+\mathcal{O}(\delta^{2})\hskip
50.00008pt(\text{For $x\neq 0,\,v>0$}),\\\
-\frac{2}{3}Qv\delta^{2}+\mathcal{O}(\delta^{3})\hskip 40.00006pt\;\,\
(\text{For $x=0,\,v>0$}),\\\ -(2x-6Q)\delta^{3}+\mathcal{O}(\delta^{4})\hskip
20.00003pt\ (\text{For $v=0$}).\end{cases}$ (82)

The leading nonvanishing term in the Taylor expansion is $\mathcal{O}(\delta)$
for $x\neq 0,v\neq 0$, $\mathcal{O}(\delta^{2})$ for $x=0,v\neq 0$ and
$\mathcal{O}(\delta^{3})$ for $x\neq 0,v=0$.

  * •

For $x\neq 0,\,v>0$, the invariant submanifold is attracting (repelling) for
$x>0$ ($x<0$), with $x=0$ acting as a separatrix.

  * •

For $x=0,v>0$, the invariant submanifold $z=0$ has different nature on both
sides of it, with $z=0$ serving as a node. Since $C_{2}>0$ in the HS model
under consideration, one recalls from Eq.(52) that the physical viability
condition $f^{\prime\prime}>0$ requires $z>0$. If one confines the attention
to the positive side of $z$, then, for $x=0,\,v>0$, the invariant submanifold
$z=0$ is attracting (repelling) for $Q>0$ ($Q<0$).

  * •

For $v=0$, the invariant submanifold is attracting (repelling) for $x>3Q$
($x<3Q$), with $x=3Q$ acting as a separatrix.

One can note that the nature of stability of the invariant submanifold $z=0$
is independent of $K$.

##  Appendix B Stability of the sheets of fixed points
$\mathcal{S}_{1},\,\mathcal{S}_{2}$

From Eqs.(47),(48) we notice that both the sheets of fixed points
$\mathcal{S}_{1},\,\mathcal{S}_{2}$ lie on the invariant submanifold $z=0$. In
fact, if we replace $Q$ by $z$ using the definition of $z$

$Q=-x\pm\sqrt{1-v-K+z},$ (83)

and use $(K,x,v,z)$ instead of $(K,x,v,Q)$ to be the set of independent
dynamical variables, then the equations for the sheets of fixed points are
greatly simplified

$\displaystyle\mathcal{S}_{1}\equiv(K,x,v,z)=\left(K,0,v,0\right),$ (84)
$\displaystyle\mathcal{S}_{2}\equiv(K,x,v,z)=\left(K,x,0,0\right).$ (85)

In general, there can be phase flows within the invariant submanifold $z=0$
(which represents the limit $f(R)\rightarrow-2\Lambda+R$), but the flow is
contained within it. There are, however, parts of $z=0$, given by
$\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, where there is no flow, as these
parts represents sheets of fixed points. In the immediate vicinity of a point
$\mathcal{P}\subset\mathcal{S}_{1}\,\text{or}\,\mathcal{S}_{2}$, the phase
flow is essentially 2-dimensional; If $\mathcal{P}\subset\mathcal{S}_{1}$, the
flow is entirely on the $x$-$z$ plane and if
$\mathcal{P}\subset\mathcal{S}_{2}$, the flow is entirely on the $v$-$z$
plane. This geometrical insight dictates that to investigate the stability of
$\mathcal{S}_{1},\,\mathcal{S}_{2}$, one is only required to considered the
phase flow in the vicinity of the submanifolds $z=0$, $x=0$, $v=0$, something
that we had already considered in Appendix A.

###  B.1 Stability of $\mathcal{S}_{1}$

The 2-dimensional sheet $\mathcal{S}_{1}$ lies at the intersection of the
3-dimensional submanifold $x=0$ and the 3-dimensional invariant submanifold
$z=0$. In the immediate vicinity of any point on $\mathcal{S}_{1}$, the phase
flow is essentially planar, lying on the $x$-$z$ plane. The nature of
stability of the invariant submanifold $z=0$ for $x=0$ has already been
discussed in the last subsection of Appendix A. In particular, for $v>0$, the
flow is towards (away from) $z=0$ for $Q>0$ ($Q<0$). Phase flow in the
vicinity of the submanifold $x=0$ for $z=0$ can be found by Taylor expanding
$dx/d\eta$ around $x=0$, setting $z=0$ and looking at the leading order
nonvanishing term

$\frac{dx}{d\eta}\bigg{|}_{x=\delta}=-\frac{1}{6}v^{3}\delta^{2}+\mathcal{O}(\delta^{3})\qquad(\text{For
{\hbox{z=0,\,v>0}}}).$ (86)

The flow is towards $x=0$ on the positive side of $x$ and away from $x=0$ on
the negative side of $x$, with $x=0$ acting as a node, an inference that holds
true for all values of $Q$. Summarizing, one can infer that the sheet of fixed
points $\mathcal{S}_{1}$ is always

  * •

attracting on the positive side of $x$.

  * •

repelling on the negative side of $x$.

Based on this fact alone, one can conclude that all the points on
$\mathcal{S}_{1}$ are saddle fixed points, with there being always a flow
component along the $x$-direction from the positive to the negative side of
$x$ past $x=0$ in the vicinity of any point
$\mathcal{P}\subset\mathcal{S}_{1}$. One can note that the stability of
$\mathcal{S}_{1}$ is completely independent of $K$ and $v$, i.e. the same
nature of stability holds everywhere on $\mathcal{S}_{1}$.

###  B.2 Stability of $\mathbb{S}_{2}$

The 2-dimensional sheet $\mathcal{S}_{2}$ lies at the intersection of two
3-dimensional invariant submanifolds $v=0$ and $z=0$. In the immediate
vicinity of any point on $\mathcal{S}_{2}$, the phase flow is essentially
planar, lying on the $v$-$z$ plane. The stability of both $v=0$ and $z=0$ have
already been discussed in Appendix A. For $z=0$, the invariant submanifold
$v=0$ is attracting (repelling) for $x<0$ ($x>0$). For $v=0$, the invariant
submanifold $z=0$ is attracting (repelling) for $x>3Q$ ($x<3Q$). Summarizing,
one can infer that the sheet of fixed points $\mathcal{S}_{2}$ is

  * •

attracting for $3Q<x<0$.

  * •

repelling for $0<x<3Q$.

  * •

saddle otherwise.

Given any point on $\mathcal{S}_{2}$, one can now infer about it’s stability.
Note that the stability of $\mathcal{S}_{2}$ depends on $x$ as well as $K$ (by
virtue of Eq.(83)).

## References

## References

  * [1] J. P. Ostriker and Paul J. Steinhardt.  Cosmic concordance.  5 1995.
  * [2] S. Perlmutter et al.  Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae.  Astrophys. J., 517:565–586, 1999.
  * [3] Adam G. Riess et al.  Observational evidence from supernovae for an accelerating universe and a cosmological constant.  Astron. J., 116:1009–1038, 1998.
  * [4] John L. Tonry et al.  Cosmological results from high-z supernovae.  Astrophys. J., 594:1–24, 2003.
  * [5] Robert A. Knop et al.  New constraints on Omega(M), Omega(lambda), and w from an independent set of eleven high-redshift supernovae observed with HST.  Astrophys. J., 598:102, 2003.
  * [6] Adam G. Riess et al.  Type Ia supernova discoveries at z $>$ 1 from the Hubble Space Telescope: Evidence for past deceleration and constraints on dark energy evolution.  Astrophys. J., 607:665–687, 2004.
  * [7] P. Astier et al.  The Supernova Legacy Survey: Measurement of $\Omega_{M}$, $\Omega_{\Lambda}$ and ${\cal w}$ from the first year data set.  Astron. Astrophys., 447:31–48, 2006.
  * [8] D. N. Spergel et al.  First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Determination of cosmological parameters.  Astrophys. J. Suppl., 148:175–194, 2003.
  * [9] D. N. Spergel et al.  Wilkinson Microwave Anisotropy Probe (WMAP) three year results: implications for cosmology.  Astrophys. J. Suppl., 170:377, 2007.
  * [10] Max Tegmark et al.  Cosmological parameters from SDSS and WMAP.  Phys. Rev. D, 69:103501, 2004.
  * [11] Uros Seljak et al.  Cosmological parameter analysis including SDSS Ly-alpha forest and galaxy bias: Constraints on the primordial spectrum of fluctuations, neutrino mass, and dark energy.  Phys. Rev. D, 71:103515, 2005.
  * [12] Shaun Cole et al.  The 2dF Galaxy Redshift Survey: Power-spectrum analysis of the final dataset and cosmological implications.  Mon. Not. Roy. Astron. Soc., 362:505–534, 2005.
  * [13] Daniel J. Eisenstein et al.  Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies.  Astrophys. J., 633:560–574, 2005.
  * [14] Chris Blake, David Parkinson, Bruce Bassett, Karl Glazebrook, Martin Kunz, and Robert C. Nichol.  Universal fitting formulae for baryon oscillation surveys.  Mon. Not. Roy. Astron. Soc., 365:255–264, 2006.
  * [15] Bhuvnesh Jain and Andy Taylor.  Cross-correlation tomography: measuring dark energy evolution with weak lensing.  Phys. Rev. Lett., 91:141302, 2003.
  * [16] Thomas P. Sotiriou and Valerio Faraoni.  f(R) Theories Of Gravity.  Rev. Mod. Phys., 82:451–497, 2010.
  * [17] Antonio De Felice and Shinji Tsujikawa.  f(R) theories.  Living Rev. Rel., 13:3, 2010.
  * [18] Salvatore Capozziello.  Curvature quintessence.  Int. J. Mod. Phys. D, 11:483–492, 2002.
  * [19] Wayne Hu and Ignacy Sawicki.  Models of f(R) Cosmic Acceleration that Evade Solar-System Tests.  Phys. Rev. D, 76:064004, 2007.
  * [20] M. Abdelwahab, S Carloni, and P K. S. Dunsby.  Cosmological dynamics of exponential gravity.  Class. Quant. Grav., 25:135002, 2008.
  * [21] Sante Carloni, Peter K. S. Dunsby, Salvatore Capozziello, and Antonio Troisi.  Cosmological dynamics of R**n gravity.  Class. Quant. Grav., 22:4839–4868, 2005.
  * [22] 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–1776, 2009.
  * [23] Sante Carloni.  A new approach to the analysis of the phase space of $f(R)$-gravity.  JCAP, 09:013, 2015.
  * [24] Mohamed Abdelwahab, Rituparno Goswami, and Peter K. S. Dunsby.  Cosmological dynamics of fourth order gravity: A compact view.  Phys. Rev. D, 85:083511, 2012.
  * [25] Naureen Goheer, Jannie A. Leach, and Peter K. S. Dunsby.  Dynamical systems analysis of anisotropic cosmologies in R**n-gravity.  Class. Quant. Grav., 24:5689–5708, 2007.
  * [26] Naureen Goheer, Jannie A. Leach, and Peter K. S. Dunsby.  Compactifying the state space for alternative theories of gravity.  Class. Quant. Grav., 25:035013, 2008.
  * [27] Jannie A. Leach, Sante Carloni, and Peter K. S. Dunsby.  Shear dynamics in Bianchi I cosmologies with R**n-gravity.  Class. Quant. Grav., 23:4915–4937, 2006.
  * [28] Deon M. Solomons, Peter Dunsby, and George Ellis.  Bounce behaviour in Kantowski-Sachs and Bianchi cosmologies.  Class. Quant. Grav., 23:6585–6597, 2006.
  * [29] Kishore N. Ananda, Sante Carloni, and Peter K. S. Dunsby.  A characteristic signature of fourth order gravity.  Springer Proc. Phys., 137:165–172, 2011.
  * [30] Kishore N. Ananda, Sante Carloni, and Peter K. S. Dunsby.  A detailed analysis of structure growth in $f(R)$ theories of gravity.  Class. Quant. Grav., 26:235018, 2009.
  * [31] S. Carloni, P. K. S. Dunsby, and A. Troisi.  The Evolution of density perturbations in f(R) gravity.  Phys. Rev. D, 77:024024, 2008.
  * [32] Amare Abebe, Alvaro de la Cruz-Dombriz, and Peter K. S. Dunsby.  Large Scale Structure Constraints for a Class of f(R) Theories of Gravity.  Phys. Rev. D, 88:044050, 2013.
  * [33] S. Capozziello, V.F. Cardone, and V. Salzano.  Cosmography of f(R) gravity.  Phys. Rev. D, 78:063504, 2008.
  * [34] A. D. Dolgov and Masahiro Kawasaki.  Can modified gravity explain accelerated cosmic expansion?  Phys. Lett. B, 573:1–4, 2003.
  * [35] Sulona Kandhai and Peter K. S. Dunsby.  Cosmological dynamics of viable f(R) theories of gravity.  10 2015.
  * [36] Saikat Chakraborty, Peter K. S. Dunsby, and Kelly Macdevette.  A note on the dynamical system formulations in $f(R)$ gravity.  In Geometric Foundations of Gravity 2021, 12 2021.
  * [37] Matt Visser.  Jerk and the cosmological equation of state.  Class. Quant. Grav., 21:2603–2616, 2004.
  * [38] Maciej Dunajski and Gary Gibbons.  Cosmic Jerk, Snap and Beyond.  Class. Quant. Grav., 25:235012, 2008.
  * [39] Michel Chevallier and David Polarski.  Accelerating universes with scaling dark matter.  International Journal of Modern Physics D, 10(02):213–223, Apr 2001.
  * [40] Eric V. Linder.  Exploring the expansion history of the universe.  Phys. Rev. Lett., 90:091301, Mar 2003.
  * [41] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, and et al.  Planck 2018 results.  Astronomy & Astrophysics, 641:A6, Sep 2020.
  * [42] Peter K. S. Dunsby, Emilo Elizalde, Rituparno Goswami, Sergei Odintsov, and Diego Saez Gomez.  On the LCDM Universe in f(R) gravity.  Phys. Rev. D, 82:023519, 2010.
  * [43] Saikat Chakraborty, Kelly MacDevette, and Peter Dunsby.  A model independent approach to the study of $f(R)$ cosmologies with expansion histories close to $\Lambda$CDM.  Phys. Rev. D, 103(12):124040, 2021.
  * [44] Santos, B., Campista, M., Santos, J., and Alcaniz, J. S.  Cosmology with hu-sawicki gravity in the palatini formalism.  A&A, 548:A31, 2012.
  * [45] Petar Pavlović and Marko Sossich.  Dynamic properties of cyclic cosmologies.  Physical Review D, 103(2), January 2021.
  * [46] Nasr Ahmed and Sultan Z. Alamri.  A cyclic universe with varying cosmological constant in f(r, t) gravity.  Canadian Journal of Physics, 2019.
  * [47] Shin'ichi Nojiri, S. D. Odintsov, and D. Saez-Gomez.  Cyclic, ekpyrotic and little rip universe in modified gravity.  In AIP Conference Proceedings. AIP, 2012.
  * [48] Ali R. Amani.  The bouncing cosmology with $F(R)$ gravity and its reconstructing.  Int. J. Mod. Phys. D, 25(06):1650071, 2016.


$\mathbb{R}_{+}[-(-i)^{m+r+1}e^{i\hat{\theta}_{top}}]$ (which is where the
path $Z_{X}(t)$ starts), it is evident by construction that after wrapping
around ${\rm ind}_{q}^{\infty}(R_{b})$ times we have

(6.14) ${\rm arg}(Z_{X}([\alpha]))<{\rm arg}(Z_{D_{\infty}}([\alpha]))<{\rm
arg}(Z_{X}([\alpha]))+\pi.$

We remark that this inequality only holds because our choice of how to lift
${\rm arg}(Z_{X}([\alpha]))$ depended on the region containing $z_{2}$.

Now we turn to $Z_{P}([\alpha])$, and see how its argument relates to
existence of a solution to (6.2). Note in the generic case there are exactly
${\rm ind}_{0}^{\infty}(R_{0})$ curves from $\mathcal{D}_{0}$ crossing the
positive $y$-axis. For a solution to exists, we need to find a level curve of
$\mathcal{C}_{0}$ connecting $z_{2}$ to the origin. As we have discussed
$z_{2}$ lies in $A_{m+r+1-{\rm ind}_{q}^{\infty}(R_{b})}$, thus if ${\rm
ind}_{q}^{\infty}(R_{b})<{\rm ind}_{0}^{\infty}(R_{0})$ any path connecting
$z_{2}$ to the origin must cross $\mathcal{D}_{0}$ in ${\mathcal{H}}_{x}$ and
no solution exists. In the generic case there are $r$ curves from
$\mathcal{D}_{0}$ emanating from the origin. Thus the condition

${\rm ind}_{0}^{\infty}(R_{0})\leq{\rm ind}_{q}^{\infty}(R_{b})\leq{\rm
ind}_{0}^{\infty}(R_{0})+r$

is equivalent to existence of a solution to (6.2) by Theorem 1 $(ii)$. Now,
${\rm ind}_{0}^{\infty}(R_{0})$ counts the number of times the path $Z_{P}(t)$
crosses the line determined by $Z_{X}([\alpha])$, while ${\rm
ind}_{q}^{\infty}(R_{b})$ counts how many times the path $Z_{X}(t)$ starting
at $-(-i)^{m+r+1}e^{i\hat{\theta}_{top}}$ travels distance $\pi$. Since ${\rm
ind}_{0}^{\infty}(R_{0})\leq{\rm ind}_{q}^{\infty}(R_{b})$ we must always have
${\rm arg}(Z_{X}([\alpha]))<{\rm arg}(Z_{P}([\alpha]))$ for a solution to
exist. Additionally, since ${\rm ind}_{0}^{\infty}(R_{0})\geq{\rm
ind}_{q}^{\infty}(R_{b})-r,$ the most $Z_{X}(t)$ can continue to wind
clockwise (lowering the angle) after the first time it passes
$Z_{P}([\alpha])$ is $(r+1)\pi$. Thus a solution to (6.2) exists if and only
if

(6.15) ${\rm arg}(Z_{X}([\alpha]))<{\rm arg}(Z_{P}([\alpha]))<{\rm
arg}(Z_{X}([\alpha]))+(r+1)\pi.$

In the non-generic case, there are $r+1$ curves emanating from the origin,
instead of $r$. In this case, Theorem 1 $(iii)$ implies that a smooth solution
to (6.2) exists if and only if

${\rm ind}_{0}^{\infty}(R_{0})+\frac{1}{2}\leq{\rm
ind}_{q}^{\infty}(R_{b})\leq{\rm ind}_{0}^{\infty}(R_{0})+r-\frac{1}{2}.$

Above we noted that in the non-generic case, the Cauchy index ${\rm
ind}_{\epsilon}^{\infty}(R_{0})={\rm ind}_{0}^{\infty}(R_{0})-\frac{1}{2}$
counts the number of times $Z_{P}(t)$ crosses the line determined by
$Z_{X}([\alpha])$. Thus we see $Z_{X}(t)$ must wind clockwise past
$Z_{P}([\alpha])$ and then travel a further distance $\pi$, implying ${\rm
arg}(Z_{X}([\alpha]))+\pi<{\rm arg}(Z_{P}([\alpha]))$. While as above,
$Z_{X}(t)$ can continue to wind clockwise (lowering the angle) after the first
time it passes $Z_{P}([\alpha])$, of a distance $(r+1)\pi$. By these two
observations we see a smooth solution to (6.2) exists if and only if

(6.16) ${\rm arg}(Z_{X}([\alpha]))+\pi<{\rm arg}(Z_{P}([\alpha]))<{\rm
arg}(Z_{X}([\alpha]))+(r+1)\pi.$

Now, dividing equations (6.14), (6.15), and (6.16) by $\pi$, and using the
definition of the grade of $V$ from the introduction

$\phi_{V}([\alpha])=\frac{1}{\pi}{\rm arg}(Z_{V}([\alpha])),$

we are able to conclude the proof of Theorem 3.

## References

  * [1] M. Atiyah and R. Bott, The Yang-Mills equations over Riemann surfaces. Phil. Trans. Roy. Soc. London A 308 (1982), 523-615.
  * [2] A. Benesyei, Peano’s Unnoticed Proof of Borel’s Theorem. Amer. Math. Monthly 121 (2014), no. 1, 69-72.
  * [3] W.M. Boothby, The topology of regular curve families with multiple saddle points. Amer. J. Math. 73 (1951), 405-438.
  * [4] T. Bridgeland Stability conditions on triangulated categories, Ann. of Math. (2) 166 (2007), no. 2, 317-345.
  * [5] E. Calabi, Extremal Kähler metrics, Seminar on Differential Geometry, Vol. 102 of Ann. Math. Studies, Princeton Univ. Press, Princeton, N.J. (1982), 259-290.
  * [6] E. Calabi, Métriques Kählériennes et fibrés holomorphes, Annales scientifiques de l’fÉ.N.S. $4^{e}$ série, tome 12, n. 2, (1979), p. 269-294.
  * [7] G. Chen, The J-equation and the supercritical deformed Hermitian-Yang-Mills equation. Invent. Math. 225 (2021), no. 2, 529-602.
  * [8] X. Chen, S. Donaldson and S. Sun, Kähler-Einstein metrics on Fano manifolds, I: approximation of metrics with cone singularities. JAMS 28 (2015), 183-197.
  * [9] X. Chen, S. Donaldson and S. Sun, Kähler-Einstein metrics on Fano manifolds, II: limits with cone angle less than $2\pi$. JAMS 28 (2015), 199-234.
  * [10] X. Chen, S. Donaldson and S. Sun, Kähler-Einstein metrics on Fano manifolds, III: limits as cone angle approaches $2\pi$ and completion of the main proof. JAMS 28 (2015), 235-278.
  * [11] T.C. Collins, A. Jacob, and S.-T. Yau, (1,1) forms with specified Lagrangian phase. Camb. J. Math. 8 (2020), no. 2, 407-452.
  * [12] T.C. Collins and G. Székelyhidi, Convergence of the J-flow on toric manifolds J. Differential Geom. 107 (2017), no. 1, 47-81.
  * [13] T.C. Collins and S.-T. Yau, Moment maps, nonlinear PDE, and stability in mirror symmetry, arXiv:1811.04824 [math.DG].
  * [14] T.C. Collins and S.-T. Yau, Moment maps, nonlinear PDE, and stability in mirror symmetry I; Geodesics, Ann. PDE 7 (2021), no. 1, Paper No. 11, 73 pp.
  * [15] J. Chu, M.-C. Lee, and R. Takahashi A Nakai-Moishezon type criterion for supercritical deformed Hermitian-Yang-Mills equation, arXiv:2105.10725.
  * [16] O. Debarre, Higher-dimensional algebraic geometry, Universitext. Springer-Verlag, New York, 2001. xiv+233 pp.
  * [17] Simon Donaldson, A new proof of a theorem of Narasimhan and Seshadri. Journal of Differential Geometry 18 (1983), 269-277.
  * [18] Simon Donaldson, Anti-self-dual Yang-Mills connections on complex algebraic surfaces and stable vector bundles. Proc. London Math. Soc. 50 (1985), 1-26.
  * [19] M. Douglas, B. Fiol, and C. Römelsberger, Stability and BPS branes, J. High Energy Phys. 2005, no. 9, 006, 15 pp.
  * [20] M. Eisermann, The fundamental theorem of algebra made effective: an elementary real-algebraic proof via Sturm chains. Amer. Math. Monthly 119 (2012), no. 9, 715-752.
  * [21] H. Fang and M. Lai, Convergence of general inverse $\sigma_{k}$-flow on Kähler manifolds with Calabi ansatz, Trans. Amer. Math. Soc. 365 (2013), no. 12, 6543-6567.
  * [22] A. Jacob and N. Sheu The deformed Hermitian-Yang-Mills equation on the blowup of $\mathbb{P}^{n}$. arXiv:2009.00651
  * [23] A. Jacob and S.-T. Yau, A special Lagrangian type equation for holomorphic line bundles, Math. Ann. 369 (2017). no 1-2, 869-898.
  * [24] R.P. Jerrard and L.A. Rubel, On the curvature of the level lines of a harmonic function. Proc. Amer. Math. Soc. 14 (1963), 29-32.
  * [25] M. Lejmi, and G. Székelyhidi, The $J$-flow and stability, Advances in Math. 274 (2015), 404-431.
  * [26] C. Leung, S.-T. Yau, and E. Zaslow, From special Lagrangian to Hermitian-Yang-Mills via Fourier-Mukai transform, Winter School on Mirror Symmetry, Vector Bundles and Lagrangian Submanifolds (Cambridge, MA, 1999), 209-225, AMS/IP Stud. Adv. Math., 23, Amer. Math. Soc., Providence, RI, 2001.
  * [27] M. Marino, R. Minasian, G. Moore, and A. Strominger, Nonlinear Instantons from Supersymmetric p-Branes, hep-th/9911206.
  * [28] N. Sheu, The deformed Hermitian-Yang-Mills equation with Calabi Ansatz, Ph.D. Thesis. UC Davis.
  * [29] J. P. Solomon The Calabi homomorphism, Lagrangian paths and special Lagrangians, Math. Ann 357 (2013), no. 4, 1389-1424.
  * [30] J. Song, Ricci flow and birational surgery, arXiv:1304.2607.
  * [31] J. Song and B. Weinkove, Contracting exceptional divisors by the Kähler-Ricci flow, Duke Math. J. 162 (2013), no. 2, 367-415.
  * [32] J. Song and B. Weinkove, Contracting exceptional divisors by the Kähler-Ricci flow II, Proc. Lond. Math. Soc. (3) 108 (2014), no. 6, 1529-1561.
  * [33] J. Song and B. Weinkove, The Kähler-Ricci flow on Hirzebruch surfaces, J. Reine Angew. Math. 659 (2011), 141-168.
  * [34] J. Song and Y. Yuan, Metric flips with Calabi ansatz, Geom. Func. Anal. 22 (2012), no. 1, 240-265.
  * [35] R. P. Thomas Moment maps, monodromy, and mirror manifolds, Symplectic geometry and mirror symmetry (Seoul, 2000), 467-498, World Sci. Publ., River Edge, NJ, 2001\.
  * [36] K. Uhlenbeck and S.-T. Yau, On the existence of Hermitian Yang-Mills connections in stable vector bundles. Communications on Pure and Applied Mathematics 39 (1986), 257-293.
  * [37] S.-T. Yau, Calabi’s conjecture and some new results in algebraic geometry. Proceedings of the National Academy of Sciences of the United States of America 74 (1977),1798-1799.