# Membrane Thickness Sensitivity
of Avian Prestin: Implications

Kuni H Iwasa

(National Institute for Deafness and Other Communication Disorders,
National Institutes of Health, Bethesda, MD)

###### Abstract

Avian prestin is sensitive to membrane thickness as much as mammalian prestin,
which undergoes conformational transitions in membrane area and thereby drives
length changes of the cylindrical cell body of outer hair cells. The membrane
thickness dependence of mammalian prestin stems from changes in hydrophobic
profile in conformational states, accompanied by changes in their membrane
area. Even though such area changes are not detected for avian prestin, it
nonetheless bends hair bundles of avian short hair cells. Here it is suggested
that the motile function of avian prestin can be based on conformational
transitions involving shearing deformation of the membrane protein, which also
leads to membrane thickness sensitivity.

## Introduction

Membrane thickness dependence of a membrane protein arises from difference in
hydrophobic profile of its conformational states. Specifically for mammalian
prestin, which undergoes changes in the surface area [1] coupled with charge
transfer $q$ across the membrane, membrane thickness dependence was expected
since the volume of the protein is conserved during conformational changes.

Indeed, that is the case with mammalian prestin: A reduction of membrane
thickness shifts the transition voltage of the protein in the positive
direction and an increase in membrane thickness has an opposite effect [2].
This observation is consistent with the expectation that a decreased
hydrophobic thickness of the protein is associated with an increase in the
surface area on hyperpolarization [3, 4].

Membrane thickness dependence can be closely associated with surface area
changes during conformational transitions, which result in changes in the
hydrophobic profile associated with thickness changes. It can be also
associated with changes in the contour length of the hydrophobic interface
between the protein and membrane lipid.

The mode of motion, with which avian prestin is associated, differs from that
of mammalian prestin. Avian prestin is associated with bending of hair bundles
of avian short hair cells [5], which is quite different from length changes of
cell body, which mammalian prestin drives [6].

How can these observations provide a physical picture? The present report is
an attempt to address this issue.?

## Membrane thickness dependence

Chicken prestin shows considerable membrane thickness dependence accompanied
by somewhat smaller charge transfer compared with mammalian prestin [7].

| $V_{\mathrm{pk}}$ | q | thickness sensitivity
---|---|---|---
| (mV) | ($e$) | (mV/%)
gerbil | $-88\pm 11$ | $0.73\pm 0.07$ | $2.7\pm 0.2$
platypus | $-56\pm 11$ | $0.79\pm 0.10$ | $4.8\pm 0.1$
chicken | $54\pm 11$ | $0.35\pm 0.12$ | $8.9\pm 0.7$
Table 1: Here, Peak potential $V_{\mathrm{pk}}$, motor charge $q$, and
thickness sensitivity. Membrane thickness sensitivity is expressed as voltage
shifts per % change in the linear capacitance. $e$: the electronic charge.
Taken from Ref. [7].

Avian prestin has motile charge $q$ smaller than that of mammalian prestin.
However, its membrane thickness dependence is larger than that of mammalian
counterpart. Another significant difference in the operating voltage $V_{pk}$
can be attributed to rather depolarized membrane potential of avian hair cells
[8], which is associated large hair bundle current due to cellular turning of
the avian ear [9, 10].

### Interpretation

Membrane thickness sensitivity results from difference in interaction energy
of a membrane protein, which undergoes conformational transitions, with
membrane lipid.

Let us assume a model membrane protein has two states $S_{0}$ and $S_{1}$.
Assume that $S_{0}$ and $C_{1}$ has, respectively, length of circumference
$L_{0}$ and $L_{1}=L_{0}+\Delta L$, at which it interacts with membrane lipid.
Let $\mathcal{E}(d-d_{0})$ be the energy per unit length due to hydrophobic
mismatch in state $S_{0}$ with membrane lipid with thickness $d$, where
$d_{0}$ is the characteristic hydrophobic thickness of this state. Let $d_{1}$
be the characteristic thickness of $S_{1}$.

The energy $E_{0}$ and $E_{1}$ due to hydrophobic mismatch respectively for
$S_{0}$ and for $S_{0}$ can be expressed by

$\displaystyle E_{0}(d)$ $\displaystyle=L_{0}\mathcal{E}(d-d_{0}),$ (1a)
$\displaystyle E_{1}(d)$ $\displaystyle=(L_{0}+\Delta L)\mathcal{E}(d-d_{1}).$
(1b)

It is likely that the function $\mathcal{E}(x)$ is dominated by even powers of
$x$ because it should increase with the absolute value of the variable $x$.
The simplest case could be $\mathcal{E}(x)=ax^{2}$, where $a$ is a constant,
but it likely includes higher-order terms. The difference in energy between
the two states due to hydrophobic incompatibility is
$E_{\mathrm{diff}}(d)=E_{1}(d)-E_{0}(d)$.

Now the effect of membrane thickness change $d\rightarrow d+\Delta d$. This
change shifts the energy difference between the two states by

$\displaystyle\Delta E=E_{\mathrm{diff}}(d+\Delta d)-E_{\mathrm{diff}}(d).$
(2)

This energy difference changes the relative weight of the conformation states,
i.e. the conformational state of lower energy is favored. That leads to a
shift of voltage dependence of prestin.

Figure 1: Two types of conformational transitions, which result in membrane
thickness dependence. Shear transitions (left) undergoe without membrane area
changes. A transition from a circle to an ellipse, for example, changes the
circumference length from $L_{0}$ to $L_{0}+\Delta L$ without changing the
area. Area transitions (right) undergo change in the membrane area, which
result in thickness change, affecting hydrophobic mismatch.

If we can assume that $\Delta d/d\ll 1$ but $\Delta L/L$ may not be small, Eq.
2 can be approximated by the first-order terms of expansion, i.e.

$\displaystyle\Delta E/\Delta
d=L_{0}(\mathcal{E}^{\prime}(d-d_{1})-\mathcal{E}^{\prime}(d-d_{0}))+\Delta
L\mathcal{E}^{\prime}(d-d_{1})$ (3)

where $\mathcal{E}(x)^{\prime}$ is the first derivative of $\mathcal{E}(x)$
with respect to $x$. In the simplest case $\mathcal{E}(x)=ax^{2}$, the first
term is proportional to $d_{1}-d_{0}$. This special case illustrates that the
first term in Eq. 3 is due to the difference in hydrophobic thickness of the
two states. The second term is due to the difference $\Delta L$ in the
circumference in the two conformational states.

For mammalian prestin, which includes membrane area changes during
conformational transitions [1], changes in the hydrophobic thickness is
expected by assuming volume conservation. In addition, changes in the
circumference can be expected as schematically illustrated in Fig. 1.

Eq. 3 also shows that membrane thickness dependence, however, can arise from
the difference $\Delta L$ in the circumference alone without changes in the
surface area. If these membrane proteins are randomly oriented, membrane share
should cancel each other and no collective motion is expected. This
expectation is consistent with the observation that plasma membrane of HEK
cells transfected with avian prestin did not show movement of cell membrane
elicited by voltage pulses, unlike mammalian prestin [11]. That is because
those avian prestin molecules expressed at low densities in the plasma
membrane of those transfected cells are likely oriented randomly. The
randomness of orientation of this protein does not lead to macroscopic
movement.

## Macroscopic shear

Mean membrane shear can be achieved by shearing mode of a membrane protein as
illustrated in Fig. 2. The main advantage of this mode compared with area
motor mode is that bending shear can be provided without geometrical asymmetry
of the cell membrane. This advantage, however, requires alignment of the
motile molecules.

Figure 2: A schematic representation of creating membrane shear. Molecular
shear produced by conformational transitions can lead to membrane shear if the
motile molecules are aligned. Two of the possible alignments are shown. The
red and blue outlines illustrate two conformational states as Fig. 1 does.

## Bending of hair bundle

It has been found that avian prestin is also associated with motility [5].
Instead of the cell body as mammalian prestin does, avian prestin drives hair
bundles of short hair cells and thereby moves the tectorial membrane
associated with those hair bundles [5].

For the mechanism of this motion, localization of avian prestin is critically
important. Histochemical study showed that prestin is located in the lateral
plasma membrane [5]. The authors point out in the appendix that avian short
hair cells prominently expand toward the apex on the neural side and the
cuticular place, which is associated with the hair bundle does not cover the
neural side of the apical membrane [5].

They propose that the neural side of the lateral membrane contracts on
depolarization, tilting the hair bundle in the neural direction [5].

## Discussion

The mechanism proposed by Beurg et al. [5] are based on the assumption that
avian prestin undergoes membrane area changes, or “area transitions,” similar
to mammalian prestin. This proposal is of interest even though the assumption
contradicts experimental observations [11].

While the structural asymmetry can be confirmed by prestin-antibody staining,
the asymmetry is not so large. It is unclear if relatively small asymmetry of
the cell morphology as revealed by these images is sufficient to produce
bending of the cell if the mode of the molecular displacement is area changes.

Another issue is that the lateral membrane is in contact with the lateral
membrane of another cell, unlike mammalian outer hair cells. A necessary
condition is that the neighboring cell is compliant in the direction of the
movement. This requirement would apply to the submembranous cytoskeletal
structure of those cells.

The neural side of the apical membrane appears to be the localization, which
is logical for driving the hair bundle. However, antibody-staining does not
show the presence of prestin. An analogous argument can be made with a
molecular motor, which produces shear instead of membrane area changes.

Located in the lateral membrane, molecular “shear generator” could bend the
hair bundle independent of the morphological asymmetry. However, generation of
regional shear of the plasma membrane by molecular shear generators requires
molecular alignment of the shear generators. Such alignment could be based on,
e.g. cortical protein organization, and is plausible in view of the polarized
organization of hair cells.

Those observations described here are of theoretical nature. To clarify the
mechanism of hair bundle bending of avian short hair cells, further studies in
both molecular and cellular levels are essential. Those studies concern
details of localization and orientation of avian prestin on the cellular
level. In addition, experimental determination of the upper bound of membrane
area changes and studies of molecular structure of avian prestin would be
critical.

## References

  * Adachi and Iwasa [1999] Adachi, M., and K. H. Iwasa, 1999.  Electrically driven motor in the outer hair cell: Effect of a mechanical constraint.  _Proc. Natl. Acad. Sci. USA_ 96:7244–7249.
  * Fang et al. [2010] Fang, J., C. Izumi, and K. H. Iwasa, 2010.  Sensitivity of prestin-based membrane motor to membrane thickness.  _Biophys J_ 98:2831–2838.
  * Bavi et al. [2021] Bavi, N., M. D. C. abd Gustavo F Contreras, R. Shen, B. G. Reddy, W. Milewski, and E. Perozo, 2021.  The conformational cycle of prestin underlies outer hair cell electromotility.  _Nature_ 600:553–558.
  * Futamata et al. [2022] Futamata, H., M. Fukuda, R. Umeda, K. Yamashita, A. Tomita, S. Takahashi, T. Shikakura, S. Hayashi, T. Kusakizako, T. Nishizawa, K. Homma, and O. Nureki, 2022.  Cryo-EM structures of thermostabilized prestin provide mechanistic insights underlying outer hair cell electromotility.  _Nat Commun_ 6208\.
  * Beurg et al. [2013] Beurg, M., X. Tan, and R. Fettiplace, 2013.  A prestin motor in chicken auditory hair cells: active force generation in a nonmammalian species.  _Neuron_ 79:69–81.
  * Ashmore [1987] Ashmore, J. F., 1987.  A fast motile response in guinea-pig outer hair cells: the molecular basis of the cochlear amplifier.  _J. Physiol._ 388:323–347.
  * Izumi et al. [2011] Izumi, C., J. E. Bird, and K. H. Iwasa, 2011.  Membrane thickness dependence of prestin orthologs: The evolution of a piezoelectric protein.  _Biophys. J_ 100:2614–2622.
  * Tan et al. [2013] Tan, X., M. Beurg, C. Hackney, S. Mahendrasingam, and R. Fettiplace, 2013.  Electrical tuning and transduction in short hair cells of the chicken auditory papilla.  _J. Neurophysiol._ 109:2007 – 2020.
  * Crawford and Fettiplace [1981] Crawford, A. C., and R. Fettiplace, 1981.  An electrical tuning mechanism in turtle cochlear hair cells.  _J. Physiol._ 312:377–412.
  * Iwasa [2022] Iwasa, K. H., 2022.  Of Mice and Chickens: Revisiting the RC Time Constant Problem.  _Hearing Res._ 423:108422.
  * Tan et al. [2011] Tan, X., J. L. Pecka, J. Tang, O. E. Okoruwa, Q. Zhang, K. W. Beisel, and D. Z. Z. He, 2011.  From zebrafish to mammal: functional evolution of prestin, the motor protein of cochlear outer hair cells.  _J Neurophysiol_ 105:36–44.


# A model independent approach to the study of $f(R)$ cosmologies with
expansion histories close to $\Lambda$CDM

Saikat Chakraborty<EMAIL_ADDRESS><EMAIL_ADDRESS>Center for
Gravitation and Cosmology, College of Physical Science and Technology,
Yangzhou University, Yangzhou 225009, China
and
International Center for Cosmology, Charusat University, Anand 388421, Gujrat,
India  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

###### Abstract

We propose a new framework for studying the cosmology of $f(R)$ gravity which
completely avoids using the reconstruction programme. This allows us to easily
obtain a qualitative feel of how much the $\Lambda$CDM model differs from
other $f(R)$ theories of gravity at the level of linear perturbation theory
for theories that share the same background dynamics. This is achieved by
using the standard model independent cosmographic parameters to develop a new
dynamical system formulation of $f(R)$ gravity which is free from the
limitation of having to first specify the functional form of $f(R)$. By
considering a set of representative trajectories, which are indistinguishable
from $\Lambda$CDM, we use purely qualitative arguments to determine the extent
to which these models deviate from the standard model by including an analysis
of the linear growth rate of density fluctuations and also whether or not they
suffer from the Dolgov-Kawasaki instability. We find that if one demands that
a late time $f(R)$ cosmology is observationally close to the $\Lambda$CDM
model, there is a higher risk that it suffers from a Dolgov-Kawasaki
instability. Conversely, the more one tries to construct a physically viable
late time $f(R)$ cosmology, the more likely it is observationally different
from the $\Lambda$CDM model.

##  I Introduction

Assuming that its geometry is well described by a Friedmann-Robertson-Walker
metric, there is now considerable observational evidence [1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13] supporting the idea that our Universe is currently
undergoing an accelerated phase of expansion. Given that gravity is attractive
in standard General Relativity (GR) in the presence of both relativistic and
non-relativistic matter components, these observations suggest two different
possibilities. The first proposes that there exists a matter component which
currently dominates the energy density of the Universe, in the presence of
which, gravity is repulsive even within standard General Relativity. Because
of the lack of a physical description of this matter component, this term in
the governing field equations is referred to as “Dark Energy”. The simplest
such model adds a _cosmological constant_ $\Lambda$ to the Einstein-Hilbert
Lagrangian, which behaves effectively like a perfect fluid with an equation of
state parameter $w=-1$, whose energy density remains constant with time.
Together with another unknown matter component called Cold Dark Matter (CDM)
this description of the Universe is known as the $\Lambda$CDM model (or
_Concordance Model_ [14]). $\Lambda$CDM appears to be almost perfectly
consistent with current observations. Other models in these category introduce
either an extra perfect fluid or an extra scalar field to produce a late time
accelerating cosmological solution (See Refs. [15, 16] for a review on various
Dark Energy models).

The second possibility is that there is no dark energy, but instead, the
acceleration of the present Universe is accounted for by a modification of how
gravity behaves at very large cosmological scales. This line of thought
differs from the first proposal in that it replaces our ignorance of what the
dominant matter component is by a lack of understanding of how gravity works
at very large distance scales. Consequently, many so called theories of
“modified gravity” have been developed in order to provide a more geometrical
explanation of the late time acceleration of the Universe [17, 18]. Almost all
such models belong to the class of scalar-tensor theories which involve a
modification of the Einstein-Hilbert Lagrangian of the following generic form:

$R\longrightarrow f(R,\phi)\,,$

while the standard matter part of the Lagrangian remains intact. Careful
analysis of the dynamics of such theories show that such theories in general
contain two extra scalar degrees of freedom: one is the scalar field $\phi$
and the other is hidden in the form of $\varphi\equiv F(R,\phi)=\frac{\partial
f}{\partial R}$. The simplest modified gravity models include only one scalar
degree of freedom, which include generalised Brans-Dicke theories ($f(\phi)R$)
and $f(R)$ theories. The first ever modified gravity model for dark energy
presented in [19] was of this type and had the following Lagrangian:

$f(R)=R-\frac{\mu^{4}}{R}\,,$ (1)

although shortly afterwards this particular $f(R)$ theory was shown to be
plagued by an instability [20]. Later attempts in this line gave birth to
several different $f(R)$ gravity models which can describe a late time
acceleration of the universe while being free from instabilities, some of the
well-known models being the Hu-Sawicki model [21], Starobinski model [22], MJW
model [23] etc. In what following we will focus our concern to the $f(R)$
class of modified gravity only.

No matter how many different $f(R)$ gravity models have been proposed as an
alternative to $\Lambda$CDM, the $\Lambda$CDM model still remains the one that
best fits current observations. The problem with the $\Lambda$CDM model is
that the only known candidate for $\Lambda$, namely the quantum vacuum energy
as calculated in quantum field theory, has a value many orders of magnitude
higher than its observed value. This is the principle motivation for looking
for alternative dark energy models. Because $f(R)$ theories present us with an
extra dynamical degree of freedom to play with, it is _in principle_ possible
most of the time to find a class of $f(R)$ theories that can produce the exact
same cosmological solution as one finds in General Relativity, _i.e.,_ the
solution is not unique to GR [24]. The systematic bottom-up approach to
finding the class of $f(R)$ gravity models that produces a desired
cosmological solution is called a reconstruction programme of $f(R)$ gravity
[25]. Even though a number of such reconstruction techniques exist, these
methods are not always helpful. More often than not a compact functional form
of $f(R)$ cannot be found and even if they can be found, they involve
functions that are too complicated for further analytical treatment. For
example, since the $\Lambda$CDM model best fits the observations, there have
been attempts to reconstruct $f(R)$ theories that can exactly mimic the
$\Lambda$CDM evolution history without invoking a $\Lambda$-term, which gave
back compact solutions in the form of Hypergeometric functions [24, 26, 27].
If one tries to analyze further into the perturbative level to look whether
perturbation-dependent observables can distinguish between $R-2\Lambda$ and
$f(R)$ models, Hypergeometric functions are not the easiest functions to deal
with.

The main purpose of this paper is to completely avoid using the reconstruction
programme, while still being able to get a qualitative feel for how much the
$\Lambda$CDM model differs at the linear perturbation level from other $f(R)$
theories of gravity, based on late time models that produce the exact same
background dynamics.

A very good approach for qualitative understanding of cosmological models is
the dynamical system approach, first developed by Collins [28] and extensively
reviewed in the book edited by Ellis and Wainwright [29] (see also [30, 31,
16]). In fact, using these techniques to study cosmological models has the
advantage of providing a relatively simple method for obtaining exact
solutions, which appear as fixed points of the system, and obtaining a global
picture of the dynamics of these models. All the autonomous dynamical system
formulations for $f(R)$ gravity that have appeared in literature up to now
require one to specify the functional form of $f(R)$ in order to be able to
write down a closed system of first order nonlinear ODEs [32, 33, 34, 35, 36].
This approach is therefore not particularly helpful for the general dynamical
analysis required in this present work, where the functional form of $f(R)$ is
not known _a-priori_ (see however [37]). We however circumvent this problem by
proposing a novel dynamical system formulation of $f(R)$ gravity that is free
from this limitation. The trick is to extend the phase space by including a
set of dimensionless cosmographic parameters. Cosmographic parameters are key
cosmological observables based on performing a Taylor expansion of the scale
factor around the present redshift and, as we will show later on with the
specific example of the $\Lambda$CDM case, the choice of a cosmological
solution places a constraint on these parameters. It is this constraint which
allows one to write down a closed system of first order nonlinear ODEs
_without_ needing to specify the functional form of $f(R)$. By its very
construction this formulation is intimately related to the reconstruction
programme and helps us in situations where the reconstruction programme fails
to provide a clear picture. To the best of our knowledge this is the first
ever proposal for a dynamical systems formulation in $f(R)$ which does not
require specifying a functional form of $f(R)$, and therefore links the
powerful dynamical system approach to the reconstruction programme.

The paper is organised as follows:

  * •

Section II consists of an overview of $f(R)$ cosmology and the reconstruction
programme for $f(R)$ gravity.

  * •

Section III introduces our novel dynamical systems formulation in terms of the
cosmographic parameters.

  * •

Section IV applies our dynamical systems formulation to $\Lambda$CDM cosmology
as an example.

  * •

Section V gives a brief overview on the calculation of the density contrast
parameter at the linear regime of cosmological perturbations.

  * •

Section VI details the qualitative comparison between the $\Lambda$CDM model
in GR and alternative late time models in $f(R)$ that are cosmographically
equivalent at the background level.

  * •

We conclude in section VII with a summary of the present work and mention some
other possible cosmologically relevant research problems that can be addressed
along the same line of thought.

Throughout the paper we use the $(-+++)$ signature and unit system with $c=1$
and $\kappa=8\pi G=1$.

##  II $f(R)$ cosmology and reconstruction method

In this section we give a very brief overview of cosmology in $f(R)$ gravity.
An interested reader is referred to the beautiful topical review on the
subject by Sotiriu and Faraoni [38] or De-Fellice and Tsujikawa [39]. $f(R)$
gravity is characterised by the existence of a propagating scalar degree of
freedom $\varphi=F(R)\equiv f^{\prime}(R)$ as apparent from the trace field
equation

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

$T$ being the trace of the energy momentum tensor. General Relativity is the
trivial case of $f(R)$ for which $F(R)=1$ and the scalar degree of freedom is
no longer propagating. This scalar degree of freedom is sometimes dubbed as
the _scalaron_ , a term that we will use hereafter. In the presence of a
perfect fluid with energy density $\rho$ and pressure $P$, the field equations
for $f(R)$ gravity can be expressed as

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

where $F\equiv\frac{df}{dR}$ and we have defined the scalaron energy density
and pressure as

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

The effective equation of state parameter of the universe is defined as

$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(H^{2}+k/a^{2})}\,,$
(5)

and the equation of state of the scalaron is

$w_{R}\equiv\frac{P_{R}}{\rho_{R}}=\frac{\ddot{F}+2H\dot{F}-\frac{1}{2}(RF-f)}{\frac{1}{2}(RF-f)-3H\dot{F}}\,.$
(6)

If the perfect fluid is barotropic, with an equation of state parameter $w$,
then $w_{\rm eff}$ and $w_{R}$ are related to each other via the relation

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

There are two important conditions for physical viability of any $f(R)$
gravity which we just mention below:

  * •

$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
perturbations in the weak gravity limit (also known as the Dolgov-Kawasaki
instability [20]). Therefore, one requires that $f^{\prime\prime}(R)>0$ at
least during the early epoch of matter domination.

Given a particular $f(R)$ theory of gravity, one can analyse the solutions
using the field equations $\eqref{fe}$. However, at times it might be
important to find what form of $f(R)$ gravity can give rise to a particular
desired solution. This problem is systematically addressed by the
_reconstruction_ methods of $f(R)$ gravity [25]. There are in fact various
reconstruction techniques [40], and it is quite possible that for some
particular desired cosmology one of the reconstruction techniques fails while
another succeeds to give back a compact form of $f(R)$. The classical and most
common reconstruction technique attempts to reconstruct an $f(R)$ from a given
solution $a(t)$. The limitation of this technique is that it relies on the
invertibility of the function

$R(t)=6(\dot{H}(t)+2H^{2}(t))\,.$ (8)

If one can find an explicit function $t=g(R)$, then $f(R)$ can be
reconstructed by solving the second order differential equation [40]

$\displaystyle
3H[g(R)]\dot{R}[g(R)]f^{\prime\prime}(R)+\left(3H^{2}[g(R)]-\frac{R}{2}\right)f^{\prime}(R)$
$\displaystyle+\frac{1}{2}f(R)=\rho(R)$ . (9)

Solutions to this equation give a class of $f(R)$ gravity theories for which
the given function $a=a(t)$ is an exact solution. If the solution is not given
in the compact functional form $a(t)$ or the function $R=R(t)$ is non-
invertible, it proves useful to adopt a different reconstruction technique.
The whole reconstruction programme relies on expressing the Ricci scalar $R$
as a function of some cosmological variable such that the function is
invertible. Alternative reconstruction techniques involve expressing the Ricci
scalar as a function of other variables _e.g._
$R=R(a),\,R(H)\,\text{or}\,R(\tau)$ where $\tau\equiv\ln a$. Even if $R(t)$ is
non-invertible, $R=R(a),\,R(H)\,\text{or}\,R(\tau)$ _etc._ might be
invertible. This is what makes one reconstruction method succeed but the
others fail in a particular case.

For example $\Lambda$CDM cosmological evolution is given by the condition [40]

$\dot{a}=\sqrt{\frac{\rho_{0}}{a}+\Lambda a^{2}}\,,$ (10)

where $\rho_{0}$ is a positive constant and $\Lambda>0$ is the cosmological
constant. In this case it proves to be much more useful to express the Ricci
scalar as either a function of the scale factor $a$ (as in [26])

$R=\frac{\rho_{0}}{a^{3}}+4\Lambda\,,$ (11)

or a function of $\tau\equiv\ln a$ (as in [27])

$R=\rho_{0}e^{-3\tau}+4\Lambda\,,$ (12)

both of which are perfectly invertible. In the first case $f(R)$ is
reconstructed from the differential equation [26]

$\displaystyle-3(R-3\Lambda)(R-4\Lambda)f^{\prime\prime}(R)+\left(\frac{R}{2}-3\Lambda\right)f^{\prime}(R)$
$\displaystyle+\frac{1}{2}f(R)=R-4\Lambda$ $\displaystyle\,.$ (13)

In the second case $f(R)$ is reconstructed from the following pair of
equations

$\displaystyle\frac{d^{2}F[R(\tau)]}{d\tau^{2}}+\left(\frac{d\ln
H}{d\tau}-1\right)\frac{dF[R(\tau)]}{d\tau}$ $\displaystyle\hskip
56.9055pt+2\frac{d\ln H}{d\tau}(F[R(\tau)]-1)=0\,,$ (14a) $\displaystyle
f(R)=R+\int(F[R(\tau)]-1)R^{\prime}(\tau)d\tau\,.$ (14b)

If it is possible to write a reconstruction differential equation, there is no
guarantee that a compact form for the general solution can be found. Even in
cases where a compact form can be found, it can involve functions that are too
complicated for any further analytical treatment (_e.g.,_ an analysis of the
cosmological perturbations). For example in the $\Lambda$CDM case both the
Refs.[26] and [27] find the general solution in terms of Hypergeometric
functions, whereas the particular solution is the usual $R-2\Lambda$. The
whole philosophy of this paper is to avoid solving the reconstruction
differential equation, but still be able to qualitatively compare the general
solution $f(R)$ against the particular solution $R-2\Lambda$.

Before moving on to the next section, let us mention here that, as pointed out
in Ref.[41], even though in most cases a compact form $f(R)$ cannot be found
as a solution of the reconstruction differential equation, a series solution
in $R$ can always be written (except in rare situations when the point $R=0$
is an irregular singular point of the differential equation). Numerical
solution can be obtained provided a pair of boundary conditions. For the
$\Lambda$CDM case, setting the boundary conditions $f(R\rightarrow
0)=-2\Lambda$, $F(R\rightarrow 0)=1$ gives back the unique solution
$f(R)=R-2\Lambda$.

##  III A generic dynamical system formulation for $f(R)$ gravity

In terms of expansion normalised dynamical dimensionless variables for $f(R)$
gravity [33, 16, 39],

$\displaystyle x=\frac{\dot{F}}{HF},\quad y=\frac{R}{6H^{2}},\quad
z=\frac{f}{6FH^{2}},$ (15) $\displaystyle\Omega_{=}\frac{\rho}{3FH^{2}},\quad
K=\frac{k}{a^{2}H^{2}}\,,$ (16)

the Friedmann constraint equation becomes

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

Choosing to eliminate $K$ using the Friedmann constraint, the dynamical system
can be expressed as

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

where $\Gamma=\Gamma(R)$ is defined as

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

Given a functional form for $f(R)$, one can _in principle_ invert the relation

$\frac{y}{z}=\frac{RF}{f}$ (20)

to determine $R=R(y/z)$ and correspondingly find $\Gamma=\Gamma(y/z)$, so as
to make the dynamical system (18) autonomous. On the other hand, from the
definition of $\Gamma$ one can write

$xy\Gamma=\frac{\ddot{H}}{H^{3}}-2(1+y)+6(x+z+\Omega)\,.$ (21)

One can keep the term $\ddot{H}/H^{3}$ as an explicit time dependant term in
the system (18), but this makes the system non-autonomous and in general fixed
point analysis cannot be applied unless this term either vanishes or is
constant. This is the approach taken in [42]. However the authors in [42] only
considered the very special cases with $\dot{H}=mH^{3}$ for a constant $m$,
for which the dynamical system becomes autonomous. An alternative approach
towards the dynamical systems formulation of $f(R)$ gravity can be given by
introducing three cosmographic parameters, namely the deceleration, jerk and
snap parameters [43]:

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

which are related to each other by

$\displaystyle j=2q^{2}+q-\frac{dq}{d\tau}\,,$ (23a) $\displaystyle
s=\frac{dj}{d\tau}-j(2+3q)\,.$ (23b)

Since $f(R)$ gravity is a fourth order theory of gravity, i.e., the field
equations contain terms including up to fourth derivatives of the metric,
higher order cosmographic parameters cannot be included in the field
equations. Now one can replace $\Gamma$ in terms of the cosmographic
parameters using

$xy\Gamma=-2(1+y)+6(x+z+\Omega)+j+3q-2\,.$ (24)

The definition of the Ricci scalar provides an additional constraint equation

$\displaystyle y=\frac{\dot{H}}{H^{2}}+2+K=1-q+K\,.$ (25)

The Friedmann constraint (17) can therefore be written as

$z=-x+\Omega-q\,.$ (26)

The two constraint equations (25) and (26) can be used to eliminate $y$ and
$z$ and write the dynamical system in terms of $x,\,\Omega,\,K,\,q,\,j$:

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

The system of equations (27) has a much simpler form as compared to the system
of equations (18) and also does not require us to explicitly specify the
functional form of the underlying $f(R)$ gravity to make the system
autonomous. The last two equations are completely kinematical in nature as
they depend only on _how_ the universe evolves and not at all on _what_
inherent dynamics cause the universe to evolve in that way. In general, a
particular type of time evolution of the universe can be specified in terms of
cosmographic parameters, which allows us to either find the underlying $f(R)$
gravity by solving the reconstruction differential equation [40, 44], or to
write a closed autonomous dynamical system without finding the explicit
solutions. Before proceeding further some comments are in order:

  * •

The present universe is accelerating, _i.e.,_ the deceleration parameter
$q<0$. On the other hand one always has $|K|\leq 1$, _i.e._ $(1+K)\geq 0$.
Therefore from the purely kinematic relation (25) we can conclude that $y>0$
in our present universe. $y=0$ is an invariant submanifold of the dynamical
system (18), _i.e.,_ it divides the whole phase space into two disjoint
regions with $y>0$ ($R>0$) and $y<0$ ($R<0$) respectively. All the physically
relevant dynamics that lead to the present day accelerating universe must take
place within the region $y>0$. Existence of an invariant submanifold also
implies there cannot be a global attractor or repeller, unless it lies on the
invariant submanifold itself.

  * •

As mentioned in the previous section, the physical viability of any $f(R)$
gravity requires $F>0$ throughout the physically relevant region of the phase
space and $F^{\prime}\geq 0$ ($F^{\prime}=0$ corresponding to the special case
$f(R)=R+\Lambda$) at least in the neighbourhood of the fixed point
corresponding to the matter dominated epoch. Eliminating $y$ and $z$ from (24)
using the constraints (25), (26) and then using the definition of $\Gamma$ one
can write

$\frac{1}{y\Gamma}=\frac{6F^{\prime}H^{2}}{F}=\frac{x}{12\Omega-2K-q+j-6}\,.$
(28)

Assuming the condition $F>0$ is met, demanding $F^{\prime}\geq 0$ puts the
following constraint on the phase phase:

$\frac{x}{12\Omega-2K-q+j-6}\geq 0\,.$ (29)

The submanifold $x=0$ correspond to the GR limit ($F^{\prime}=0$).

For any fixed point $P$ given by the coordinates

$P\equiv(x^{*},y^{*},z^{*},\Omega^{*},K^{*},q^{*},j^{*})\,,$

the effective equation of state parameter (5) can be expressed in terms of the
dynamical quantities as

$w_{\rm eff}=\frac{2q^{*}}{3(1+K^{*})}-\frac{1}{3}\,\qquad(K^{*}\neq-1)\,.$
(30)

The relation between different equation of state parameters ((7)) can also be
expressed in terms of the dynamical quantities

$w_{\rm
eff}=w\frac{\Omega}{1+K^{*}}+w_{R}\frac{y^{*}-z^{*}-x^{*}}{1+K^{*}}=w_{R}+(w-w_{R})\frac{\Omega^{*}}{1+K^{*}}\,,$
(31)

where in the last step we have used the Friedmann constraint (26). The case
$K=-1$ corresponds to a Milne solution. In the vicinity of the fixed point one
can _approximately_ reconstruct a form of $f(R)$ from the relation

$\dot{H}=h(H)\equiv-(1+q^{*})H^{2}$ (32)

using the method detailed in Ref.[40]. It should however be kept in mind that
the $f(R)$ form so obtained is valid only during the epoch represented by the
fixed point $P$, and not throughout the whole evolution history that we are
considering.

##  IV Application to $\Lambda$CDM cosmology

As a simple example let us take the observationally successful $\Lambda$CDM
cosmological model, where acceleration of the universe is due to a positive
cosmological constant $\Lambda$ ($w=-1$) and cold dark matter (CDM) is
modelled by a dust-like fluid ($w=0$). The equations of motion are

$\displaystyle 3\left(H^{2}+\frac{k}{a^{2}}\right)=\rho+\Lambda\,,$ (33a)
$\displaystyle 2\dot{H}+3H^{2}+\frac{k}{a^{2}}=\Lambda\,,$ (33b)
$\displaystyle\dot{\rho}+3H\rho=0\,.$ (33c)

As the first and simplest model for late time cosmology, $\Lambda$CDM dynamics
are very well studied ([16, 15, 45]). In terms of the expansion normalised
dimensionless dynamical quantities

$\Omega_{m}=\frac{\rho}{3H^{2}}\,,\qquad\Omega_{\Lambda}=\frac{\Lambda}{3H^{2}}\,,\qquad
K=\frac{k}{a^{2}H^{2}}\,,$ (34)

so that the constraint equation becomes

$\Omega_{m}-K=1-\Omega_{\Lambda}\,.$ (35)

Using the constraint equation, we choose to eliminate $\Omega_{\Lambda}$. The
dynamical system for the $\Lambda$CDM model can then be written as

$\displaystyle\frac{d\Omega_{m}}{d\tau}=-\Omega_{m}(2K-3\Omega_{m}+3)\,,$
(36a) $\displaystyle\frac{dK}{d\tau}=-K(2K-3\Omega_{m}+2)\,.$ (36b)

For any fixed point $Q$ given by the coordinates

$Q\equiv(\Omega_{\Lambda}^{*},\Omega_{m}^{*},K^{*})\,,$

the effective equation of state parameter is

$w_{\rm
eff}=-\frac{\Omega_{\Lambda}^{*}}{1+K^{*}}=\frac{\Omega_{m}^{*}}{1+K^{*}}-1\,\qquad(K\neq-1)\,.$
(37)

Fixed points of the above system, along with their nature under a linear
stability analysis, are listed in Table 1.

Fixed | Coordinates | Stability | Cosmological solution
---|---|---|---
Point | $(\Omega_{m}^{*},K^{*})$ | Nature |
$Q_{1}$ | $(0,0)$ | Attractor |  $\Lambda$-dominated De-Sitter
|  |  | ($w_{\rm eff}=w_{\Lambda}=-1,\,H=const.$)
$Q_{2}$ | $(1,0)$ | Repeller | Matter dominated power law
|  |  | ($w_{\rm eff}=w_{m}=0,\,a~{}t^{2/3}$)
$Q_{3}$ | $(0,-1)$ | Saddle | Milne solution ($a~{}t$)
Table 1: Fixed points of the dynamical system 36.

The whole $\Lambda$CDM cosmology can be specified by the cosmographic
requirement [43, 44]

$j=K+1\,.$ (38)

Instead of attempting to reconstruct an $f(R)$ theory that can exactly mimic
the $\Lambda$CDM evolution history at the observational level, we rather
investigate some of the generic dynamical features of such $f(R)$ theories,
and look for any possible deviations from the $\Lambda$CDM cosmological model.
Using the observational requirement (38) the dynamical system for any generic
value of the spatial curvature parameter $k$ becomes

$\displaystyle\frac{dx}{d\tau}=-x(x-q)+2(x+K+q)-3\Omega+2\,,$ (39a)
$\displaystyle\frac{d\Omega}{d\tau}=-\Omega(x-2q+1)\,,$ (39b)
$\displaystyle\frac{dK}{d\tau}=2qK\,,$ (39c)
$\displaystyle\frac{dq}{d\tau}=2q^{2}+q-K-1\,.$ (39d)

Since CDM is modelled by a dust fluid ($w=0$), Eq.(31) simplifies to

$w_{\rm eff}=w_{R}\left(1-\frac{\Omega^{*}}{1+K^{*}}\right)$ (40)

where a ‘*’ denotes the value of the corresponding quantity at a fixed point.
Fixed points of the above system, along with their nature under a linear
stability analysis, are listed in Table 2.

Fixed | Coordinates | Stability | Cosmological solution
---|---|---|---
Point | $(x^{*},\Omega^{*},K^{*},q^{*})$ | Nature |
$P_{1}$ | $(0,0,0,-1)$ | Saddle | scalaron dominated De-Sitter
|  |  | ($w_{\rm eff}=w_{R}=-1,\,H=const.$)
$P_{2}$ | $(1,0,0,-1)$ | Attractor | scalaron dominated De-Sitter
|  |  | ($w_{\rm eff}=w_{R}=-1,\,H=const.$)
$P_{3}$ | $(0,1,0,\frac{1}{2})$ | Saddle | Matter dominated power law
|  |  | ($w_{\rm eff}=w=0,\,a~{}t^{2/3}$)
$P_{4}$ | $(\frac{5-\sqrt{73}}{4},0,0,\frac{1}{2})$ | Repeller | scalaron dominated power law
|  |  | ($w_{\rm eff}=w_{R}=0,\,a~{}t^{2/3}$)
$P_{5}$ | $(\frac{5+\sqrt{73}}{4},0,0,\frac{1}{2})$ | Saddle | scalaron dominated power law
|  |  | ($w_{\rm eff}=w_{R}=0,\,a~{}t^{2/3}$)
$P_{6}$ | $(0,0,-1,0)$ | Saddle | Milne solution ($a~{}t$)
$P_{7}$ | $(2,0,-1,0)$ | Saddle | Milne solution ($a~{}t$)
Table 2: Fixed points of the dynamical system (39).

$K=0$ is an invariant submanifold of the system (39) on which the spatially
flat dynamics take place. Fixed points that reside on the $K=0$ submanifold
are $P_{1}$, $P_{2}$, $P_{3}$, $P_{4}$ and $P_{5}$. On this invariant
submanifold one can reduce the phase space even more

$\displaystyle\frac{dx}{d\tau}=-x(x-q)+2(x+q)-3\Omega+2\,,$ (41a)
$\displaystyle\frac{d\Omega}{d\tau}=-\Omega(x-2q+1)\,,$ (41b)
$\displaystyle\frac{dq}{d\tau}=2q^{2}+q-1\,.$ (41c)

The system (41) represents an autonomous dynamical system corresponding to all
possible solutions of the reconstruction equation (II) or (14). This allows us
to qualitatively compare the $R-2\Lambda$ gravity against other possible
$f(R)$ theories that produce the same dynamics. We notice that for the
spatially flat case the $q$-equation decouples which leads to two new
invariant submanifolds: a submanifold $q=-1$ consisting of accelerated
cosmological solutions and a submanifold $q=\frac{1}{2}$ consisting of
decelerated cosmological solutions. The fixed points $P_{1}$ and $P_{2}$
reside on the $q=-1$ submanifold and the fixed points $P_{3}$, $P_{4}$ and
$P_{5}$ reside on the $q=\frac{1}{2}$ submanifold. Linear stability analysis
reveals that the “deceleration submanifold” $q=\frac{1}{2}$ is a repelling one
while the “acceleration submanifold” $q=-1$ is an attracting one, which is
consistent with the fact that $P_{2}$ is an attractor while $P_{4}$ is a
repeller.

Another important thing to notice is that unlike the $\Lambda$CDM model where
the matter dominated fixed point is a past attractor, in this case the matter
dominated fixed point is only a saddle, _i.e._ an intermediate epoch. The true
past attractor in this case is the fixed point $P_{4}$ which represents a
scalaron dominated epoch but, surprisingly, the scalaron itself behaves like a
dust as far as cosmological dynamics is concerned.

Given the coordinates of a fixed point, it is a straightforward exercise to
check whether it satisfies the condition (29). In Table 2 the fixed points
$P_{1}$, $P_{3}$ and $P_{6}$ lie on the submanifold $x=0$ that represents the
GR limit. Therefore these three points satisfy the condition (29) trivially.
It can be easily checked that out of the other four fixed points only the past
attractor $P_{4}$ satisfies the constraint. This fixed point corresponds to a
scalaron dominated cosmology but, surprisingly, the scalaron itself behaves
like a dust fluid ($w_{\rm eff}=w_{R}=0$). Therefore this fixed point produces
a time evolution exactly the same as that of a matter dominated epoch in GR
($a\sim t^{2/3}$). In particular it should be noted that the scalaron
dominated De-Sitter future attractor $P_{2}$ does not satisfy the condition
(29). This allows us to conclude that even if there exist possible $f(R)$
models which are able to give rise to cosmological dynamics that are
observationally indistinguishable from $\Lambda$CDM dynamics (at least at the
background level), the $f(R)$-dynamics will inevitably lead to an epoch where
the condition (29) is not met.

Before leaving this section it is worth mentioning that there are indeed
perfectly viable $f(R)$ gravity models that can successfully reproduce
transition from a matter dominated power law evolution epoch to a scalaron
dominated De-Sitter epoch, while always satisfying the condition (29). Two
such examples are _e.g._ $e^{\Lambda R}$($\Lambda>0$) [46] and $R+\alpha
R^{n}$($\alpha,n>0$) [47]. The point we make here is that none of them can
mimic the exact evolution history as produced by the $\Lambda$CDM model, which
till now remains the most observationally fitted model for late time
cosmology. On the other hand if one demands an $f(R)$ model of late time
cosmology to reproduce the same expansion history as that of the $\Lambda$CDM
model, one is bound to end up in a region where the condition (29) is no
longer satisfied.

##  V Behaviour of cosmological perturbations

It has been suggested [39, 48] that one can distinguish between the
$\Lambda$CDM model and an equivalent scalar model by taking into account
observables that depend on cosmological perturbations. In this section we
briefly review this point. For the sake of simplicity we constrain ourselves
to the spatially flat case in this section. Given the fact that present day
observation suggests that our universe is very close to being spatially flat,
it is worthwhile considering this case a little deeper. The perturbation
quantity that is of observational interest in late time cosmology is the
matter density contrast $\delta\equiv\frac{\delta\rho}{\rho}$. The evolution
of this quantity at the sub-Hubble limit is approximately governed by the
second order differential equation [39, 48]

$\frac{d^{2}\delta}{d\tau^{2}}+\left(2+\frac{\dot{H}}{H^{2}}\right)\frac{d\delta}{d\tau}=\frac{2}{3}\frac{\rho\delta}{FH^{2}}\left(\frac{1+\frac{a^{2}}{4k^{2}}\frac{F}{F^{\prime}}}{1+\frac{a^{2}}{3k^{2}}\frac{F}{F^{\prime}}}\right)\,,$
(42)

with $k$ here being the wavenumber of a particular Fourier mode of $\delta$,
not to be confused with the spatial curvature parameter. In terms of the
dynamical variables (15) and utilising the relation (24), the term on the
right hand side of the above equation can be expressed as

$\displaystyle\frac{2}{3}\frac{\rho\delta}{FH^{2}}\left(\frac{1+\frac{a^{2}}{4k^{2}}\frac{F}{F^{\prime}}}{1+\frac{a^{2}}{3k^{2}}\frac{F}{F^{\prime}}}\right)$
$\displaystyle=$ $\displaystyle
2\Omega\left(\frac{x+\frac{3}{2}xy\Gamma\left(\frac{aH}{k}\right)^{2}}{x+2xy\Gamma\left(\frac{aH}{k}\right)^{2}}\right)\delta$
$\displaystyle=$
$\displaystyle\Omega\left(\frac{2x+3(12\Omega-q-5)\left(\frac{aH}{k}\right)^{2}}{x+2(12\Omega-q-5)\left(\frac{aH}{k}\right)^{2}}\right)\delta\,.$

Therefore the sub-horizon perturbation equation can be written as

$\frac{d^{2}\delta}{d\tau^{2}}+(1-q)\frac{d\delta}{d\tau}-\Omega\left(\frac{2x+3(12\Omega-q-5)\left(\frac{aH}{k}\right)^{2}}{x+2(12\Omega-q-5)\left(\frac{aH}{k}\right)^{2}}\right)\delta=0\,.$
(43)

Keeping in mind that $k\gg aH$ in sub-horizon limit, Eq.(43) can be solved at
two different regimes of interest:

  * •

_GR regime_ : $0\lesssim|x|\ll\left(\frac{aH}{k}\right)^{2}$: This corresponds
to the limit when the modification of gravity theory can be safely ignored. In
this limit Eq.(43) reduces to

$\frac{d^{2}\delta}{d\tau^{2}}+(1-q)\frac{d\delta}{d\tau}-\frac{3}{2}\Omega\delta=0\,.$
(44)

  * •

_$f(R)$ regime_: $0\lesssim\left(\frac{aH}{k}\right)^{2}\ll|x|$: This
corresponds to a limit when the effect of the gravity modification cannot be
ignored. In this limit Eq.(43) reduces to

$\frac{d^{2}\delta}{d\tau^{2}}+(1-q)\frac{d\delta}{d\tau}-2\Omega\delta=0\,.$
(45)

A matter dominated power law evolution epoch is characterised by
$q=\frac{1}{2}$ and $\Omega_{m}=1$. For the $\Lambda$CDM model only the GR
regime is relevant. Solving Eq.(44) in this regime we get two modes:

$\displaystyle\delta\sim e^{\tau}\sim a\qquad(\text{Growing mode})$ (46)
$\displaystyle\delta\sim e^{-\frac{3}{2}\tau}\sim
a^{-\frac{3}{2}}\qquad(\text{Decaying mode})$ (47)

For $f(R)$ models however both these regimes are possible. In particular,
depending on the $f(R)$ theory, it is possible for wavenumbers $k$ relevant to
large scale structure observations to transit from the GR regime into the
$f(R)$ regime within the matter dominated regime. Solving Eq.(44) in regime
$II$, we get two modes:

$\displaystyle\delta\sim e^{\frac{1}{4}(\sqrt{33}-1)\tau}\sim
a^{\frac{1}{4}(\sqrt{33}-1)}\quad(\text{Growing mode})\,,$ (48)
$\displaystyle\delta\sim e^{-\frac{1}{4}(\sqrt{33}+1)\tau}\sim
a^{-\frac{1}{4}(\sqrt{33}+1)}\quad(\text{Decaying mode})\,.$ (49)

The growing mode solution is related to the growth rate of large scale
structures in the universe, which is a cosmological observable. Clearly there
is a difference between the growing mode solutions in the GR regime and $f(R)$
regime. In particular the $f(R)$ regime leads to a faster growth of
structures. It is precisely the existence of this regime that leads to an
observable difference between the $\Lambda$CDM model and late time $f(R)$
models.

For a matter perturbation mode of wavelength $\lambda\sim 1/k$, transition
from the GR regime to the $f(R)$ regime occurs at a time which is
approximately given by

$|x|\simeq\left(\frac{aH}{k}\right)^{2}~{}\left(\frac{\lambda}{R_{c}}\right)^{2}\,,$
(50)

where $R_{c}=\frac{1}{aH}$ is the comoving Hubble horizon. The behaviour of
linear matter perturbations therefore depends on the background cosmological
evolution. From the definition of $x$ (15) it is clear that this transition
scale depends on the form of $f(R)$ as well. The observable difference in the
growth rate of structures between $\Lambda$CDM and $f(R)$ based models depends
on this transition scale for a typical perturbation mode $k$ relevant to the
large scale structures of the Universe. In the phase space picture, one can
argue that if this transition occurs very far from the matter dominated fixed
point, then the observable difference from the $\Lambda$CDM model will be
negligible. On the other hand if the transition occurs very near to the matter
dominated fixed point, then the observable difference from the $\Lambda$CDM
model may be significant.

##  VI Comparison between $\Lambda$CDM model and equivalent late time $f(R)$
models

(a)

(b)

(c)

(d)

Figure 1: Projection of the phase portrait on (a) the invariant submanifold
$q=\frac{1}{2}$, (b) the invariant submanifold $q=-1$, (c) the slice
$\Omega=1$, (d) the slice $\Omega=0$. The shaded region represents the region
in which the condition $F^{\prime}>0$ is satisfied (provided $F$ is already
positive).

The phase space corresponding to the dynamical system (41) is 3-dimensional.
One can get a fairly good idea of the phase space dynamics by considering
“projections” of the phase portrait on different planes. Fig.1 shows four sets
of such projections, along with the labelled fixed points. Figs.1(a) and 1(b)
show the projections of the phase portrait on the “deceleration submanifold”
$q=\frac{1}{2}$ and the acceleration submanifold $q=-1$ respectively, whereas
Figs.1(c) and 1(d) show the projections of the phase portrait on the slices
$\Omega=1$ and $\Omega=0$ respectively. In the figures, the line $x=0$
represents the GR limit and the region in which the condition (29) is
satisfied is shaded. As one clearly sees, both the fixed points that
correspond to accelerating epochs lie outside the shaded region. The most
important figure in the context of assessing the observational difference
between $\Lambda$CDM and $f(R)$ models is Fig.1(c), which shows trajectories
emanating from the matter dominated fixed point $P_{3}$ and moving towards the
acceleration submanifold $q=-1$. These trajectories represent possible
evolution routes for the transition from a matter dominated decelerating epoch
to a late time accelerating epoch.

Let us analyse Fig.1(c) in a little more detail. It is clear from the figure
that in the vicinity of the matter dominated fixed point $P_{3}$, the
evolution of the phase trajectories is such that $|x|$ increases with time. On
the other hand, it is straightforward to calculate that

$\frac{d(1/R_{c}^{2})}{d\tau}=-\frac{2q}{R_{c}^{2}}\,,$ (51)

so that $(1/R_{c}^{2})$ is a decreasing function of time near the fixed point
$P_{3}$. This means that even if a relevant perturbation mode of wavelength
$\lambda$ is a-priori within the GR regime ($|x|\ll(\lambda/R_{c}^{2})$), it
is possible during the course of cosmic evolution to achieve the transition
scale given by Eq.(50), beyond which it enters the $f(R)$ regime
($|x|\gg(\lambda/R_{c}^{2})$). Moreover, one can derive from a straightforward
calculation that

$\frac{d}{d\tau}\ln\left(|x|R_{c}^{2}\right)=2q+\frac{(1+q)(1-2q)}{|x|}\frac{d|x|}{d(-q)}\,.$
(52)

Within the region between the submanifolds $q=\frac{1}{2}$ and $q=-1$, $q$ is
monotonically decreasing and the quantity $(1+q)(1-2q)$ is positive
throughout. Moreover, near the fixed point $P_{3}$ both $q$ and
$\frac{d|x|}{d(-q)}$ are positive. Therefore from Eq.(52) one can conclude
that the quantity $|x|R_{c}^{2}$ is increasing with time near the point
$P_{3}$, _i.e.,_ the perturbation modes are moving from the GR regime towards
the $f(R)$ regime. It is also clearly seen from Eq.(52) that the rate of
change of the quantity $|x|R_{c}^{2}$ along a phase trajectory nearby $P_{3}$
directly depends on the slope of the phase trajectory $\frac{d|x|}{d(-q)}$.
The more is the slope $\frac{d|x|}{d(-q)}$, the faster is the increment of the
quantity $|x|R_{c}^{2}$ and consequently quicker is the transition from GR to
$f(R)$ regimes for the perturbation modes. In the paragraph below we discuss
this point by taking as example four characteristic phase trajectories.

In the figure have highlighted four characteristic phase trajectories
emanating from $P_{3}$ in colours red, orange, brown and green respectively.
These trajectories represent segments of four possible cosmic evolutions, each
of which goes though a matter dominated decelerating epoch into an
accelerating phase, and is observationally indistinguishable from the
$\Lambda$CDM model at the background level. How much they will observationally
deviate from the $\Lambda$CDM model at the perturbative level and how
physically viable they will be can be qualitatively assessed by carefully
examining the phase portrait. From our discussion in the the last paragraph we
know that the faster a phase trajectory moves away from the $x=0$ line (_i.e._
has a steeper slope $\frac{d|x|}{d(-q)}$), the quicker is the transition from
GR regime to $f(R)$ regime for some particular characteristic perturbation
mode of wavelength $\lambda\sim 1/k$, and moreover the corresponding cosmic
evolution is expected to observationally deviate from $\Lambda$CDM. On the
other hand since the shaded region represents the region where the condition
(29) is satisfied, the faster a phase trajectory goes out of this region, the
more the corresponding cosmic evolution is expected to encounter the Dolgov-
Kawasaki instability, hence being more physically non-viable.

We note the following:

  * •

The leftmost highlighted trajectory (in red) moves away from both the $x=0$
line much faster compared to the other three highlighted trajectories, and
does not stay within the shaded region at all in the vicinity of the matter
dominated fixed point. Therefore this and nearby trajectories represent a
class of cosmic evolutions that is expected to show significant observational
deviation from the $\Lambda$CDM, while also being severely plagued by the
Dolgov-Kawasaki instability. We can therefore rule out such cosmic evolutions
from being physically viable.

  * •

The highlighted trajectory second from the left (in orange) stays near the
$x=0$ line a little longer compared to the leftmost red trajectory, but spends
only a short time within the shaded region. This and nearby trajectories
represent a class of cosmic evolutions which are observationally closer to the
$\Lambda$CDM model (compared to the red trajectory), but we still cannot
characterise them as physically viable enough.

  * •

The highlighted trajectory second from the right (in brown) always stays
closer to the $x=0$ line, while also being within the shaded region compared
to most other trajectories. This and nearby trajectories represent the most
optimum cosmic evolutions one can get because of their observable closeness
with $\Lambda$CDM and that they avoid the Dolgov-Kawasaki instability.

  * •

Finally, the rightmost highlighted trajectory (in green) stays within the
shaded region longer than all the other three trajectories, but also moves
away from the $x=0$ line almost as fast as the red trajectory. Therefore this
trajectory and those near it represent a class of cosmic evolutions that can
be safely assumed to be free from the Dolgov-Kawasaki instability, but are
expected to show significant observational deviation from $\Lambda$CDM.

Figs. 1(a) and 1(d) show that all trajectories that pass near the saddle point
$P_{3}$ must end up at the future attractor $P_{5}$. Since $P_{5}$ lies
outside the shaded region, as is clear from Figs.1(b) and 1(d), all such
trajectories necessarily end up in a region where the condition (29) is
violated. The trajectories emanating from a region close to $P_{3}$ that do
not remain in the shaded region at all (_e.g._ the red one in Fig.1(c)) can be
immediatly discarded as being physically non-viable. Among the other
trajectories in Fig.1(c) we notice two competing tendencies. The more a
trajectory wants to spend time within the shaded region, the faster it has to
deviate from the $x=0$ line (_e.g._ compare the brown and the green
trajectory). This result, albeit a qualitative one, is nonetheless quite
interesting. Physically this means that the more one demands that a late time
$f(R)$ cosmology should be observationally close to the $\Lambda$CDM model,
the higher the risk is that it is physically non-viable. On the other hand,
the more one tries to construct a physically viable late time $f(R)$
cosmology, the higher the risk of it being observationally different from the
$\Lambda$CDM model. This statement is very generic as it is independent of the
functional form of $f(R)$; the only constraint being it should be
observationally indistinguishable from $\Lambda$CDM at the background level.
That such a very generic statement can be made from a purely qualitative phase
space analysis is truly remarkable.

##  VII Conclusion

In this paper we developed a new dynamical systems framework for studying the
cosmology of $f(R)$ gravity which completely circumvents the reconstruction
programme. This is achieved by using cosmographic parameters to write $f(R)$
cosmology in such a way that it is theory independent. The use of cosmographic
parameters as dynamical variables gives rise to a set of algebraic constraints
on the phase space which are fixed by observations. All earlier autonomous
dynamical system formulations of $f(R)$ gravity require one to specify the
form of $f(R)$ to close the system of dynamical equations. To the best of our
knowledge, this is the first time an autonomous dynamical system formulation
of $f(R)$ gravity is presented that is _model independent_.

By considering the qualitative properties of the resulting phase space and the
growth rate of matter perturbations, we found that models that are
observationally close to $\Lambda$CDM suffer from a higher risk that they
encounter a Dolgov-Kawasaki instability in their future. On the other hand,
demanding that such instabilities should not occur leads to trajectories which
are very different from the standard model.

Other well known examples, such as the Hu and Sawicki $f(R)$ theories of
gravity [21] can also be studied using this approach. As far as we are aware
only a special case of this theory has been considered using a dynamical
systems approach [49], for which it is possible to write the function $\Gamma$
in terms of the dynamical systems variables. This method should allow for a
much more general analysis of the background dynamics of such models and their
parameter space. In Ref.[49] it was found for example that great care must be
taken when fixing the initial conditions. In situations where the cosmological
parameters are chosen to exactly coincide with a LCDM cosmology at $z=0$, the
high redshift behaviour deviated greatly from LCDM and in fact corresponded to
a model dominated by dark-radiation (the effective equation of state was equal
to $1/3$). It was found that viable cosmological evolutions were found if one
rather fixed the initial conditions to coincide with LCDM at high redshift and
evolve the model towards the present time ($z=0$). Care also needs to be taken
to avoid sudden curvature singularities, where the cosmographic parameters
diverge at finite redshift [50]. It will be possible to explore in detail
where in parameter space these pathologies occur using our new Dynamical
Systems approach. All these issues will be addressed in a future paper.

Finally it is worth mentioning that the framework could also be used to
compare different inflationary models that produce a scale-invariant power
spectrum. Moreover, it is also possible to use the same approach to perform
dynamical analysis of other modified gravity theories, _e.g.,_ scalar-tensor
theories.

###### Acknowledgements.

PKSD thanks the First Rand Bank for financial support. 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 Numbers: 123055).

## References

  * [1] S. Perlmutter et al.  Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae.  Astrophys. J., 517:565–586, 1999.
  * [2] Adam G. Riess et al.  Observational evidence from supernovae for an accelerating universe and a cosmological constant.  Astron. J., 116:1009–1038, 1998.
  * [3] John L. Tonry et al.  Cosmological results from high-z supernovae.  Astrophys. J., 594:1–24, 2003.
  * [4] 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.
  * [5] 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.
  * [6] 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.
  * [7] D. N. Spergel et al.  First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Determination of cosmological parameters.  Astrophys. J. Suppl., 148:175–194, 2003.
  * [8] Max Tegmark et al.  Cosmological parameters from SDSS and WMAP.  Phys. Rev. D, 69:103501, 2004.
  * [9] 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.
  * [10] 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.
  * [11] 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.
  * [12] 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.
  * [13] Bhuvnesh Jain and Andy Taylor.  Cross-correlation tomography: measuring dark energy evolution with weak lensing.  Phys. Rev. Lett., 91:141302, 2003.
  * [14] J. P. Ostriker and Paul J. Steinhardt.  Cosmic concordance.  5 1995.
  * [15] Edmund J. Copeland, M. Sami, and Shinji Tsujikawa.  Dynamics of dark energy.  Int. J. Mod. Phys. D, 15:1753–1936, 2006.
  * [16] Sebastian Bahamonde, Christian G. Böhmer, Sante Carloni, Edmund J. Copeland, Wei Fang, and Nicola Tamanini.  Dynamical systems applied to cosmology: dark energy and modified gravity.  Phys. Rept., 775-777:1–122, 2018.
  * [17] Shinji Tsujikawa.  Modified gravity models of dark energy.  Lect. Notes Phys., 800:99–145, 2010.
  * [18] Eleftherios Papantonopoulos (Ed.).  Modifications of Einstein’s Theory of Gravity at Large Distances.  Springer, 2015.
  * [19] Sean M. Carroll, Vikram Duvvuri, Mark Trodden, and Michael S. Turner.  Is cosmic speed - up due to new gravitational physics?  Phys. Rev. D, 70:043528, 2004.
  * [20] A.D. Dolgov and Masahiro Kawasaki.  Can modified gravity explain accelerated cosmic expansion?  Phys. Lett. B, 573:1–4, 2003.
  * [21] Wayne Hu and Ignacy Sawicki.  Models of f(R) Cosmic Acceleration that Evade Solar-System Tests.  Phys. Rev. D, 76:064004, 2007.
  * [22] Alexei A. Starobinsky.  Disappearing cosmological constant in f(R) gravity.  JETP Lett., 86:157–163, 2007.
  * [23] Vinicius Miranda, Sergio E. Joras, Ioav Waga, and Miguel Quartin.  Viable Singularity-Free f(R) Gravity Without a Cosmological Constant.  Phys. Rev. Lett., 102:221101, 2009.
  * [24] T. Multamaki and Iiro Vilja.  Cosmological expansion and the uniqueness of gravitational action.  Phys. Rev. D, 73:024018, 2006.
  * [25] Shin’ichi Nojiri, Sergei D. Odintsov, and Diego Saez-Gomez.  Cosmological reconstruction of realistic modified F(R) gravities.  Phys. Lett. B, 681:74–80, 2009.
  * [26] 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.
  * [27] Jian-hua He and Bin Wang.  Revisiting $f(R)$ gravity models that reproduce $\Lambda$CDM expansion.  Phys. Rev. D, 87(2):023508, 2013.
  * [28] C. B. Collins.  Tilting at cosmological singularities.  Communications in Mathematical Physics, 39(2):131 – 151, 1974.
  * [29] Wainwright J and Ellis G F.  Dynamical System in Cosmology.  Cambridge Univ. Press, 1997.
  * [30] A.A.Coley.  Dynamical Systems and Cosmology.  Springer, 2003.
  * [31] Christian G. Boehmer and Nyein Chan.  Dynamical systems in cosmology.  9 2014.
  * [32] Luca Amendola, Radouane Gannouji, David Polarski, and Shinji Tsujikawa.  Conditions for the cosmological viability of f(R) dark energy models.  Phys. Rev. D, 75:083504, 2007.
  * [33] 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.
  * [34] Jun-Qi Guo and Andrei V. Frolov.  Cosmological dynamics in f(R) gravity.  Phys. Rev. D, 88(12):124036, 2013.
  * [35] Sante Carloni.  A new approach to the analysis of the phase space of $f(R)$-gravity.  JCAP, 09:013, 2015.
  * [36] Artur Alho, Sante Carloni, and Claes Uggla.  On dynamical systems approaches and methods in $f(R)$ cosmology.  JCAP, 08:064, 2016.
  * [37] Jose C. C. de Souza and Valerio Faraoni.  The Phase space view of f(R) gravity.  Class. Quant. Grav., 24:3637–3648, 2007.
  * [38] Thomas P. Sotiriou and Valerio Faraoni.  f(R) Theories Of Gravity.  Rev. Mod. Phys., 82:451–497, 2010.
  * [39] Antonio De Felice and Shinji Tsujikawa.  f(R) theories.  Living Rev. Rel., 13:3, 2010.
  * [40] Sante Carloni, Rituparno Goswami, and Peter K.S. Dunsby.  A new approach to reconstruction methods in $f(R)$ gravity.  Class. Quant. Grav., 29:135012, 2012.
  * [41] Saikat Chakraborty.  Reconstruction method of $f(R)$ gravity for isotropic and anisotropic spacetimes.  Phys. Rev. D, 98(2):024009, 2018.
  * [42] S.D. Odintsov and V.K. Oikonomou.  Autonomous dynamical system approach for $f(R)$ gravity.  Phys. Rev. D, 96(10):104049, 2017.
  * [43] Maciej Dunajski and Gary Gibbons.  Cosmic Jerk, Snap and Beyond.  Class. Quant. Grav., 25:235012, 2008.
  * [44] S. Capozziello, V.F. Cardone, and V. Salzano.  Cosmography of f(R) gravity.  Phys. Rev. D, 78:063504, 2008.
  * [45] Martin Goliath and George F. R. Ellis.  Homogeneous cosmologies with cosmological constant.  Phys. Rev. D, 60:023502, 1999.
  * [46] M. Abdelwahab, S Carloni, and P K.S. Dunsby.  Cosmological dynamics of exponential gravity.  Class. Quant. Grav., 25:135002, 2008.
  * [47] Mohamed Abdelwahab, Rituparno Goswami, and Peter K.S. Dunsby.  Cosmological dynamics of fourth order gravity: A compact view.  Phys. Rev. D, 85:083511, 2012.
  * [48] Seokcheon Lee.  Reconstruction of f(R) gravity models from observations.  Phys. Dark Univ., 25:100305, 2019.
  * [49] Sulona Kandhai and Peter K. S. Dunsby.  Cosmological dynamics of viable f(R) theories of gravity.  10 2015.
  * [50] Andrei V. Frolov.  A Singularity Problem with f(R) Dark Energy.  Phys. Rev. Lett., 101:061103, 2008.