#
Dynamical analysis of coupled curvature-matter scenario in viable $f(R)$ dark
energy models at de-sitter phase

Anirban Chatterjee<EMAIL_ADDRESS>&<EMAIL_ADDRESS>Ramakrishna Mission Vivekananda Educational and Research Institute, Belur
Math, Howrah 711202, West-Bengal, India

(August 31, 2024)

###### Abstract

We explore the interaction between dark matter and curvature-driven dark
energy within viable $f(R)$ gravity models, employing the phase-space analysis
approach of linear stability theory. By incorporating an interacting term,
denoted as $\mathcal{Q}=\alpha H\tilde{\rho}_{\rm
m}\left(\frac{\kappa^{2}}{3H^{2}}\rho_{\rm curv}+1\right)$, into the
continuity equations of both sectors, we examine dynamics of two $f(R)$
gravity models that adhere to local gravity constraints and fulfill
cosmological viability criteria. In de-sitter phase, our investigation reveals
modifications of critical points compared to the conventional form, attributed
to the introduced interaction. Through a comprehensive phase-space analysis,
we illustrate the trajectories near critical points and outline the
constraints on the phase space based on cosmological and cosmographic
parameters.

##  1 Introduction

Observational evidence of redshift and luminosity distance for type Ia
supernovae events [1, 2] has been instrumental in validating the transition of
the universe from a decelerated phase to its presently accelerated phase
during late-time evolution. These evidences find support in the examination of
temperature anisotropies within the cosmic microwave background data from the
WMAP mission [3, 4] and the identification of baryon acoustic oscillations
[5]. In accordance with the general theory of relativity, the mechanism
propelling accelerated expansion is attributed to a theoretical energy density
component characterized by negative pressure, commonly known as dark energy.
The inquiry into the origin and properties of dark energy, responsible for the
ongoing cosmic acceleration, persists as a significant unsolved enigma in
contemporary cosmology. Late-time acceleration of cosmos eludes explanation
through the standard equation of state, denoted as $\omega=p/\rho$, where $p$
and $\rho$ signify the pressure and energy density of conventional universe
constituents like radiation and matter. Evidently, there exists a necessity
for an as-yet-unidentified component marked by negative pressure, yielding an
equation of state $\omega<-1/3$, to elucidate the observed accelerated
expansion.

Einstein initially introduced the $\Lambda g_{\mu\nu}$ cosmological constant
term in his general relativity equations for a static universe. However, he
later discarded it in favor of Hubble’s evidence of an expanding universe. In
the late 20th century, the term regained significance due to its potential to
explain late-time cosmic acceleration, leading to the development of the
$\Lambda$-CDM model, incorporating cold dark matter. Unfortunately, this model
faces challenges such as the coincidence [6] and fine-tuning [7] problems,
prompting the exploration of alternative dark energy models. One model
category involves field-theoretic approaches to dark energy, modifying
Einstein’s field equations by introducing a scalar field as an additional
universe component, distinct from matter and radiation. This category includes
quintessence [8, 9, 10, 11, 12, 13, 14, 15, 16] and $k-$essence [17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29] models. Another model class focuses on
altering the geometric aspects of Einstein’s equations, particularly the
Einstein-Hilbert action, to address late-time cosmic acceleration. These
models, detailed in references [30, 31, 32, 33, 34, 35, 36], encompass $f(R)$
gravity, scalar-tensor theories, Gauss-Bonnet gravity, and braneworld models
of dark energy. In many scenarios, the universe is assumed to be homogeneous
and isotropic on large scales, as described by the Friedmann-Lemaître-
Robertson-Walker (FLRW) metric with a time-dependent scale factor, denoted as
$a(t)$, and a curvature constant. Some approaches also consider an
inhomogeneous universe, using a perturbed FLRW metric to represent the
background spacetime. Dynamical systems techniques are powerful for studying
cosmic evolution in both generic cosmological models and specific solutions.
Recent applications include assessing the stability of scalar field dark
energy models [37, 38, 39, 40, 41] and exploring $f(R)$ gravity [42, 43, 44,
45]. Comprehensive investigations into dynamical systems within $f(R)$ gravity
theories are documented in various research endeavors, such as [45, 46, 47,
48, 49]. These studies aimed to categorize accurate cosmological models among
different $f(R)$ gravity models. They contributed to developing
‘cosmologically viable’ $f(R)$ gravity models that align with observed cosmic
phenomena and meet local gravity constraints [42, 43, 44, 45]. These models
offer cosmic acceleration without needing a cosmological constant or
introducing dark energy as an additional field, relying on the curvature-
driven dynamics described by the modified field equations with the function
$f(R)$ of the Ricci scalar.

Here, we have investigated the impact of interacting curvature-driven dark
energy and (dark) matter in viable $f(R)$ gravity models at de-sitter phase.
Theoretical framework results from a modified geometric part and a matter
component, ensuring minimal coupling with the geometry. The modified field
equations, resembling standard field equations, involve the Einstein tensor
and a total energy-momentum tensor. The latter comprises a curvature component
determined by $f(R)$ and its derivatives, and a modified matter component
derived by scaling the usual energy-momentum tensor from Einstein’s equation.
In FLRW spacetime background, conservation of total energy-momentum tensor
leads to a continuity equation for the comprehensive fluid, including matter
and curvature components. We introduce interactions between matter and
curvature components with a source term $\mathcal{Q}$ in their continuity
equations, representing the rate of energy exchange between them. The
interaction term has been chosen for this work is, $\mathcal{Q}=\alpha
H\tilde{\rho}_{\rm m}\left(\frac{\kappa^{2}}{3H^{2}}\rho_{\rm curv}+1\right)$.
In previous articles, authors have opted for alternative interaction terms,
such as those exclusively related to matter or a combination of matter and
curvature. But here, the motivation behind choosing this particular form of
interaction is to study the system’s dynamics in the presence of both additive
and multiplicative types of matter, matter-curvature interactions. Here,
$\alpha$ indicates the coupling strength of the interaction term, and
$\tilde{\rho}_{\rm m}$ and $\rho_{\rm curv}$ represent their respective energy
densities.

Incorporating interactions between matter and curvature introduces the
parameter $\alpha$ into this analytical framework, intertwined with parameters
defining the $f(R)$ function. This coupling parameter significantly influences
evolutionary dynamics in the matter-curvature coupled scenario. We use
dynamical analysis to explore this interaction in viable $f(R)$ gravity
models, setting up autonomous equations with dynamical variables. Identifying
fixed points and assessing their stability through linear stability analysis
is crucial for understanding aspects of cosmic evolution influenced by
curvature-matter interactions. We focus on two $f(R)$ models, the generalized
$\Lambda$-CDM model and the power-law model, both capable of driving cosmic
acceleration while remaining cosmologically viable. To grasp the dynamics of
the curvature-matter system in two $f(R)$ gravity models within the de-sitter
phase, we streamlined the approach from 3D to 2D variables and analyzed them
using linear stability theory. From these phase-space plots, we can conclude
the allowed region of each model at the de-sitter phase in terms of stability,
cosmographic parameters, and the viability condition of critical energy
density of matter.

The paper is organized as follows: In sec. [2], we have established a
theoretical framework to address the interaction between the curvature and
matter sectors in the background of a flat FLRW metric. In sec. [3], we have
formulated the autonomous equations of the dynamical system, incorporating the
interaction between matter and curvature-driven dark energy in the framework
of $f(R)$ gravity. We have also discussed the two specific $f(R)$ models
considered for this study. In sec. [4], we have examined this interacting
scenario in the context of two different types of viable $f(R)$ models through
the dynamical analysis approach and put some extra constraints from stability
and other criteria. Finally, We have summarized our conclusions in sec. [5].

##  2 Framework of curvature-matter interaction in $f(R)$ gravity

Action of minimally coupled curvature-matter sector [30, 31, 45] in the
context of $f(R)$ theory of gravity can be written as,

$\displaystyle\mathcal{S}=\frac{1}{2\kappa^{2}}\int d^{4}x\sqrt{-g}f(R)+\int
d^{4}x\leavevmode\nobreak\ \mathcal{L}_{M}(g_{\mu\nu},\phi_{M})\,,$ (1)

Here, $\kappa^{2}=8\pi G$ and $f(R)$ is any function of the Ricci scalar $R$,
and ${\cal L}_{M}$ is the Lagrangian density for matter component of the
universe, with $g$ as the determinant of the metric tensor $g{\mu\nu}$. In
metric formalism, the modified field equation is obtained by varying this
action with respect to the field $g_{\mu\nu}$ as,

$\displaystyle F(R)R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}f(R)+g_{\mu\nu}\square
F(R)-\nabla_{\mu}\nabla_{\nu}F(R)=\kappa^{2}\tilde{T}_{\mu\nu}^{(\rm M)}\,.$
(2)

Here, $F(R)$ is the derivative of $f$ with respect to $R$, and
$\tilde{T}_{\mu\nu}^{(\rm M)}$ is the stress-energy tensor for matter
components given by

$\displaystyle\tilde{T}_{\mu\nu}^{(\rm
M)}\equiv-{2\over\sqrt{-g}}{\delta(\sqrt{-g}{\cal
L}_{M})\over\delta(g^{\mu\nu})}\,.$ (3)

Eq. (2) may be reconstructed as [30, 45]

$\displaystyle G_{\mu\nu}\equiv
R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}=\kappa^{2}\left(T_{\mu\nu}^{(\rm
M)}+T_{\mu\nu}^{(\rm curv)}\right)\equiv\kappa^{2}T_{\mu\nu}^{(\rm tot)}\,,$
(4)

$G_{\mu\nu}$ depicts the Einstein tensor and

$\displaystyle\kappa^{2}T_{\mu\nu}^{(\rm curv)}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{F}\left[\frac{1}{2}(f-RF)g_{\mu\nu}+\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\right)F\right]$
(5) $\displaystyle\mbox{and}\quad T_{\mu\nu}^{(\rm M)}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{F}\tilde{T}_{\mu\nu}^{(\rm M)}$ (6)

The reformulated field equation Eq. (4) describes the FLRW universe with a
fluid represented by the comprehensive energy-momentum tensor
$T_{\mu\nu}^{(\rm tot)}\equiv T_{\mu\nu}^{(\rm curv)}+T_{\mu\nu}^{(\rm M)}$.
The influence of the function $f(R)$ is evident in the distinct components
$T_{\mu\nu}^{(\rm curv)}$ and $T_{\mu\nu}^{(\rm M)}$, forming
$T_{\mu\nu}^{(\rm tot)}$. The former, determined by $f(R)$ and its higher
derivatives, vanishes for $f(R)=R$, acting as a fluid-equivalent portrayal of
the curvature function. The latter modifies $\tilde{T}_{\mu\nu}^{(\rm M)}$
using a functional multiplier of $1/F$, reflecting gravity modifications from
$f(R)$ in the revised stress-energy tensor for matter component. In a flat
FLRW spacetime background, the unmodified stress-energy tensor
$\tilde{T}_{\mu\nu}^{(\rm M)}$, representing matter as a perfect fluid,
becomes diagonal with entries derived from combined energy densities
$\tilde{\rho}_{\rm m}$. The matter component is considered non-relativistic.
Consequently, the modified stress-energy tensor
$T^{(M)}_{\mu\nu}=(1/F)\tilde{T}_{\mu\nu}^{(\rm M)}$ corresponds to a fluid
with energy density $\rho_{\rm m}$ and zero pressure, where $\rho_{\rm
m}\equiv\tilde{\rho}_{\rm m}/F$. Positivity of $\rho_{m}$ is maintained due to
the inherent positivity of $\tilde{\rho_{m}}$ and $F$. The ‘00’ and ‘$ii$’
components of the modified field equation Eq. (4) yield modified versions of
the Friedmann equations under these considerations;

$\displaystyle H^{2}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{3}\Big{[}\rho_{\rm m}+\rho_{\rm
curv}\Big{]}\,,$ (7) $\displaystyle\dot{H}+H^{2}$ $\displaystyle=$
$\displaystyle-\frac{\kappa^{2}}{6}\Big{[}\rho_{\rm m}+\rho_{\rm curv}+3P_{\rm
curv}\Big{]}\,,$ (8)

Here, $H$ denotes the Hubble parameter, $\rho_{\rm curv}$ and $P_{\rm curv}$
can be expressed as,

$\displaystyle\rho_{\rm curv}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\kappa^{2}F}\left(\frac{RF-f}{2}-3H\dot{R}F^{\prime}\right)$
(9) $\displaystyle P_{\rm curv}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\kappa^{2}F}\left(\dot{R}^{2}F^{\prime\prime}+2H\dot{R}F^{\prime}+\ddot{R}F^{\prime}+\frac{1}{2}(f-RF)\right)\,,$
(10)

Where, ′ signifies the derivative with respect to $R$. In a flat FLRW
spacetime, the stress-energy tensor $T_{\mu\nu}^{\rm{curv}}$ for curvature-
fluid is equivalent to that of an ideal fluid, with energy density and
pressure represented by $\rho_{\rm curv}$ and $P_{\rm curv}$ as defined in
Eqs. (9) and (10) from Eq. (5). We can also express Eq. (7) as

$\displaystyle\Omega_{\rm m}+\Omega_{\rm curv}$ $\displaystyle=$
$\displaystyle 1$ (11)

Here $\Omega_{\rm m}$, and $\Omega_{\rm curv}$ are redefined density
parameters by the following equations;

$\displaystyle\Omega_{\rm m}\equiv\frac{\kappa^{2}\rho_{\rm
m}}{3H^{2}}=\frac{\kappa^{2}\tilde{\rho}_{\rm m}}{3FH^{2}}\,,\quad\Omega_{\rm
curv}\equiv\frac{\kappa^{2}\rho_{\rm curv}}{3H^{2}}$ (12)

Taking divergence in both sides of (4) and applying Bianchi’s identity yields
the conservation equation for the total stress-energy tensor as,

$\displaystyle\nabla^{\mu}T_{\mu\nu}^{\rm tot}$ $\displaystyle=$
$\displaystyle\nabla^{\mu}\left(T_{\mu\nu}^{(\rm M)}+T_{\mu\nu}^{(\rm
curv)}\right)=0$ (13)

In FLRW spacetime, (13) becomes a continuity equation for the inclusive fluid
with energy density. ${\rho}_{\rm tot}\equiv\rho_{\rm m}+\rho_{\rm curv}$ and
pressure $P_{\rm tot}\equiv P_{\rm curv}$ as

$\displaystyle\dot{\rho}_{\rm tot}+3H{\rho}_{\rm tot}(1+\omega_{\rm
tot})=0\,,$ (14)

Here, $\omega_{\rm tot}$ denotes the grand equation of state (EoS) parameter
for the combined (matter and curvature) fluid. While the total stress-energy
tensor $T^{(\rm tot)}{\mu\nu}$ remains conserved, individual components
$T^{(\rm curv)}_{\mu\nu}$ and $T^{(M)}_{\mu\nu}$ may not be individually
conserved. This allows for potential interactions between the matter and
curvature sectors.

With these considerations, we express the conservation Eq. (13) as
$\nabla^{\mu}T_{\mu\nu}^{(\rm m)}=-\nabla^{\mu}T_{\mu\nu}^{(\rm
curv)}\equiv-\mathcal{J}_{\nu}\neq 0$. In the context of FLRW space-time, non-
conserving continuity equation of curvature-matter sectors can be written as,

$\displaystyle\dot{\rho}_{\rm curv}+3H\left(\rho_{\rm curv}+P_{\rm
curv}\right)$ $\displaystyle=$ $\displaystyle\mathcal{Q}$ (15)
$\displaystyle\dot{\rho}_{\rm m}+3H\rho_{\rm m}$ $\displaystyle=$
$\displaystyle-\mathcal{Q}$ (16)

Where source term $\mathcal{Q}$ represents a time-varying function indicating
the instantaneous rate of energy exchange between the curvature and matter
sectors. Prior research has explored related avenues [46, 45], adopting
interaction forms based solely on matter components or the multiplicative
nature of both curvature and matter sectors. However, our study uniquely
delves into both forms of interactions, $\mathcal{Q}=\alpha H\tilde{\rho}_{\rm
m}\left(\frac{\kappa^{2}}{3H^{2}}\rho_{\rm curv}+1\right)$, employing a
comprehensive phase-space analysis method. We restrict our study only to the
de-sitter phase, with matter components minimally coupled to curvature sectors
under a flat FLRW metric background. Here, $\alpha$ is the coupling strength
of the interactions. This specific type of interaction incorporates degrees of
freedom from both sectors and also ensures equality between matter and
curvature in terms of critical density (where, $\frac{\Omega_{\rm
m}}{\Omega_{\rm curv}}=r_{\rm mc}\sim 1$), depending on the strength of the
coupling parameter ($\alpha$).

Imposing $\rho_{\rm curv}>0$ assigns the concept of ‘energy density’ to it,
placing constraints on $f(R)$ models to satisfy $6H\dot{R}F^{\prime}<RF-f$ as
seen in Eq. (9). This constraint also impacts the coupling parameter $\alpha$
through its connection with source term $Q$ in Eq. (15). Under this condition,
each modified density parameter $\Omega_{\rm m}$, and $\Omega_{\rm curv}$
stays positive and collectively follows the constraint in Eq. (11). This
allows us to explore the parameter space for curvature-matter interaction by
varying $\Omega_{\rm m}$ within the range $0\leqslant\Omega_{\rm m}\leqslant
1$, leading to constraints on the model parameters. Additional constraints
arise from the non-phantom and accelerating nature of the grand equation of
state (EoS) parameter, along with cosmographic parameters (deceleration and
jerk) and the critical energy density ratio of matter-to-curvature components.
Further details on these constraints will be provided in the following
section.

##  3 Dynamical analysis of interacting matter-curvature scenario in two types
of $f(R)$ models

Dynamical analysis is employed to investigate curvature-matter interaction in
cosmologically viable $f(R)$ gravity models. We can express autonomous
equations in terms of dynamically relevant variables, laying the groundwork
for applying dynamical analysis techniques in our context. Fundamental
dynamical variables are,

$\displaystyle X_{1}=-\frac{\dot{F}}{HF},\quad X_{2}=-\frac{f}{6FH^{2}},\quad
X_{3}=\frac{R}{6H^{2}}$ (17)

Introducing the dimensionless parameter $N=\ln a$ and using the defined
dynamical variables ($X_{i}$’s), Eqs. (9), (15), and (16) become a set of four
autonomous equations, forming a 3-D dynamical system.

$\displaystyle\frac{dX_{1}}{dN}$ $\displaystyle=$
$\displaystyle-1-X_{1}X_{3}-X_{3}-3X_{2}+X_{1}^{2}$ (18)
$\displaystyle\qquad+\
\alpha\left(1-X_{1}-X_{2}-X_{3}\right)\left(1+X_{1}+X_{2}+X_{3}\right)$
$\displaystyle\frac{dX_{2}}{dN}$ $\displaystyle=$
$\displaystyle\frac{X_{1}X_{3}}{m}-X_{2}\left(2X_{3}-4-X_{1}\right)$ (19)
$\displaystyle\frac{dX_{3}}{dN}$ $\displaystyle=$
$\displaystyle-\frac{X_{1}X_{3}}{m}-2X_{3}\left(X_{3}-2\right)\,,$ (20)

Here, $m$ and $r$ parameters can be defined as,

$\displaystyle m$ $\displaystyle\equiv$ $\displaystyle\frac{d\ln F}{d\ln
R}=\frac{RF^{\prime}}{F}$ $\displaystyle r$ $\displaystyle\equiv$
$\displaystyle-\frac{d\ln f}{d\ln R}=-\frac{RF}{f}=\frac{X_{3}}{X_{2}}\,.$
(21)

The dimensionless parameters $m$ and $r$ depend on $R$ through $f(R)$,
enabling us to express $m$ as a function of $r$, denoted as $m=m(r)$. Each
specific functional relation $m=m(r)$ corresponds to a unique class of $f(R)$
models. With $\alpha$ sets to zero, the equations reduce to autonomous
equations for a scenario without interactions between curvature and matter
sectors.

Eqs. (11) and (12) impose constraints on the dynamical variables in Eq. (17)
as,

$\displaystyle\Omega_{\rm m}=1-X_{1}-X_{2}-X_{3}$ (22)

We can further express $\Omega_{\rm curv}$ as,

$\displaystyle\Omega_{\rm curv}=1-\Omega_{\rm m}=X_{1}+X_{2}+X_{3}\,.$ (23)

Furthermore, the grand EoS parameter ($\omega_{\rm tot}$) can be written as,

$\displaystyle\omega_{\rm tot}$ $\displaystyle=$
$\displaystyle-1-\frac{2\dot{H}}{3H^{2}}=-\frac{1}{3}(2X_{3}-1)$ (24)

The ratio of critical energy density of matter to curvature, denoted as
$r_{mc}\equiv\Omega_{\rm m}/\Omega_{\rm curv}$, along with two cosmographic
parameters - deceleration $(q\equiv-a\ddot{a}/\dot{a}^{2})$ and jerk
$(j\equiv\dddot{a}/aH^{3})$ \- serve as valuable indicators for assessing
various aspects of cosmic dynamics, mainly in late-time scenario. In this
matter-curvature interaction scenario, the expressions for the three
parameters ($r_{mc},q,j$) can be expressed in terms of the dynamical variables
($X_{i}$’s) as,

$\displaystyle r_{mc}$ $\displaystyle\equiv$ $\displaystyle\frac{\Omega_{\rm
m}}{\Omega_{\rm curv}}=\frac{1-X_{1}-X_{2}-X_{3}}{X_{1}+X_{2}+X_{3}}$ (25)
$\displaystyle q$ $\displaystyle\equiv$
$\displaystyle-\frac{a\ddot{a}}{\dot{a}^{2}}=-1-\frac{\dot{H}}{H^{2}}=1-X_{3}$
(26) $\displaystyle j$ $\displaystyle\equiv$
$\displaystyle\frac{\dddot{a}}{aH^{3}}=-\frac{X_{1}X_{3}}{m}+2(1-X_{3})^{2}+(1-X_{3})-2X_{3}(X_{3}-2)$
(27)

The fixed points in the 3-D dynamical system, where $dX_{i}/dN=0$ ($i=1,2,3$),
are crucial for understanding cosmic evolution driven by curvature-matter
interactions. Linear stability analysis, involving a first-order Taylor
expansion, is used to assess stability around these points. The Jacobian
matrix, $\rm J=||{\partial f_{i}}/{\partial x_{j}}||$ aids in determining
stability, with asymptotic stability (unstable) indicated when all eigenvalues
have negative (positive) real parts. A saddle point is identified if any pair
of eigenvalues has a relative opposite sign in their real parts. If any
eigenvalue approaches zero, center manifold theory is required for a more
insightful exploration of these critical points characteristics.

This entire framework of dynamical analysis is based on $f(R)$-driven dark
energy models, considering interactions between curvature and matter sectors.
We selected specific models capable of inducing cosmic acceleration while
remaining cosmologically viable. In the metric formalism, any viable $f(R)$
function must meet stringent conditions discussed in [50, 51]. Below, we list
and provide the rationale for these conditions. We examine the autonomous
system, focusing on two viable scenarios: The first involves a $f(R)$ yielding
a constant $m$, while the second pertains to a $f(R)$ where $m$ is a function
of $r$, denoted as $m(r)$. These scenarios are labeled as (I) and (II), with
their specifications discussed below.

  * (I)

In the modified gravity model, $f(R)=(R^{b}-\Lambda)^{c}$ with $(c\geqslant
1,bc\approx 1)$, a constant $m$ scenario arises. Explored without considering
matter-curvature interactions in [44, 49], this model extends the
$\Lambda$-CDM model, termed ‘generalized $\Lambda$-CDM model’. Constrained to
$c\geqslant 1$ and $bc\approx 1$, it converges to the $\Lambda$-CDM model with
$m=0$, ensuring viable cosmological evolution. The parameters $m$ and $r$ are
given by $m=\frac{(bc-1)R^{b}-b\Lambda+\Lambda}{R^{b}-\Lambda}$,
$r=-\frac{bcR^{b}}{R^{b}-\Lambda}$, and $m$ can be expressed as
$m(r)=\left(\frac{1-c}{c}\right)r+b-1$. A dynamical analysis (see [44, 49] for
details) indicates stable points where $r=-1-m$, resulting in $m=-1+bc$, a
constant value.

  * (II)

Power law model considers $m$ as a specific function of $r$ given by
$m(r)=\frac{n(1+r)}{r}$, attainable with the power law form $f(R)=R-\gamma
R^{n}$ ($\gamma>0$, $0<n<1$). This scenario accurately describes cosmic
evolution within the framework of non-interacting curvature-matter scenarios,
meeting the conditions for cosmological viability within the specified range
of $\gamma$ and $n$, as extensively discussed in [44, 52]. The corresponding
$f(R)$ form yields $m=\frac{\gamma(n-1)nR^{n}}{R-\gamma nR^{n}}$ and
$r=\frac{R-\gamma nR^{n}}{R-\gamma R^{n}}$, with the elimination of $\gamma$
leading to the form $m(r)=\frac{n(1+r)}{r}$.

In our study, the selection of the two cases: (I) $m=-1+bc=$ constant and (II)
$m=m(r)=\frac{n(1+r)}{r}$ aims to investigate implications of the mentioned
$f(R)$-models in the presence of matter-curvature interactions at the de-
sitter phase.

From the formulation of the dynamical variables in Eq. (17), it is evident
that both $X_{1}$ and $X_{2}$ are contingent on the chosen $f(R)$ models. In
contrast, $X_{3}$ remains independent of the $f(R)$ selection and assumes a
pivotal role in defining the grand equation of state parameter. Given our
primary objective of investigating the dynamics of curvature-matter
interaction in viable $f(R)$ dark energy models specifically during the de-
sitter phase, we have kept the variable $X_{3}$ fixed. This action reduces the
phase-space dimension from 3-D to 2-D, focusing exclusively on a constant
$X_{3}$ plane. The grand Equation of State (EoS) parameter at the de-sitter
phase, as indicated in Eq. (24), yields a value of $-1$, which consequently
fixes the $X_{3}$ variables at 2. Simultaneously, the deceleration parameter
($q$) in Eq. (26) is also negative unity at this point. Subsequently, we have
established a 2-D phase space by varying $X_{1}$ and $X_{2}$ for specific
benchmark values of $X_{3}$ (set to 2), where the system behaves akin to a de-
sitter solution.

##  4 Phase-space analysis of matter-curvature coupling in modified $f(R)$
gravity models

In both models, we identify fixed points of the autonomous system (see in sec.
[3]) governing curvature-matter interactions in $f(R)$ gravity at de-sitter
phase. Obtained by solving $dX_{k}/dN=0$ ($k=1,2$), using $m=-1+bc$ for model
I and $m=\frac{n(1+r)}{r}=\frac{n(X_{2}+2)}{2}$ for model II (where, $X_{3}=2$
at de-sitter phase). Coupling parameter $\alpha$ with $(b,c)$ for model I and
$(n)$ for model B are denoted as the free parameters of this interacting
picture. Critical density parameters ($\Omega_{\rm m}$, $\Omega_{\rm curv}$)
at fixed points $(X_{i})$ are evaluated using Eqs. (22), (23). Positive
modified density parameters, subject to the constraint in Eq. (11), determine
cosmological significance. Describing fixed points and their stability
references relevant constraints within the model parameter space. Value of
density parameters at each fixed point provides insights into characterizing
the associated de-sitter phase. In the following subsections, we have
thoroughly examined interacting curvature-matter framework in the context of
two types of $f(R)$ models at de-sitter phase through linear stability
analysis theory. By fixing $X_{3}=2$, we’ve analyzed the fixed points behavior
with variables $X_{1}$ and $X_{2}$ and illustrated their dynamics in phase
space.

###  4.1 Analysis of generalized $\Lambda$-CDM model:

We have initiate our investigation by employing the generalized $\Lambda$-CDM
model to scrutinize the dynamic behavior of the coupled system. By exploiting
the set of autonomous eqautions (18), (19) and also setting $X_{3}=2$, a total
of four critical points have been found. We have tabulated the critical points
in tab.[1].

Critical points (at $X_{3}=2$)
---
Critical Points | $X_{1}$ | $X_{2}$
$P_{1}$ | $0$ | $-1$
$P_{2}$ | $0$ | $-\frac{3(\alpha+1)}{\alpha}$
$P_{3\mp}$ | $\frac{-8\alpha-4\alpha\text{bc}^{2}-2\text{bc}^{2}+12\alpha\text{bc}+4\text{bc}-2}{(2(-1+2bc-bc^{2}+\alpha-2bc\alpha+bc^{2}\alpha))}$ | $-\frac{2}{\text{bc}-1}$
| $\mp\frac{1}{2\left(\alpha+\alpha\text{bc}^{2}-\text{bc}^{2}-2\alpha\text{bc}+2\text{bc}-1\right)}\Big{[}-4\left(\alpha+\alpha\text{bc}^{2}-\text{bc}^{2}-2\alpha\text{bc}+2\text{bc}-1\right)$ |
| $\left(15\alpha+3\alpha\text{bc}^{2}+3\text{bc}^{2}-14\alpha\text{bc}-12\text{bc}+9\right)$ |
| $+\left(8\alpha+4\alpha\text{bc}^{2}+2\text{bc}^{2}-12\alpha\text{bc}-4\text{bc}+2\right)^{2}\Big{]}^{1/2}$ |
Table 1: Critical points of generalized $\Lambda$-CDM model at de-sitter phase

Further details about the attributes of these critical points are discussed
below,

  * •

$P_{1}$ Point: This critical point remains model and coupling-parameter
independent. The critical matter density ($\Omega_{m}$) at this point is
consistently zero, making it a curvature-dominant critical point.
Additionally, the ratio of critical matter to curvature density at this
specific point is always zero. Although it exhibits a positive unit jerk in
late times, the positive eigenvalues render this point consistently unstable.

  * •

$P_{2}$ Point: The characteristics of this critical point solely depend on the
coupling parameter. The critical matter density ($\Omega_{m}$) and the ratio
of critical matter to curvature density ($r_{\rm mc}$) at this point can be
expressed in terms of the coupling parameter $\alpha$ as
$\frac{3(\alpha+1)}{\alpha}-1$ and
$\frac{\frac{3(\alpha+1)}{\alpha}-1}{2-\frac{3(\alpha+1)}{\alpha}}$,
respectively. Similar to the previous critical point, this point consistently
exhibits a positive jerk of unity. The conditions for stability, critical
matter density, and a positive ratio of critical matter to curvature density
impose a common constraint in the parameter space spanned by $bc$ and
$\alpha$, as illustrated in fig. [1].

Figure 1: Parameter space of generalized $\Lambda$-CDM model at $P_{2}$ in de-
sitter phase

Figure 2: Phase space trajectories of generalized $\Lambda$-CDM model at
different benchmark points. Critical points $P_{1}$ and $P_{2}$ are identified
as repeller and attractor, respectively. The magenta region signifies positive
jerk, and the yellow band represents $0\leq\Omega_{m}\leq 1$ with $r_{\rm
mc}\geq 0$. The central green line within the yellow band indicates the point
where $r_{\rm mc}=1$.

From the shared permissible region defined by stability conditions ($r_{\rm
mc}\geq 0$, $0\leq\Omega_{m}\leq 1$), we selected three specific values for
the parameters $bc$ and $\alpha$ to generate a phase-space plot involving
variables $X_{1}$ and $X_{2}$. As shown in fig.[2], we illustrate the behavior
of both critical points, $P_{1}$ (independent of model parameters) and $P_{2}$
(dependent on model parameters). Across all three plots, the phase-space
behavior of $P_{1}$ exhibits a repeller-type nature, with trajectories being
repelled from $P_{1}$ and attracted towards $P_{2}$, establishing $P_{2}$ as a
global attractor point. The emergence of $P_{2}$ is attributed solely to the
presence of interaction. In the phase plot, the $X_{2}$ coordinate of $P_{2}$
is observed to be dependent on the coupling parameter ($\alpha$), shifting
towards the negative direction for large negative values of the coupling
parameter. The yellow band in the phase space indicates critical matter
density between 0 and 1, with a positive jerk depicted in the magenta-marked
region. The central green line within the yellow band signifies the epoch when
critical matter and curvature density are the same. At this point, the $r_{\rm
mc}$ parameter attains a value of one. In the context of a significantly
negative coupling parameter ($\alpha\sim-3$), the attractor point $P_{2}$ is
situated at the border of the yellow band. The analysis of the phase space at
the de-sitter point leads to the conclusion that curvature-matter interaction
results in a critical point that satisfies all conditions and yields a stable
attractor critical point for the system.

  * •

$P_{3\mp}$ Points: Only these points can be expressed regarding both model and
coupling parameters. However, they fail to meet the specified physicality
conditions and yield unstable critical points through linear stability
analysis. The energy density and cosmographic parameters result in non-
physical values at these points within the accessible region of model
parameters identified from the $P_{2}$ point. Due to this restriction, these
points cannot be included in the phase space diagram in fig.[2].

###  4.2 Analysis of power-law model:

We commence our exploration by utilizing the power-law model to examine the
dynamic characteristics of the autonomous system at de-sitter phase. By
analyzing the set of autonomous equations (18) and (20) while fixing $X_{3}$
to 2, a total of four fixed points have been identified, which are given in
tab.[2].

Critical points (at $X_{3}=2$)
---
Critical Points | $X_{1}$ | $X_{2}$
$P_{1}$ | $0$ | $-1$
$P_{2}$ | $0$ | $-\frac{3(\alpha+1)}{\alpha}$
$P_{3\mp}$ | $2(\alpha n-n)\Big{[}-2\alpha\left(-\sqrt{n^{2}-4n}-n\right)-4\alpha n-2n$ | $\frac{-\sqrt{n^{2}-4n}-n}{n}$
| $\mp 2\sqrt{n^{2}-3n\sqrt{n^{2}-4n}-4\alpha n+3\alpha n^{2}+\alpha^{2}\left(n^{2}-4n\right)-\alpha n\sqrt{n^{2}-4n}+4\alpha^{2}n}\Big{]}$ |
$P_{4\mp}$ | $2(\alpha n-n)\Big{[}-2\alpha\left(-\sqrt{n^{2}-4n}-n\right)-4\alpha n-2n$ | $\frac{-\sqrt{n^{2}-4n}-n}{n}$
| $\mp 2\sqrt{n^{2}-3n\sqrt{n^{2}-4n}-4\alpha n+3\alpha n^{2}+\alpha^{2}\left(n^{2}-4n\right)+\alpha n\sqrt{n^{2}-4n}+4\alpha^{2}n}\Big{]}$ |
Table 2: Critical Points of power-law model at de-sitter phase

Like the previous model, here we also employ linear stability theory to
analyze the dynamical behavior of power-law model. This involves determining
the Jacobian matrix and subsequently evaluating the stability of critical
points. Further details regarding the properties of these critical points are
discussed below,

  * •

$P_{1}$ Point: This specific critical point remains independent of model
parameters. The critical matter density is consistently zero at this juncture,
indicating the dominance of curvature-driven dark energy. The stability
condition of eigenvalues, $r_{\rm mc}\geq 0$ and $0\leq\Omega_{m}\leq 1$
establish an allowed region, which is expressed in terms of the coupling
parameter $\alpha$ and model parameter $n$. In fig. [3], we illustrate the
permissible region of this power-law model during the de-sitter phase
corresponding to the $P_{1}$ point. The parameter space nature is divided into
two regions based on the model parameter $n$, namely (i) $0<n<1$ and (ii)
$n>1$. In the first region, the accessible range of coupling strength
($\alpha$) initiates from a numerical value of approximately $5.5$, while in
the second case, the allowed range of coupling parameter ($\alpha$) starts
from a small negative value ($\sim-0.5$). The coupling parameter $\alpha$
plays a crucial role in constraining $n$ in terms of stability and other
physical conditions, as mentioned earlier.

Figure 3: Parameter space of $n$ vs. $\alpha$ of power-law model at $P_{1}$ in
de-sitter phase

  * •

$P_{2}$ Point: While this critical point is expressed as a function of the
coupling strength $\alpha$, we are unable to identify any stable critical
point. Due to its tendency to exhibit unstable behavior, we refrain from
presenting detailed information on energy density and other parameters for
this specific point.

We have selected two benchmark values for model parameters ($n,\alpha$) from
each region in the parameter space shown in fig. [3] and examined the 2-D
phase-space behavior, as depicted in fig. [4]. In the first region, we have
chosen lower values for both coupling strength ($\alpha$) and model parameter
($n$) to observe the phase-space trajectories in the $X_{1}$-$X_{2}$ plane. In
the power-law model, we have found that the value of $m=n(X_{2}+2)/2$, and
when $X_{2}=-X_{3}=-2$, $m$ becomes zero, rendering dynamical Eq. (19)
undefined at the de-sitter solution. For this issue, we have identified a
disjoint phase-space behavior marked by a black dotted line, where the
dynamical variable $X_{2}=-X_{3}=-2$. We have chosen benchmark values from the
accessible range of $n,\alpha$ for the model parameter-independent critical
points $P_{1}$, which also produced values for point $P_{2}$, plotted in the
same phase space. All trajectories are repelled from the point $P_{2}$ and
attracted towards $P_{1}$ making it into the attractor. Higher positive values
of both parameters ($n,\alpha$) influence the movement of the critical point
$P_{2}$ along the negative $X_{2}$-axis. Similar to the previous model, we
marked a magenta region for a positive jerk parameter and a yellow band for
the region where critical matter density is between 0 and 1, as well as the
ratio of critical matter to curvature density being greater than zero. The
middle green line in the yellow band indicates the time when critical matter
and curvature density are equal. At this juncture, the $r_{\rm mc}$ parameter
reaches a value of one. In all plots, the attractor point $P_{1}$ satisfies
all conditions and lies inside the common intersecting region of the magenta
and yellow bands. Similar behavior has been observed in the other two cases
where $n>1$.

Figure 4: Phase space trajectories of power-law model at different benchmark
points. Critical points $P_{1}$ and $P_{2}$ are designated as a repeller and
an attractor, respectively. The magenta region signifies positive jerk, and
the yellow band represents $0\leq\Omega_{m}\leq 1$ with $r_{\rm mc}\geq 0$.
The central green line within the yellow band indicates the point where
$r_{\rm mc}=1$.

  * •

$P_{3\mp}$ and $P_{4\mp}$ Points: While both of these critical points can be
defined in terms of model and coupling parameters, they both fall short of
offering stable solutions during the de-sitter phase. Additionally, they
exhibit non-physical behavior in terms of critical matter density at the
specific benchmark point chosen to analyze the phase-space trajectories of the
2-D system. Due to this non-physical characteristic observed in the density
parameters at points $P_{3\mp}$ and $P_{4\mp}$ for the benchmark values of $n$
and $\alpha$, which have been chosen from the accessible range of $P_{1}$. We
omitted these points from the phase space plot in fig.[4].

##  5 Conclusion

In this article, we have investigated the potential incorporation of
interactions between curvature and matter within the context of viable $f(R)$
dark-energy models. Our approach involves introducing a minimally coupled
interacting term at the action level, allowing us to analyze the overall
dynamics using a flat FLRW metric as the background. The selected form of
interaction is designed to leave the total conserving continuity equation for
the curvature and matter sectors unaffected. However, at the individual sector
level, the presence of a source term ($\mathcal{Q}$) becomes crucial in
influencing the dynamics. The chosen interaction term ($\mathcal{Q}$)
effectively combines the influences of both matter and curvature sectors,
determining the rate of energy transfer between them. This transfer rate
reaches either a maximum or minimum based on the nature of the coupling
strength ($\alpha$).

This study’s crucial aspect involves simplifying the 3-D dynamical framework
to a 2-D framework by keeping one of the $f(R)$ model-independent dynamical
variables constant. Subsequently, it examines the dynamics of the curvature-
matter interaction during the de-sitter phase using the linear stability
approach. The inclusion of the interaction term $\mathcal{Q}$ results in
modifying the autonomous equations within the dynamical system compared to
non-interacting scenarios. We identified critical points in the dynamical
system and assessed their stability during the de-sitter phase, separately
considering two viable $f(R)$ models, namely, the generalized $\Lambda$-CDM
model and the power-law model. The study thoroughly investigates how
interaction parameters impact the stability characteristics of various fixed
points within the dynamical system.

The introduction of a novel curvature-matter coupling in the context of $f(R)$
gravity results in a modification of critical points during the de-sitter
phase. Common acceptable regions, satisfying stability, critical matter
density, and the ratio of matter to curvature density conditions, have been
identified for both models based on coupling strength and model parameters.
Different benchmark points within these regions have been selected to
illustrate the behavior of critical points in 2-D phase-space plots in the
presence of interaction. In each model, an attractor-type critical point has
been identified that met all specified criteria. The nature of the attractor
point has been categorized based on the $f(R)$ model. In the generalized
$\Lambda$-CDM model, the attractor point depends on both model and coupling
parameters, while in the power-law model, it is entirely independent of these
parameters. For the stable attractor point in the generalized $\Lambda$-CDM
model, the allowed coupling range is $\alpha<-2.15$, with the model parameter
$bc$ ranging between 1.85 and 2.25, as seen in fig.[1]. In contrast, for the
power-law model, the acceptable coupling range is $\alpha>-1$ and the model
parameter $n$ falls between $0.35-4$, as depicted in fig.[3]. All these
parameter ranges are accompanied by the de-sitter phase.

Since our study primarily revolves around investigating matter-curvature
interaction in viable $f(R)$ dark energy models during the de-sitter phase,
this specific phase constrains the grand equation of state parameter
($\omega_{\rm tot.}$) to be equal to -1. This constraint, in turn, imposes
limitations on the selection of the dynamical variable $X_{3}$ at 2. With this
consideration, we have analyzed the 2-D phase space defined by variables
$X_{1}$ and $X_{2}$ in the presence of interacting curvature-matter
interaction for both the generalized $\Lambda$-CDM model and the power-law
model. In the phase portrait of each model, we have included all critical
points that exhibit physical behavior concerning critical matter density at
selected benchmark values of model parameters ($bc$ or $n$) and coupling
strength ($\alpha$). In the generalized $\Lambda$-CDM scenario, we have
identified a model parameter-dependent critical point $P_{2}$ as an attractor,
while in the power-law scenario, a model parameter-independent critical point
$P_{1}$ has been identified as an attractor. Both attractors satisfy all
physicality criteria and are situated within the shared region of magenta and
yellow in the phase plot. The cosmological coincidence between the curvature-
driven dark energy and dark matter is highlighted by the green line, where the
energy densities of both sectors are equal. The attractor points are also near
this line, indicating that the cosmological coincidence between dark energy
and dark matter is achieved by introducing this specific form of interaction
during the de-sitter phase.

In conclusion, this study advances our comprehension of the intricate dynamics
of interacting curvature-matter scenarios within the framework of modified
$f(R)$ gravity theories. The identification of stable attractors during the
de-sitter phase, nuanced stability features, and phase trajectory patterns
contribute to a comprehensive understanding of the complex behavior exhibited
by the system. The insights presented not only deepen our knowledge of the
dynamics in these systems but also provide new opportunities for exploration
in cosmology, offering valuable perspectives into the fundamental dynamics of
the universe.

#### Acknowledgement

Author would like to convey appreciation to the Department of Physics at
RKMVERI for their warm hospitality during the author’s tenure as a visiting
fellow. Special thanks are extended to Abhijit Bandyopadhyay and Rahul Roy for
their enriching and insightful discussions.

## References

  * [1] A. G. Riess et al. [Supernova Search Team], Astron. J. 116, 1009 (1998) doi:10.1086/300499 [astro-ph/9805201].
  * [2] S. Perlmutter et al. [Supernova Cosmology Project Collaboration], Astrophys. J. 517, 565 (1999) doi:10.1086/307221 [astro-ph/9812133].
  * [3] D. N. Spergel et al. [WMAP], Astrophys. J. Suppl. 148 (2003), 175-194 doi:10.1086/377226 [arXiv:astro-ph/0302209 [astro-ph]]
  * [4] G. Hinshaw et al. [WMAP], Astrophys. J. Suppl. 180 (2009), 225-245 doi:10.1088/0067-0049/180/2/225 [arXiv:0803.0732 [astro-ph]].
  * [5] D. J. Eisenstein et al. [SDSS], Astrophys. J. 633 (2005), 560-574 doi:10.1086/466512 [arXiv:astro-ph/0501171 [astro-ph]]
  * [6] I. Zlatev, L. M. Wang and P. J. Steinhardt, Phys. Rev. Lett. 82 (1999), 896-899 doi:10.1103/PhysRevLett.82.896 [arXiv:astro-ph/9807002 [astro-ph]]
  * [7] J. Martin, Comptes Rendus Physique 13 (2012), 566-665 doi:10.1016/j.crhy.2012.04.008 [arXiv:1205.3365 [astro-ph.CO]].
  * [8] R. D. Peccei, J. Sola and C. Wetterich, Phys. Lett. B 195 (1987), 183-190 doi:10.1016/0370-2693(87)91191-9
  * [9] L. H. Ford, Phys. Rev. D 35 (1987), 2339 doi:10.1103/PhysRevD.35.2339
  * [10] P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75 (2003), 559-606 doi:10.1103/RevModPhys.75.559 [arXiv:astro-ph/0207347 [astro-ph]].
  * [11] T. Nishioka and Y. Fujii, Phys. Rev. D 45 (1992), 2140-2143 doi:10.1103/PhysRevD.45.2140
  * [12] P. G. Ferreira and M. Joyce, Phys. Rev. Lett. 79 (1997), 4740-4743 doi:10.1103/PhysRevLett.79.4740 [arXiv:astro-ph/9707286 [astro-ph]].
  * [13] P. G. Ferreira and M. Joyce, Phys. Rev. D 58 (1998), 023503 doi:10.1103/PhysRevD.58.023503 [arXiv:astro-ph/9711102 [astro-ph]]
  * [14] R. R. Caldwell, R. Dave and P. J. Steinhardt, Phys. Rev. Lett. 80 (1998), 1582-1585 doi:10.1103/PhysRevLett.80.1582 [arXiv:astro-ph/9708069 [astro-ph]].
  * [15] S. M. Carroll, Phys. Rev. Lett. 81 (1998), 3067-3070 doi:10.1103/PhysRevLett.81.3067 [arXiv:astro-ph/9806099 [astro-ph]]
  * [16] E. J. Copeland, A. R. Liddle and D. Wands, Phys. Rev. D 57 (1998), 4686-4690 doi:10.1103/PhysRevD.57.4686 [arXiv:gr-qc/9711068 [gr-qc]].
  * [17] W. Fang, H. Tu, Y. Li, J. Huang and C. Shu, Phys. Rev. D 89 (2014) no.12, 123514 doi:10.1103/PhysRevD.89.123514 [arXiv:1406.0128 [gr-qc]].
  * [18] C. Armendariz-Picon, T. Damour and V. F. Mukhanov, Phys. Lett. B 458 (1999), 209-218 doi:10.1016/S0370-2693(99)00603-6 [arXiv:hep-th/9904075 [hep-th]]
  * [19] C. Armendariz-Picon, V. F. Mukhanov and P. J. Steinhardt, Phys. Rev. D 63 (2001), 103510 doi:10.1103/PhysRevD.63.103510 [arXiv:astro-ph/0006373 [astro-ph]]
  * [20] C. Armendariz-Picon, V. F. Mukhanov and P. J. Steinhardt, Phys. Rev. Lett. 85 (2000), 4438-4441 doi:10.1103/PhysRevLett.85.4438 [arXiv:astro-ph/0004134 [astro-ph]]
  * [21] C. Armendariz-Picon and E. A. Lim, JCAP 08 (2005), 007 doi:10.1088/1475-7516/2005/08/007 [arXiv:astro-ph/0505207 [astro-ph]]
  * [22] T. Chiba, T. Okabe and M. Yamaguchi, Phys. Rev. D 62 (2000), 023511 doi:10.1103/PhysRevD.62.023511 [arXiv:astro-ph/9912463 [astro-ph]]
  * [23] N. Arkani-Hamed, H. C. Cheng, M. A. Luty and S. Mukohyama, JHEP 05 (2004), 074 doi:10.1088/1126-6708/2004/05/074 [arXiv:hep-th/0312099 [hep-th]]
  * [24] R. R. Caldwell, Phys. Lett. B 545 (2002), 23-29 doi:10.1016/S0370-2693(02)02589-3 [arXiv:astro-ph/9908168 [astro-ph]]
  * [25] A. Bandyopadhyay and A. Chatterjee, Mod. Phys. Lett. A 34 (2019) no.27, 1950219 doi:10.1142/S0217732319502195 [arXiv:1709.04334 [gr-qc]]
  * [26] A. Bandyopadhyay and A. Chatterjee, Eur. Phys. J. Plus 134 (2019) no.4, 174 doi:10.1140/epjp/i2019-12587-0 [arXiv:1808.05259 [gr-qc]]
  * [27] A. Bandyopadhyay and A. Chatterjee, Eur. Phys. J. Plus 135 (2020) no.2, 181 doi:10.1140/epjp/s13360-020-00161-w [arXiv:1902.04315 [gr-qc]]
  * [28] A. Bandyopadhyay and A. Chatterjee, Res. Astron. Astrophys. 21 (2021) no.1, 002 doi:10.1088/1674-4527/21/1/2 [arXiv:1910.10423 [gr-qc]]
  * [29] A. Chatterjee, B. Jana and A. Bandyopadhyay, Eur. Phys. J. Plus 137 (2022) no.11, 1271 doi:10.1140/epjp/s13360-022-03476-y [arXiv:2207.00888 [gr-qc]]
  * [30] S. Capozziello, Int. J. Mod. Phys. D 11 (2002), 483-492 doi:10.1142/S0218271802002025 [arXiv:gr-qc/0201033 [gr-qc]]
  * [31] S. Capozziello, V. F. Cardone, S. Carloni and A. Troisi, Int. J. Mod. Phys. D 12 (2003), 1969-1982 doi:10.1142/S0218271803004407 [arXiv:astro-ph/0307018 [astro-ph]]
  * [32] S. Nojiri and S. D. Odintsov, Phys. Rev. D 68 (2003), 123512 doi:10.1103/PhysRevD.68.123512 [arXiv:hep-th/0307288 [hep-th]]
  * [33] S. Nojiri and S. D. Odintsov, Phys. Rept. 505 (2011), 59-144 doi:10.1016/j.physrep.2011.04.001 [arXiv:1011.0544 [gr-qc]]
  * [34] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Phys. Rept. 692 (2017), 1-104 doi:10.1016/j.physrep.2017.06.001 [arXiv:1705.11098 [gr-qc]]
  * [35] B. Jana, A. Chatterjee, K. Ravi and A. Bandyopadhyay, Class. Quant. Grav. 40 (2023) no.19, 195023 doi:10.1088/1361-6382/acf554 [arXiv:2303.06961 [gr-qc]]
  * [36] K. Ravi, A. Chatterjee, B. Jana and A. Bandyopadhyay, MNRAS stad3705, 2023 doi:10.1093/mnras/stad3705 [arXiv:2306.12585 [astro-ph.CO]]
  * [37] A. Chatterjee, S. Hussain and K. Bhattacharya, Phys. Rev. D 104 (2021) no.10, 2021 doi:10.1103/PhysRevD.104.103505 [arXiv:2105.00361 [gr-qc]]
  * [38] A. Chatterjee, A. Bandyopadhyay and B. Jana, Eur. Phys. J. Plus 137 (2022) no.4, 518 doi:10.1140/epjp/s13360-022-02747-y [arXiv:2108.12186 [gr-qc]]
  * [39] S. Hussain, A. Chatterjee and K. Bhattacharya, Universe 9 (2023) no.2, 65 doi:10.3390/universe9020065 [arXiv:2203.10607 [gr-qc]]
  * [40] K. Bhattacharya, A. Chatterjee and S. Hussain, Eur. Phys. J. C 83 (2023) no.6, 488 doi:10.1140/epjc/s10052-023-11666-w [arXiv:2206.12398 [gr-qc]]
  * [41] S. Hussain, A. Chatterjee and K. Bhattacharya, Phys. Rev. D 108 (2023) no.10, 103502 doi:10.1103/PhysRevD.108.103502 [arXiv:2305.19062 [gr-qc]].
  * [42] A. De Felice and S. Tsujikawa, Living Rev. Rel. 13 (2010), 3 doi:10.12942/lrr-2010-3 [arXiv:1002.4928 [gr-qc]]
  * [43] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82 (2010), 451-497 doi:10.1103/RevModPhys.82.451 [arXiv:0805.1726 [gr-qc]]
  * [44] L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, Phys. Rev. D 75 (2007), 083504 doi:10.1103/PhysRevD.75.083504 [arXiv:gr-qc/0612180 [gr-qc]]
  * [45] A. Chatterjee, R. Roy, S. Dey and A. Bandyopadhyay, [arXiv:2310.05578 [gr-qc]]
  * [46] D. Samart, B. Silasan and P. Channuie, Phys. Rev. D 104 (2021) no.6, 063517 doi:10.1103/PhysRevD.104.063517 [arXiv:2104.12687 [gr-qc]]
  * [47] L. Amendola, D. Polarski and S. Tsujikawa, Phys. Rev. Lett. 98 (2007), 131302 doi:10.1103/PhysRevLett.98.131302 [arXiv:astro-ph/0603703 [astro-ph]].
  * [48] L. Amendola, D. Polarski and S. Tsujikawa, Int. J. Mod. Phys. D 16 (2007), 1555-1561 doi:10.1142/S0218271807010936 [arXiv:astro-ph/0605384 [astro-ph]]
  * [49] L. Amendola and S. Tsujikawa, Phys. Lett. B 660 (2008), 125-132 doi:10.1016/j.physletb.2007.12.041 [arXiv:0705.0396 [astro-ph]]
  * [50] V. Faraoni, [arXiv:0810.2602 [gr-qc]]
  * [51] S. Tsujikawa, doi:10.1007/978-90-481-8685-3_8 [arXiv:1004.1493 [astro-ph.CO]]
  * [52] B. Li and J. D. Barrow, Phys. Rev. D 75 (2007), 084010 doi:10.1103/PhysRevD.75.084010 [arXiv:gr-qc/0701111 [gr-qc]]