# Reconstruction of the singularity-free $f(\mathcal{R})$ gravity via
Raychaudhuri equations

Gaurav N. Gadbail 0000-0003-0684-9702<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Simran Arora 0000-0003-0326-8945
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute of
Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. P. K. Sahoo 0000-0003-2130-8832
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute
of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Kazuharu Bamba 0000-0001-9720-8817
<EMAIL_ADDRESS>Faculty of Symbiotic Systems Science, Fukushima
University, Fukushima 960-1296, Japan

###### Abstract

We study the bounce cosmology to construct a singularity-free $f(\mathcal{R})$
model using the reconstruction technique. The formulation of the
$f(\mathcal{R})$ model is based on the Raychaudhari equation, a key element
employed in reconstructed models to eliminate singularities. We explore the
feasibility of obtaining stable gravitational Lagrangians, adhering to the
conditions $f_{\mathcal{R}}>0$ and $f_{\mathcal{R}\mathcal{R}}>0$.
Consequently, both models demonstrate stability, effectively avoiding the
Dolgov-Kawasaki instability. Our assessment extends to testing the
reconstructed model using energy conditions and the effective equation-of-
state (EoS). Our findings indicate that the reconstructed super-bounce model
facilitates the examination of a singularity-free accelerating universe for
both phantom and non-phantom phases. However, in the case of the reconstructed
oscillatory bounce model, two scenarios are considered with $\omega=-1/3$ and
$\omega=-2/3$. While the model proves suitable for studying a singular-free
accelerating universe in the $\omega=-1/3$ case, it fails to demonstrate such
behavior under energy conditions for the $\omega=-2/3$ scenario. The
reconstructed models accommodate early-time bouncing behavior and late-time
cosmic acceleration within a unified framework.

Keywords: Reconstruction; $f(\mathcal{R})$ models; Singularity; Super-bounce;
Oscillatory-bounce; Raychaudhari equation

##  I Introduction

Understanding cosmic acceleration is crucial for gaining insights into the
nature of the universe. Several proposed scenarios aim to clarify this
remarkable discovery, falling into two main categories. In one approach,
introducing an exotic component known as dark energy is suggested to be the
driving force behind cosmic acceleration. Up to the present, every independent
cosmological observation agrees with the existence of either vacuum energy
density or a non-zero cosmological constant $\Lambda$ [1]. However, the
straightforward cosmological constant model still faces challenges,
particularly the coincidence problem and the fine-tuning problem [2]. These
issues have prompted extensive efforts to grasp the early stages of the
evolution of the universe. On the other side, achieving accelerated expansion
is possible by modifying the gravitational aspect of the Einstein-Hilbert
action [3]. In line with this approach, numerous extended theories of gravity
have been suggested in recent decades, including a nonminimally coupled scalar
field theory [4], scalar-tensor theory [5], and $f(\mathcal{R})$ gravity [6].
Particularly, $f(\mathcal{R})$ gravity theories, where $f(\mathcal{R})$ is a
function of the Ricci scalar $\mathcal{R}$, are widely utilized in
contemporary cosmology. Higher powers of $\mathcal{R}$ in these models can
explain early-stage accelerated expansion known as inflation, while exploring
negative powers of $\mathcal{R}$ may offer insights into late-time
acceleration in the universe [7]. To validate these theoretical models,
rigorous testing against observational data and various criteria, including
solution stability and local astronomical phenomena, is essential [8]. For a
comprehensive review of $f(\mathcal{R})$ gravity models, covering their
compatibility with observations and theoretical considerations, check the
following references [9]. In summary, $f(\mathcal{R})$ gravity models present
an avenue for alternative gravity theories, necessitating thorough examination
against observational and theoretical benchmarks for validation. They are
appreciated for their capacity to portray the distribution of large-scale
structures in the universe and provide explanations for the observed cosmic
acceleration during later cosmic epochs [10] and the solar system predictions
[11].
The intricate nature of the fourth-order field equations complicates efforts
to comprehend the physics of these theories, making it challenging to obtain
both exact and numerical solutions for comparison with observations.
Nonetheless, recent advancements have tackled this complexity using various
effective techniques. A notable approach, grounded in the theory of dynamical
systems, has shown considerable success [12]. This method offers a simplified
way to obtain exact solutions and provides a qualitative overview of the
overall dynamics of these models. Another interesting approach is the
reconstruction method. In this method, one starts with the assumption that the
expansion history of the universe is precisely understood. The field equations
are then inverted to identify the particular category of $f(\mathcal{R})$
theories associated with a given FLRW model. This method demonstrates numerous
instances where the reconstruction of $f(\mathcal{R})$ gravity aligns with
well-known cosmological phases. These instances encompass the $\Lambda$CDM
epoch, epochs of deceleration/acceleration resembling the presence of phantom
and non-phantom matter, late-time acceleration involving the crossing of the
phantom-divide line, a transient phantom epoch, and scenarios depicting an
oscillating universe [13].

The Raychaudhuri equation is commonly seen as a highly refined equation that
has significantly contributed to exploring Einstein’s general theory of
relativity. Raychaudhuri formulated his well-known equation, characterized by
its pure geometric nature and absence of a direct association with the
Einstein field equation [14]. The Raychaudhuri equation addresses the concept
of geodesic congruence, aiming to explore and understand the temporal
evolution of a congruence of geodesics as perceived by an adjacent member
within the exact congruence. In essence, the equation elucidates the
relationship between the rate of change in the scalar expansion of a test
particle and the gravitational impact associated explicitly with the energy
momentum of that test particle.

Considering the geometric perspective, we aim to create a mathematically
versatile analytical solution for different epochs within the $f(\mathcal{R})$
gravity model. Since the Raychaudhuri equation is entirely based on geometry,
we find it applicable to our case. We have opted for two types of cosmological
bouncing scale factors: super bounce and oscillatory bounce-to establish the
functional form of $f(\mathcal{R})$ gravity. Here are some advantages
associated with the super-bounce and oscillatory bouncing models: the
avoidance of singularities is a key benefit; the cyclic nature of the model
provides an alternative to the standard cosmological paradigm; compatibility
with modified gravity theories contributes to the avoidance of singularities;
the cyclic evolution may naturally address cosmological problems such as the
horizon and flatness problems; and the super-bounce and oscillatory bounce
model provides an alternative to inflationary cosmology. Utilizing the
Raychaudhuri equation, we examine effective energy conditions and phases
dominated by dark energy, aiming to formulate a general yet practical
$f(\mathcal{R})$ gravity model as an analytical solution. It is important to
note that we are aware of the prior derivation of the Raychaudhuri equation in
the context of $f(\mathcal{R})$ gravity using a different approach.

In recent years, there has been a notable interest in cosmological bouncing
solutions as an alternative to the conventional Big Bang singularity, aiming
to address the unnaturalness associated with the origin of the universe [15,
16]. Nonsingular bouncing cosmologies, in particular, have gained attention
for their potential to resolve the singularity problem linked to the prevalent
inflationary scenario in early universe cosmology [17]. Achieving a bouncing
solution for the spatially flat FLRW metric within general relativity (GR)
typically involves introducing matter components that violate the null energy
condition (NEC), such as ghost fields, ghost condensates, or Galileons [18].
However, for those seeking to avoid exotic matter components and still realize
a bouncing solution in the spatially flat FLRW metric, modified gravity
emerges as a viable option [19]. It is reasonable to anticipate modifications
to general relativity in the high-curvature regime near a curvature
singularity. Various attempts have been made to realize a bouncing scenario
within the $f(\mathcal{R})$ theories of gravity [20].
In subsequent work, we use a reconstruction technique to analyze different
cosmological bounce scenarios within the framework of $f(\mathcal{R})$
gravity. In particular, we develop $f(\mathcal{R})$ gravity models that
exhibit late-time cosmic acceleration and a bounce in the early universe. We
look into these models’ stability. We investigate two kinds of bounces: super-
bounce and oscillatory bounce. The reconstruction approach first uses the
explicit form of the scale factor related to a specific bounce type, enabling
the calculation of the Hubble parameter $H(t)$. Next, we directly apply the
Raychaudhuri equation to determine the $f(\mathcal{R})$ gravity model for the
two cases. In the first case, we get a combination of powers of the Ricci
scalar $\mathcal{R}$, which represents an eternally accelerated model.
Combinations of hypergeometric functions appear in the second scenario. We
also assess the viability of the $f(\mathcal{R})$ models against various
cosmological requirements.
The paper is structured as follows: In Section II, a concise overview of
$f(\mathcal{R})$ gravity theory is presented. Moving to Section III, we derive
the Raychaudhuri equation for $f(\mathcal{R})$ gravity within the framework of
the FLRW metric. Section IV is dedicated to the reconstruction of suitable
gravitational Lagrangians $f(\mathcal{R})$ that can accurately replicate the
cosmic evolution observed in super-bounce and oscillatory bounce models.
Additionally, we assess the validity criteria of these $f(\mathcal{R})$
Lagrangians. The cosmological aspects of the reconstructed models are examined
in Section V. Finally, Section VI provides the conclusions drawn from the
study.

##  II $f(\mathcal{R})$ gravity

We begin by recalling the straightforward modification of the Einstein-Hilbert
action, where $\mathcal{R}$ is replaced with a general function of the Ricci
scalar, namely [6]

$S_{\mathcal{R}}=\int\sqrt{-g}\left(f(\mathcal{R})+\mathcal{L}_{m}\right)d^{4}x,$
(1)

where $f(\mathcal{R})$ is an analytic function of the Ricci scalar
$\mathcal{R}$ and $\mathcal{L}_{m}$ is the matter Lagrangian. Here we take
$8\pi\,G=1$. The governing field equations for a general $f(\mathcal{R})$
theory is derived by varying the action (1) with respect to the metric alone,
yielding the following field equations

$f_{\mathcal{R}}R_{\mu\nu}-\frac{f}{2}g_{\mu\nu}-\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\right)f_{\mathcal{R}}=T_{\mu\nu},$
(2)

where $f_{\mathcal{R}}=df/d\mathcal{R}$ and
$T_{\mu\nu}\equiv-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g})\mathcal{L}_{m}}{\delta
g^{\mu\nu}}$ is the stress-energy tensor, $R_{\mu\nu}$ is the Ricci tensor,
and $\square=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the d’Alembert operator.

In $f(\mathcal{R})$ gravity, the field equations can be rearranged in terms of
the Einstein tensor, comprising the higher-order curvature contributions

$G_{\mu\nu}=\frac{1}{f_{\mathcal{R}}}\left(T_{\mu\nu}+T^{eff}_{\mu\nu}\right),$
(3)

where

$T^{eff}_{\mu\nu}=\frac{f-\mathcal{R}f_{\mathcal{R}}}{2}g_{\mu\nu}+\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\right)f_{\mathcal{R}},$
(4)

where the effective energy-momentum tensor $T^{eff}_{\mu\nu}$ behaves as a
source for the resulting geometry. This source term is usually referred to as
the curvature fluid. Henceforth, when the stress-energy tensor for the
curvature fluid is exactly zero, we recover GR. The difference between this
equation and Einstein’s equations, at least technically, is that the presence
of $f$ denotes a non-minimal coupling. In this formulation, the effective
gravitational coupling will not be a constant. To ensure the classical and
quantum stability of a theory in the physically relevant domain, the following
requirements on the derivatives of $f$ are crucial, that is,
$f_{\mathcal{R}}>0$ and $f_{\mathcal{RR}}>0$. Indeed, the two conditions
guarantee that gravity is attractive and that the effective gravitational
constant is positive, respectively. It also ensures that the graviton is not a
ghost and ensures the avoidance of the Dolgov-Kawasaki instability [21, 8].

##  III Raychaudhuri equations

In this section, we study the positive contributions of the Raychaudhuri
equations for timelike geodesics, which guarantee the attractive character of
the gravitational interaction in $f(\mathcal{R})$ theories. For a timelike
congruence with velocity vector $u^{\mu}$, the Raychaudhuri equation is
written as [14]

$\frac{d\theta}{d\tau}=-\frac{1}{3}\theta^{2}+\nabla_{\mu}a^{\mu}-\sigma_{\mu\nu}\sigma^{\mu\nu}+w_{\mu\nu}w^{\mu\nu}-R_{\mu\nu}u^{\mu}u^{\nu},$
(5)

where $\tau$ is the affine parameter, $\theta=\nabla_{\mu}u^{\mu}$ is the
expansion scalar, $u_{\mu}$ is the timelike velocity vector,
$\sigma_{\mu\nu}=\nabla_{(_{\nu}u_{\mu})}-\frac{1}{3}h_{\mu\nu}\,\theta+a_{(_{\nu}u_{\mu})}$
is the shear, $h_{\mu\nu}$ is the spatial metric,
$w_{\mu\nu}=\nabla_{[_{\nu}u_{\mu}]}-a_{[_{\nu}u_{\mu}]}$ is the rotation of a
congruence of timelike geodesics, and $a^{\mu}=u^{\nu}\nabla_{\nu}\,u^{\mu}$
is the acceleration vector.

The last term on the right side of the Raychaudhuri equation, i.e.,
$-R_{\mu\nu}u^{\mu}u^{\nu}$, emphasizes the role of space-time geometry and is
independent of the vector field’s derivative. Using the field equations (2)
for $f(\mathcal{R})$ theory, this term may be expressed as follows:

$R_{\mu\nu}u^{\mu}u^{\nu}=\frac{1}{f_{\mathcal{R}}}\left(T_{\mu\nu}+\frac{f}{2}g_{\mu\nu}+\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\right)f_{\mathcal{R}}\right)u^{\mu}u^{\nu}.$
(6)

Consequently, this term has more general importance than the further terms in
(5). This term may be understood geometrically as a mean curvature in the
direction of $u$. The Raychaudhuri equation becomes more straightforward if we
make the following assumptions: (i) Geodesics are the congruence of time-like
curves. Then, $a_{\mu}=0$ is found along the geodesics. (ii) Timelike geodesic
congruence to be hyper-surface orthogonal entails zero vorticity, i.e., the
fourth term of the right-hand side of (5) vanishes by virtue of the Frobenius
Theorem. The Raychaudhuri equation (5) becomes

$\frac{d\theta}{d\tau}=-\frac{1}{3}\theta^{2}-\sigma_{\mu\nu}\sigma^{\mu\nu}-R_{\mu\nu}u^{\mu}u^{\nu}.$
(7)

We now assume the Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetime with a
metric of the type that is homogeneous, isotropic, and spatially flat at the
background level

$ds^{2}=-dt^{2}+a^{2}(t)(dx^{2}+dy^{2}+dz^{2}),$ (8)

where $a(t)$ is the scale factor that is a function of cosmic time $t$. For
such a metric and a matter distribution of a perfect fluid, given by
$T_{\mu\nu}=(p+\rho)u_{\mu}u_{\nu}+p\,g_{\mu\nu}$, the Raychaudhuri equation
(5) takes the form,

$\frac{\ddot{a}}{a}=\frac{1}{f_{\mathcal{R}}}\left(\frac{f(\mathcal{R})}{6}+H\,\dot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}-\frac{\rho}{3}\right).$
(9)

So far, we have refrained from making any presumptions regarding the equation
of state for the fluid distribution. Instead, the Raychaudhuri equation
eliminates the fluid pressure $p$ using field equations (3).

##  IV Bouncing solution in $f(\mathcal{R})$ gravity

In this section, we examine the feasibility of getting suitable gravitational
Lagrangians $f(\mathcal{R})$ capable of reproducing the cosmic evolution given
by super-bounce and oscillatory bounce models. For a review of these bounces,
one can check References [22].

###  IV.1 Super-bounce

The super-bouncing cosmology is distinguished by a power-law scale factor,
which is first explored in [23]. This bouncing model is used to construct a
universe that collapses and rebirths without a singularity [24]. The Super-
bounce scale factor and corresponding Hubble parameter are written as

$a(t)\sim(t_{s}-t)^{2/c^{2}},\,\,\,\,\,\,\,H(t)=-\frac{2}{c^{2}(t_{s}-t)},$
(10)

where $t_{s}$ represents the time at which the bounce occurs.
Thus, we express the quantity $\frac{\ddot{a}}{a}$ and the Ricci scalar
$\mathcal{R}$ as

$\displaystyle\frac{\ddot{a}(t)}{a(t)}=\frac{2}{c^{2}}\left(\frac{2}{c^{2}}-1\right)\frac{1}{(t_{s}-t)^{2}},$
(11)
$\displaystyle\mathcal{R}=6\frac{\dot{a}^{2}}{a^{2}}+6\frac{\ddot{a}}{a}=\frac{12}{c^{2}}\left(\frac{4}{c^{2}}-1\right)\frac{1}{(t_{s}-t)^{2}}.$
(12)

Further, rewriting the equation (12) gives

$\frac{1}{(t_{s}-t)^{2}}=\frac{c^{4}}{12(4-c^{2})}\mathcal{R}.$ (13)

By using the above expressions, the equation (9) gives the following
expression

$\mathcal{R}^{2}f_{\mathcal{R}\mathcal{R}}+\frac{2-c^{2}}{2c^{2}}\mathcal{R}f_{\mathcal{R}}-\frac{4-c^{2}}{2c^{2}}f=-\frac{4-c^{2}}{c^{2}}\rho.$
(14)

Assuming the barotropic equation of state $p=\omega\rho$ yields the energy
density using the conservation equation as $\rho=\rho_{0}\,a^{-3(1+\omega)}$.
Henceforth, we obtain the energy density in terms of $\mathcal{R}$ (using Eqs.
(10) and (12)) given by

$\rho=\rho_{0}\,\left(\frac{c^{4}\,\mathcal{R}}{12(4-c^{2})}\right)^{\frac{3(1+\omega)}{c^{2}}}.$
(15)

Consequently, the equation (14) takes the form

$\mathcal{R}^{2}f_{\mathcal{R}\mathcal{R}}+\frac{2-c^{2}}{2c^{2}}\mathcal{R}f_{\mathcal{R}}-\frac{4-c^{2}}{2c^{2}}f=\mathcal{X}_{1}\,\mathcal{R}^{\gamma},$
(16)

where
$\mathcal{X}_{1}=-\frac{\rho_{0}\,(4-c^{2})}{c^{2}}\left(\frac{c^{4}}{12(4-c^{2})}\right)^{\frac{3(1+\omega)}{c^{2}}}$.
In fact, one can now obtain the solution to the above differential equation
given by

$f(\mathcal{R})=c_{1}\,\mathcal{R}^{\alpha}+c_{2}\,\mathcal{R}^{\beta}+\mathcal{X}_{2}\,\mathcal{R}^{\gamma},$
(17)

where $c_{1}$ and $c_{2}$ are integrating constants and

$\displaystyle\alpha=-\frac{2-3c^{2}+\sqrt{4+20c^{2}+c^{4}}}{4c^{2}},$
$\displaystyle\beta=\frac{-2+3c^{2}+\sqrt{4+20c^{2}+c^{4}}}{4c^{2}},$
$\displaystyle\gamma=\frac{3(1+\omega)}{c^{2}},$
$\displaystyle\mathcal{X}_{2}=\frac{2\mathcal{X}_{1}\,c^{4}}{c^{4}-c^{2}(13+9\omega)+6(4+7\omega+3\omega^{2})}.$

Figure 1: Variation of model parameters $\alpha$, $\beta$, $\gamma$ with
respect to $c^{2}$.

One can notice that the solution for $f(\mathcal{R})$ incorporates three
distinct powers of $\mathcal{R}$. The behavior of $\alpha$, $\beta$, and
$\gamma$ with respect to $c^{2}$ is depicted in FIG. 1. Notably, for
$0<c^{2}<4$, $\alpha$ consistently is negative, while $\beta$ remains
positive. However, the value of $\gamma$ is contingent upon the equation of
state parameter $\omega$. In regions where $\omega>-1$ (non-phantom phase),
$\gamma$ invariably assumes positive values, whereas in the phantom phase,
where $\omega<-1$, $\gamma$ consistently adopts negative values.

####  IV.1.1 Strong Energy Condition

The possibility of accelerating expansion without the requirement for dark
energy (or an inflaton field) is a significant part of the motivation for
$f(\mathcal{R})$ gravity. A straightforward way to see this is to define the
effective pressure and effective energy density of the geometry as [9]

$\displaystyle H^{2}=\frac{\rho_{eff}}{3},$ (18)
$\displaystyle\frac{\ddot{a}}{a}=-\frac{\rho_{eff}+3p_{eff}}{6},$ (19)

where $p_{eff}$ and $\rho_{eff}$ are the effective effective pressure and
energy density, respectively.

In this context, we may also examine the effective energy condition by
calculating the quantity $\rho_{eff}+3p_{eff}$, which in this case is given by

$\rho_{eff}+3p_{eff}=\frac{12}{c^{4}}(c^{2}-4)\frac{1}{(t_{s}-t)^{2}}.$ (20)

The negative nature of the solution for $f(\mathcal{R})$ within the range
$0<c^{2}<4$ signifies a violation of the strong energy condition. This aligns
with expectations for an eternally accelerated model.

####  IV.1.2 Viability Criteria

We aim to verify the theoretical consistency of the relevant $f(\mathcal{R})$
gravity model. Any viable $f(\mathcal{R})$ model must adhere to the conditions
$f_{\mathcal{R}}>0$ and $f_{\mathcal{R}\mathcal{R}}>0$ [25, 26, 9]. The
stability of the model hinges on the latter requirement, while the former
ensures the positivity of the effective gravitational constant. Our obtained
model satisfies both conditions $f_{\mathcal{R}}>0$ and
$f_{\mathcal{R}\mathcal{R}}>0$ (evidenced in FIG. 2).

We find the viability of the model at low and high curvature for two different
phases.

  * •

For non-phantom phase :- In the expression of $f(\mathcal{R})$, the first term
dominates and $\mathcal{R}$ approaches $0$ as $\alpha$ is negative. Thus, for
the model to be valid at low curvature, it is necessary to have $c_{1}=0$,
ensuring the viability. Additionally, the model remains viable at high
curvatures.

  * •

For phantom phase :- In this case, the dominance shifts towards the first and
third terms as $\mathcal{R}\to 0$ due to negative values of $\alpha$ and
$\gamma$. As a result, for the model to be viable at low curvature, we must
have $c_{1}=0$, and $\mathcal{X}_{2}=0$. Further, the model remains viable
when curvature is high.

  * •

For $\Lambda$CDM ($\omega=-1$) :- According to the functional form obtained
for $f(\mathcal{R})$, the third term remains constant as
$\mathcal{X}_{2}=\frac{2\rho_{0}}{c^{2}}$. Further, $c_{1}=0$ satisfy the
viability criteria at low curvature. This simplifies our model to
$f(\mathcal{R})=c_{2}\,\mathcal{R}^{\beta}+\Lambda(=\mathcal{X}_{2})$.
Additionally, in Fig. 1, as $c^{2}$ increases, the parameter $\beta$ tend
towards $1$, aligning our model with the $\Lambda$CDM
($f(\mathcal{R})=\mathcal{R}+\Lambda$).

Figure 2: Behavior of $f_{\mathcal{R}}$ and $f_{\mathcal{R}\mathcal{R}}$ vs.
$\mathcal{R}$ for different values of $c^{2}$.

###  IV.2 Oscillatory bounce

A study on an oscillatory bouncing model in cosmology has gained attention for
its potential advantages and alternative approach to understanding the
universe. The Oscillatory bouncing cosmology is distinguished by a periodic
scale factor and the corresponding Hubble parameter, given by

$a(t)=A\,\sin^{2}\left(\frac{Bt}{t_{*}}\right),\,\,\,\,\,\,\,H(t)=\frac{2B}{t_{*}}\,\cot\left(\frac{Bt}{t_{*}}\right).$
(21)

where $t_{*}>0$ is some reference time, $A>0$ and $B>0$ are dimensionless
constants. This bouncing model characterizes a cyclic universe as a series of
contractions and expansions (periodically). With a choice of scale factor, the
model displays two different bouncing patterns. Each cycle begins with a “Big
Bang”, finishes with a ”big crunch”, and then begins again with a “Big Bang”.

The first scenario involves a singularity occurring regularly within each
cycle, characterized by the scale factor reaching zero and the Hubble
parameter becoming singular, indicating a Big Crunch or Big Bang singularity.
The second bounce occurs when the universe reaches its maximum considerable
extent, ceasing expansion and initiating its contracting toward Big Crunch
singularity. For more details, one can review these References [27, 28, 29,
22].

We define the quantity $\frac{\ddot{a}}{a}$ and the Ricci scalar $\mathcal{R}$
as

$\displaystyle\frac{\ddot{a}}{a}=\frac{2B^{2}}{t^{2}_{*}}\left[\cot^{2}\left(\frac{Bt}{t_{*}}\right)-1\right],$
(22)
$\displaystyle\mathcal{R}=6\frac{\dot{a}^{2}}{a^{2}}+6\frac{\ddot{a}}{a}=\frac{12B^{2}}{t^{2}_{*}}\left[3\,\cot^{2}\left(\frac{Bt}{t_{*}}\right)-1\right].$
(23)

One can rewrite Eq. (21) as

$\sin^{2}\left(\frac{Bt}{t_{*}}\right)=\frac{a}{A},\,\,\,\,\,\cot^{2}\left(\frac{Bt}{t_{*}}\right)=\frac{A}{a}-1.$
(24)

Furthermore, using Eqs.(23) and (24), we obtain

$\frac{A}{a}=\frac{t_{*}^{2}}{36\,B^{2}}\mathcal{R}+\frac{4}{3}.$ (25)

Hence, introducing all the above expressions, Eq. (9) turns out to be

$\left(\mathcal{R}+\frac{12\,B^{2}}{t_{*}^{2}}\right)\left(\mathcal{R}+\frac{48\,B^{2}}{t_{*}^{2}}\right)\,f_{\mathcal{R}\mathcal{R}}\\\
+\frac{1}{2}\left(\mathcal{R}-\frac{24\,B^{2}}{t_{*}^{2}}\right)\,f_{\mathcal{R}}-\frac{3}{2}\,f=-3\,\rho.$
(26)

We further use the barotropic equation of state $p=\omega\rho$ to obtain the
energy density as $\rho=\rho_{0}\,a^{-3(1+\omega)}$. Now, by using Eqs.(10)
and (12), we express the energy density in terms of $\mathcal{R}$ as

$\rho=\rho_{0}\left(\frac{t_{*}^{2}}{36A\,B^{2}}\right)^{3(1+\omega)}\,\left(\mathcal{R}+\frac{48\,B^{2}}{t_{*}^{2}}\right)^{3(1+\omega)}.$
(27)

The equation (26) takes the form

$\left(\mathcal{R}+\frac{12\,B^{2}}{t_{*}^{2}}\right)\left(\mathcal{R}+\frac{48\,B^{2}}{t_{*}^{2}}\right)\,f_{\mathcal{R}\mathcal{R}}+\\\
\frac{1}{2}\left(\mathcal{R}-\frac{24\,B^{2}}{t_{*}^{2}}\right)\,f_{\mathcal{R}}-\frac{3}{2}\,f=\mathcal{K}\,\left(\mathcal{R}+\frac{48\,B^{2}}{t_{*}^{2}}\right)^{3(1+\omega)},$
(28)

where
$\mathcal{K}=-3\rho_{0}\left(\frac{t_{*}^{2}}{36A\,B^{2}}\right)^{3(1+\omega)}$.
To simplify the above differential equation, we put
$\lambda=\frac{12B^{2}}{t_{*}^{2}}$ to obtain

$(\mathcal{R}+\lambda)(\mathcal{R}+4\lambda)\,f_{\mathcal{R}\mathcal{R}}+\frac{1}{2}(\mathcal{R}-2\lambda)\,f_{\mathcal{R}}-\frac{3}{2}\,f\\\
=\mathcal{K}\,\left(\mathcal{R}+4\lambda\right)^{3(1+\omega)}.$ (29)

The corresponding homogeneous equation above can be transformed into the
hypergeometric differential equation by the substitution
$x=-\frac{1}{3\lambda}(\mathcal{R}+\lambda)$, as

$x(1-x)\frac{d^{2}f}{dx^{2}}-\frac{1}{2}(x+1)\frac{df}{dx}+\frac{3}{2}f=0.$
(30)

The solution of the homogeneous differential equation is

$f_{H}(x)=c_{1}\,{}_{2}F_{1}\left(1,-\frac{3}{2},-\frac{1}{2},x\right)+c_{2}\,x^{3/2},$
(31)

where ${}_{2}F_{1}$ is a hypergeometric function. To find the particular
solutions, we categorize them based on the EoS parameter $\omega$. In these
cases, we observe that the real-valued particular solutions are only obtained
when $\omega>-1$ (non-phantom phase), while the solutions become imaginary
when $\omega<-1$ (phantom phase).

Table 1: The particular integral for differential Eq. (29) for different values of the EoS parameter $\omega$. 3$(1+\omega)$ | EoS parameter | Particular Integral | solutions
---|---|---|---
$3(1+\omega)=0$ | $\omega=-1$ | $f_{p}(\mathcal{R})=2\rho_{0}(=\Lambda)$ | Accepted
$3(1+\omega)=1$ | $\omega=-\frac{2}{3}$ | $f_{p}(\mathcal{R})=-\mathcal{K}(\mathcal{R}+2\lambda)$ | Accepted
$3(1+\omega)=2$ | $\omega=-\frac{1}{3}$ | $f_{p}(\mathcal{R})=\mathcal{K}\left(6\mathcal{R}^{2}-\frac{200}{27}\lambda\,\mathcal{R}-16\lambda^{2}\right)-\frac{128\,\mathcal{K}}{81}\sqrt{3\lambda}\,(\mathcal{R}+\lambda)^{3/2}\,\tan^{-1}\left(\sqrt{\frac{\mathcal{R}+\lambda}{3\lambda}}\right)$ | Accepted
$3(1+\omega)<0$ | $\omega<-1$ | Imaginary solutions | Rejected

####  IV.2.1 Strong Energy Condition

We can additionally analyze the effective energy condition to assess an
oscillatory bounce. This entails evaluating the quantity
$\rho_{eff}+3p_{eff}$, which for this case is

$\rho_{eff}+3p_{eff}=-\frac{12B^{2}}{t_{*}}\left[\cot^{2}\left(\frac{B\,t}{t_{*}}\right)-1\right].$
(32)

The expression above is negative for
$\cot^{2}\left(\frac{B\,t}{t_{*}}\right)>1$, i.e.,
$-\frac{(2n+1)\pi}{4}<\frac{B\,t}{t_{*}}<\frac{(2n+1)\pi}{4}$ and thus
violates the strong energy condition, which is expected for an eternally
accelerated model.

####  IV.2.2  Viability Criteria

We seek to verify the theoretical consistency of the relevant $f(\mathcal{R})$
gravity model. Our derived models in all instances satisfy both the conditions
$f_{\mathcal{R}}>0$ and $f_{\mathcal{R}\mathcal{R}}>0$, as depicted in Figs.
3). We assess the viability of the model at low and high curvature regimes for
various phases as follows:

  * •

For $3(1+\omega)=0$ :- For this scenario, $f(R)\to 0$ as $R\to 0$ when we have
$c_{1}\lambda^{3/2}+\frac{\pi\,c_{2}}{27\sqrt{3}\lambda}=\frac{2\rho_{0}}{3}$.
As a result, at low curvature, we must have
$c_{1}=\frac{54\lambda\,\rho_{0}-\sqrt{3}\pi\,c_{2}}{81\lambda^{5/2}}$ to meet
the viability criterion, and hence, the model will consistently remain viable
at large curvature.

  * •

For $3(1+\omega)=1$ :- Here, we obtain that, $f(R)\to 0$ as $R\to 0$ when
$c_{1}\,\lambda^{3/2}+\frac{\pi\,c_{2}}{27\sqrt{3}\lambda}=2\lambda\mathcal{K}$.
Consequently, we must have
$c_{1}=\frac{162\mathcal{K}\lambda^{2}-\sqrt{3}\pi\,c_{2}}{81\lambda^{5/2}}$
at low curvature to be viable. Similarly, the model is viable at large
curvature also.

  * •

For $3(1+\omega)=2$ :- In this scenario, $f(R)\to 0$ as $R\to 0$ for
$c_{1}\lambda^{3/2}-\frac{\pi\left(64\lambda^{3}\,\mathcal{K}-3\,c_{2}\right)}{81\,\sqrt{3}\lambda}=\frac{16}{3}\lambda^{2}\,\mathcal{K}$.
So, at low curvature, we get
$c_{1}=\frac{-3\pi\sqrt{3}\,c_{2}+1296\mathcal{K}\lambda^{3}+64\pi\sqrt{3}\mathcal{K}\lambda^{3}}{243\lambda^{5/2}}$
for the viability requirement. However, the model will consistently remain
viable at high curvature.

(a) (b)

(c) (d)

(e) (f)

Figure 3: Plots (a) and (d) illustrate the behavior of $f_{\mathcal{R}}$ and
$f_{\mathcal{R}\mathcal{R}}$ for $\omega=-1$. Plots (b) and (e) depict the
behavior of $f_{\mathcal{R}}$ and $f_{\mathcal{R}\mathcal{R}}$ for
$\omega=-2/3$. Plots (c) and (f) exhibit the behavior of $f_{\mathcal{R}}$ and
$f_{\mathcal{R}\mathcal{R}}$ for $\omega=-1/3$.

##  V Cosmological aspects of reconstructed models

Returning to our discussion, assuming that the stress-energy tensor is
$T_{\mu\nu}=(p+\rho)u_{\mu}u_{\nu}+p\,g_{\mu\nu}$, and inserting the flat FLRW
metric (8) in the field equations (2) yields [9, 30]

$\displaystyle 3H^{2}$ $\displaystyle=$
$\displaystyle\frac{\kappa}{f_{\mathcal{R}}}\left[\rho+\frac{Rf_{\mathcal{R}}-f}{2}-3H\dot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}\right],$
(33) $\displaystyle 2\dot{H}+3H^{2}$ $\displaystyle=$
$\displaystyle-\frac{\kappa}{f_{\mathcal{R}}}\left[p+\dot{\mathcal{R}}^{2}f_{\mathcal{R}\mathcal{R}\mathcal{R}}+2H\dot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}\right.$
(35)
$\displaystyle\left.+\ddot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}+\frac{1}{2}(f-\mathcal{R}f_{\mathcal{R}})\right].$

An important part of the motivation for $f(\mathcal{R})$ gravity is that it
can result in accelerated expansion without the need for dark energy. To
exemplify this, let us define the effective energy density and pressure as
follows

$\displaystyle\rho_{eff}$ $\displaystyle=$
$\displaystyle\frac{Rf_{\mathcal{R}}-f}{2f_{\mathcal{R}}}-\frac{3H\dot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}}{f_{\mathcal{R}}},$
$\displaystyle p_{eff}$ $\displaystyle=$
$\displaystyle\frac{\dot{\mathcal{R}}^{2}f_{\mathcal{R}\mathcal{R}\mathcal{R}}+2H\dot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}+\ddot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}+\frac{1}{2}(f-\mathcal{R}f_{\mathcal{R}})}{f_{\mathcal{R}}},$

where $\rho_{eff}$ has to be non-negative in a spatially flat FLRW spacetime.

The effective equation of state parameter
$\omega_{eff}=\frac{p_{eff}}{\rho_{eff}}$ for $f(\mathcal{R})$ can be
represented as

$\omega_{eff}=\frac{\dot{\mathcal{R}}^{2}f_{\mathcal{R}\mathcal{R}\mathcal{R}}+2H\dot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}+\ddot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}+\frac{1}{2}(f-\mathcal{R}f_{\mathcal{R}})}{\frac{Rf_{\mathcal{R}}-f}{2}-3H\dot{\mathcal{R}}f_{\mathcal{R}\mathcal{R}}}.$
(36)

Energy conditions (ECs) serve as fundamental tools for understanding the
geodesics of the universe. These conditions can be derived using the well-
known Raychaudhury equations, as described in Eq. (5). The ECs are
characterized as [30, 31]

  * •

Strong Energy Conditions (SEC) : $\rho_{eff}+3p_{eff}\geq 0$

  * •

Weak Energy Conditions (WEC) : $\rho_{eff}\geq 0$,
$p_{eff}+\rho_{eff}\geq 0$

  * •

Null Energy Conditions (NEC) : $p_{eff}+\rho_{eff}\geq 0$

  * •

Dominant Energy Conditions (DEC) : If $\rho_{eff}\geq 0$, $\rvert
p_{eff}\rvert\leq\rho_{eff}$.

One idea is to let us think about the singularity-free universal beginning,
which is the universal bouncing scenario. Many suggestions in the literature
suggested avoiding this singularity via quantum features, but these lack the
dependability to integrate into gravitational theory. So, at this stage,
gravitational theories provide a specific technique for testing the validity
of both the bounce model and its own. The null energy condition (NEC) is a
technique for accomplishing the objective. Furthermore, it has been
demonstrated that in the setting of GR, violating NEC is exceedingly difficult
for local-field models. As a result, effective field theories allow for the
detection of NEC violations and the possibility of a non-singular bounce [32].
However, in standard general relativity, the strong energy condition is
closely tied to the idea that gravity is always attractive. If this condition
is violated in modified gravity, it could lead to repulsive gravitational
effects, potentially contributing to accelerated cosmic expansion. The
dominant energy conditions are linked to the stability of solutions in
gravitational theories, and its validation suggests that the spacetime
configurations and solutions remain stable and do not exhibit instabilities
that could lead to unexpected behavior or violations of causality.

Figure 4: Energy conditions with respect to time using $\beta=0.98$,
$\gamma=0.85$ ($\gamma>0$ represent the non-phantom case), $c_{2}=-1.75$ and
$\mathcal{X}_{2}=0.5$ for the Super-bounce model with $0<c^{2}<4$. Figure 5:
Energy conditions with respect to time using $\beta=0.2$, $\gamma=-0.7$
($\gamma<0$ represent the phantom case), and $c_{2}=-0.75$ for the Super-
bounce model with $0<c^{2}<4$.

We test our reconstructed models through the energy conditions. The energy
conditions for the reconstructed from super-bounce are displayed in Figs. 4
(non-phantom case) and 5 (phantom case). In both figures, we can see that the
NEC is violated near the bouncing point and satisfied after a certain point.
The violation of NEC near the singular point derives the bouncing nature of
the universe. Also, the SEC is violated after a certain point toward the late-
time universe, and this could mimic the effects of dark energy, the mysterious
force believed to be driving the accelerated expansion of the universe.
Furthermore, the DEC is satisfied after a certain point toward the late-time
universe. Hence, we conclude that our reconstructed model is useful for
studying the singular free accelerating universe.

Figure 6: Behavior of the effective EoS parameter versus cosmic time $t$ for
the Super-bounce model.

The effective equation of state parameter shows the phantom dark energy
behavior after a certain point toward the late-time universe (see Fig. 6) for
the same value of the model parameters. Phantom dark energy, characterized by
$\omega_{eff}<-1$, leads to accelerated cosmic expansion. The negative
pressure associated with phantom dark energy causes a repulsive gravitational
effect, driving the universe to expand at an increasing rate.

The energy conditions for the oscillatory bounce model are plotted in Figs. 8
(for $\omega=-1/3$) and 9 (for $\omega=-2/3$). In Fig. 8, only the strong
energy condition (SEC) is violated, and the other two NEC and WEC are
satisfied. Hence, the results of oscillatory bounce model for the case
$\omega=-2/3$ tested through the energy conditions show the acceleration of
the universe, but it fails to show the singular free universe. In Fig. 9, we
can see that the NEC is violated near the bouncing point and satisfies its
condition after a certain point. The violation of NEC near the singular point
derives the bouncing nature of the universe. In this figure, the SEC is also
violated, and it can show the negative pressure regime of the universe.
Furthermore, the DEC is satisfied at late times. Hence, we conclude that our
reconstructed oscillatory model for the case $\omega=-1/3$ is useful for
studying the singular free accelerating universe. The effective equation of
state parameter shows the phantom dark energy behavior (see Fig. 7) for the
same values of the model parameters.

Figure 7: Behavior of the effective EoS parameter with respect to time for the
oscillatory bounce. Figure 8: Energy conditions with respect to time for the
oscillatory bounce for $\omega=-2/3$ using the value of parameters $\lambda=1$
(or $2$ or $3$), $\mathcal{K}=1$, and $c_{2}=-1$.  Figure 9: Energy conditions
with respect to time for the oscillatory bounce for $\omega=-1/3$ using the
value of parameters $\lambda=1$ (or $2$ or $3$), $\mathcal{K}=-0.005$, and
$c_{2}=-1$.

##  VI conclusion

The conventional cosmological model posits that the universe originated from a
singularity, commonly referred to as the Big Bang. Nonetheless, this model
abounds with several challenges, including the horizon problem, flatness
problem, original structure problem, entropy problem, and the persisting issue
of a singularity. In order to surmount these obstacles, a rapid phase of
expansion immediately following the Big Bang was deemed necessary to establish
a universe characterized by uniformity, flatness, and smoothness. Addressing
these concerns, Alan Guth formulated an inflationary theory, which has proven
remarkably successful in explaining various observational features of the
universe. Despite these challenges, bouncing cosmology offers a promising
avenue for exploring alternative scenarios to the conventional Big Bang model
and has the potential to address some of the longstanding puzzles in
cosmology. Ongoing theoretical and observational efforts are aimed at further
elucidating the feasibility and implications of bouncing cosmology in
understanding the origin and evolution of the universe.

In this study, we have formulated $f(\mathcal{R})$ models specifically
developed for addressing the challenge of the initial singularity in both
super-bounce and oscillatory bounce scenarios. The construction of these
$f(\mathcal{R})$ models are rooted in the Raychaudhari equation, a key element
incorporated in reconstructed models designed to circumvent singularities.
Furthermore, we examined the feasibility of getting suitable gravitational
Lagrangians for $f(\mathcal{R})$ capable of reproducing the cosmic evolution
given by super-bounce and oscillatory bounce models.

Super-bounce :- We have obtained the power-law type $f(\mathcal{R})$ model for
this bouncing scenario. We have investigated this model for two phases: non-
phantom and phantom. Our analysis of the super-bounce model involves
evaluating the quantity $\rho_{eff}+3p_{eff}$, leading to the constraint
$0<c^{2}<4$ (as expected for an accelerated model). We find the viability of
the model at low and high curvature for both non-phantom and phantom cases,
imposing certain constraints on model parameters. Furthermore, we have
verified the viability of our obtained model through two conditions: firstly,
ensuring $f_{\mathcal{R}}>0$ to maintain a positive effective gravitational
constant, which also ensures that the graviton is not a ghost [21], and
secondly, confirming $f_{\mathcal{R}\mathcal{R}}>0$ to ensure stability, i.e.,
the avoidance of the Dolgov-Kawasaki instability.

Oscillatory bounce :- In this analysis, we can see in Table 1 that our
particular solution varies for different phases of the universe, suggesting
the importance of employing corresponding models to study each phase
effectively. We further investigate the effective strong energy condition by
evaluating the quantity $\rho_{eff}+3p_{eff}$. The effective strong energy
condition is violated when $\cot^{2}\left(\frac{B\,t}{t_{*}}\right)>1$, i.e.,
$-\frac{(2n+1)\pi}{4}<\frac{B\,t}{t_{*}}<\frac{(2n+1)\pi}{4}$, aligning with
expectations for an endlessly accelerating model. Similar to the previous
model, we found the viability of our model at low and high curvature and
confirmed that the obtained model can be stable.

Our reconstructed models underwent rigorous testing based on energy
conditions, yielding noteworthy results. These include the avoidance of
classical cosmological singularities linked to the Big Bang and the Big
Crunch. Additionally, our model demonstrates the potential to drive the
accelerated expansion of the universe. We observed the violation of NEC near
the bouncing point, followed by the validation after a certain threshold. This
violation of the NEC near the singularity underscored the bouncing nature of
the universe. Moreover, the SEC is violated towards the late-time universe,
potentially depicting the accelerated expansion. The validation of DEC
indicated that the spacetime configurations and solutions maintain stability,
avoiding instabilities that might result in unexpected behavior or violations
of causality. Furthermore, the effective EoS shows the negative pressure
associated with phantom dark energy causes a repulsive gravitational effect,
driving the universe to expand at an increasing rate.

In this study, a notable advancement is the explicit analytical derivation of
the $f(\mathcal{R})$ gravity model. The models accommodate early-time bouncing
behavior and late-time cosmic acceleration within a unified framework. Such
findings provide valuable insight into constructing models that depict the
early universe and offer a pathway for understanding the cosmological
mechanism driving the ongoing accelerated expansion of the universe. It would
be interesting to find further motivated reconstructed models in Palatini’s
theories and other aspects of quantum gravity phenomenology using the same
approach in future works. We also aim to construct a viable non-singular
bounce in $f(\mathcal{R})$ gravity using the above technique in an
inflationary context.

## Data Availability

There are no new data associated with this article.

## Acknowledgments

GNG acknowledges University Grants Commission (UGC), New Delhi, India for
awarding Junior Research Fellowship (UGC-Ref. No.: 201610122060). PKS
acknowledges Science and Engineering Research Board, Department of Science and
Technology, Government of India for financial support to carry out Research
project No.: CRG/2022/001847 and IUCAA, Pune, India for providing support
through the visiting Associateship program. The work of KB is supported in
part by the JSPS KAKENHI Grant Number 21K03547 and 23KF0008.

## References

  * [1] A.G. Riess et al., Astron. J., 116, 1009-1038 (1998); E. Komatsu et al., Astrophys. J. Suppl., 192, 18 (2011); S. Cao and Z-H. Zhu, Phys. Rev. D, 90, 083006 (2014); S. Cao, Z-H. Zhu, R. Zhao, Phys. Rev. D, 84, 023005 (2011); E. R. Peterson et al., Astrophys. J., 938, 112 (2022); A. G. Riess et al., Astrophys. J., 938, 36 (2022).
  * [2] B.P. Abbott et al., Natur., 551, 85, (2017); X-G. Zheng et al., Astrophys. J., 892, 103 (2020); V. Sahni and A. Starobinsky, Int. J. Mod. Phys. D, 09, 373-443 (2000).
  * [3] S.M. Carroll et al., Phys. Rev. D, 71, 063513 (2005).
  * [4] O. Bertolami and P.J. Martins, Phys. Rev. D, 61, 064007 (2000); N. Banerjee and D. Pavon, Phys. Rev. D, 63, 043504 (2001).
  * [5] D. Alonso et al., Phys. Rev. D, 95, 063502 (2017); L. Heisenberg, R. Kase, S. Tsujikawa, Phys. Rev. D, 98, 024038 (2018).
  * [6] S. Capozziello, Int. J. Mod. Phys. D, 11, 483 (2002); S. Capozziello, S. Carloni, and A. Troisi, Recent Res. Dev. Astron. Astrophys., 1, 625 (2003); S. Capozziello, V. F. Cardone, S. Carloni, and A. Troisi, Int. J. Mod. Phys. D, 12, 1969 (2003); S. Capozziello, V. F. Cardone, and A. Trosi, Phys. Rev. D, 71, 043503 (2005); S. Carloni, P. K. S. Dunsby, S. Capozziello, and A. Troisi, Class. Quantum Grav., 22, 4839 (2005); S. Nojiri and S.D. Odintsov, Phys. Rev. D, 68, 123512 (2003).
  * [7] A.A. Starobinsky, Phys. Lett. B, 91, 99 (1980).
  * [8] A.D. Dolgov and M. Kawasaki, Phys. Lett. B, 573, 1 (2003); S. Nojiri and S.D. Odintsov, Phys. Rev. D, 74, 086005 (2006).
  * [9] S. Nojiri and S.D. Odintsov, Int. J. Geom. Meth. in Mod. Phys., 4, 115-145 (2007); T.P. Sotiriou and V. Faraoni, Rev. Mod. Phys., 82, 451 (2010); A. De Felice and S. Tsujikawa, Living Rev. Rel., 13, 3 (2010); S. Nojiri and S. D. Odintsov, Phys. Rept., 505, 59-144 (2011); S. Capozziello and M. De Laurentis, Phys. Rept., 509, 167-321 (2011).
  * [10] Y-S Song, W. Hu, I. Sawicki, Phys. Rev. D, 75, 044004 (2007); A. Abebe, A. de la Cruz-Dombriz, P. K. S. Dunsby, Phys. Rev. D, 88, 044050 (2013); S. Capozziello, V. F. Cardone, V. Salzano, Phys. Rev. D, 78, 063504 (2008); T. Koivisto, Phys. Rev. D, 73, 083517 (2006).
  * [11] T. Chiba, Phys. Lett. B, 575, 1–3, (2003); T. Chiba, T. L. Smith, and A. L. Erickcek, Phys. Rev. D, 75, 124014 (2007); W. Hu and I. Sawicki, Phys. Rev. D, 76, 064004 (2007); S. A. Appleby and R. A. Battye, Phys. Lett. B, 654, 7–12 (2007); S. Tsujikawa, Phys. Rev. D, 77, 023507 (2008); G. Cognola et al., Phys. Rev. D, 77, 046009 (2008); S. Capozziello and S. Tsujikawa, Phys. Rev. D, 77, 107501 (2008); C. Negrelli et al., Phys. Rev. D, 101, 064005 (2020).
  * [12] S.D. Odintsov and V. K. Oikonomou, Phys. Rev. D, 96, 104049 (2017); S. Carloni, J. Cosmol. Astropart. Phys., 09, 013 (2015); S.D. Odintsov and V.K. Oikonomou, Phys. Rev. D, 98, 024013 (2018).
  * [13] S. Carloni et al., Class. Quantum Grav., 29, 135012 (2012); S.D. Odintsov and V.K. Oikonomou, Ann. Phys., 388, 267-275 (2018); S. Chakraborty, Phys. Rev. D, 98, 024009 (2018); S. Nojiri, S.D. Odintsov, D. S-Gomez, Phys. Lett. B, 681, 74-80 (2009).
  * [14] A. Raychaudhuri, Phys. Rev., 98, 1123 (1955); J. Ehlers, Akad. Wiss. Lit. Abhandl. Math.-Nat.Kl., 11, 793 (1961); J. Ehlers, Gen. Relativ. Gravit., 25, 1225 (1993).
  * [15] S. Nojiri, S.D. Odintsov, V.K. Oikonomou, Phys. Rep., 692, 1-104 (2017); M. Novello and S.E.P. Bergliaffa, Phys. Rept., 463, 127 (2008).
  * [16] J. de Haro et al., Phys. Rep, 1034, 1-114 (2023).
  * [17] V. Mukhanov, Physical foundations of cosmology (Cambridge University Press, Oxford), (2005); K. Bamba, S.D. Odintsov, Symmetry, 7, 220-240 (2015); S.D. Odintsov, V.K. Oikonomou, I. Giannakoudi, F.P. Fronimos and E.C. Lymperiadou, Symmetry, 15(9), 1701 (2023).
  * [18] C. Deffayet, G.E-Farese, A. Vikman, Phys. Rev. D, 79, 084003 (2009); T. Kobayashi, Phys. Rev. D, 94, 043511 (2016).
  * [19] S.D. Odintsov and V.K. Oikonomou, Phys. Rev. D, 92, 024016 (2015); V.K. Oikonomou, Phys. Rev. D, 92, 124027 (2015); K. Bamba A.N. Makarenko, A.N. Myagky and S.D. Odintsov, Phys. Lett. B, 732, 349-355 (2014); M.A. Skugoreva, A.V. Toporensky, Eur. Phys. J. C, 80, 1054 (2020); L. Raul Abramo, I. Yasuda, P. Peter, Phys. Rev. D, 81, 023511 (2010); J.K. Singh et al., JHEP 03, 191 (2023); K. Bamba, A.N. Makarenko, A.N. Myagky, S. Nojiri and S.D. Odintsov, J. Cosmol. Astropart. Phys., 01, 008 (2014).
  * [20] E.N. Saridakis, S. Banerjee, R. Myrzakulov, Phys. Rev. D, 98, 063513 (2018); K. Bhattacharya and S. Chakraborty, J. Cosmol. Astropart. Phys., 2016, 030 (2016); S.D Odintsov, V.K. Oikonomou, T. Paul, Class. Quantum Grav., 37, 235005 (2020); R. Daniel, M. Campbell, C. van de Bruck and P. Dunsby, J. Cosmol. Astropart. Phys., 06, 030 (2023).
  * [21] S.A. Appleby, R.A. Battye, A.A. Starobinsky, J. Cosmol. Astropart. Phys., 1006, 005 (2010).
  * [22] M. Caruana, G. Farrugia, J. L. Said, Eur. Phys. J. C, 80, 640 (2020); A. de la C-Dombriz et al., Phys. Rev. D, 97, 104040 (2018).
  * [23] M. Koehn, J. Lehners, and B.A. Ovrut, Phys. Rev. D, 90(2), 025005 (2014).
  * [24] V.K. Oikonomou, Astrophys. Space Sci., 359(1), 30 (2015).
  * [25] I. Sawicki, W. Hu, Phys. Rev. D, 75, 127502 (2007).
  * [26] A. Silvestri, M. Trodden, Rep. Prog. Phys., 72, 096901 (2009).
  * [27] R. Tolman, Relativity, Thermodyanamics and Cosmology (Dover publications, New York, 1934)
  * [28] P. J. Steinhardt and N. Turok, Phys. Rev. D, 65, 126003 (2002).
  * [29] J. Khoury, P. J. Steinhardt, and N. Turok, Phys. Rev. Lett., 92, 031302 (2004).
  * [30] J. Santos et al., Phys. Rev. D, 76, 083513 (2007); S. Capozziello, F. S. N. Lobo, J. P. Mimoso, Phys. Lett. B, 730, 280-283 (2014).
  * [31] S. M. Carroll, Addison Wesley, 101, 102 (2004); S. Kar and S. Sengupta, Pramana, 69, 49–76 (2007).
  * [32] D. Larson et al., Astrophys. J. Suppl. Ser., 192, 16 (2011); R.R. Caldwell, Phys. Lett. B, 545, 23-29 (2002); U. Alam et al., Mon. Not. Roy. Astron. Soc., 354, 275 (2004); V. Onemli and R. Woodard, Phys. Rev. D, 70, 107301 (2004).