Dynamics of viable f ( R ) dark energy models in the presence of Curvature-Matter interactions

,


Introduction
The measurement of redshift and luminosity distance for type Ia supernovae events [1,2], played a pivotal role in confirming the transition of the universe from a decelerated phase to its current accelerated phase during its late-time evolution.These findings are substantiated by the analysis of temperature anisotropies in the cosmic microwave background data from the WMAP mission [3,4] and the detection of baryon acoustic oscillations [5].The general theory of relativity posits that the accelerated expansion is driven by a hypothetical energy density component characterized by negative pressure, commonly referred to as dark energy.Investigating the origin and nature of dark energy responsible for the ongoing cosmic acceleration remains a prominent unsolved mystery in contemporary cosmology.The present-day cosmic acceleration cannot be accounted for by the standard equation of state, expressed as ω = p/ρ, where p and ρ represent the pressure and energy density of the conventional universe contents, such as radiation and matter.In fact, we need a 'yet-to-be-identified' component characterized by negative pressure, resulting in an equation of state ω < −1/3, to provide an explanation for the accelerated expansion.
Initially, Einstein introduced a cosmological constant term Λg µν into his field equations of general relativity to support his static universe theory.However, he later abandoned it in favor of Hubble's observations of an expanding universe.Nevertheless, in the latter part of the last century, as late-time cosmic acceleration was discovered, the cosmological constant term regained attention due to its potential to offer a straightforward solution to Einstein's field equations, resulting in an accelerated expansion.This gave rise to the Λ-CDM model, where 'CDM' refers to cold dark matter.Unfortunately, this model grapples with the coincidence problem [6] and the fine-tuning problem [7].These issues drive the exploration of alternative models to account for dark energy.
Dynamical systems techniques are powerful tools for exploring cosmic evolution in generic cosmological models, as well as for investigating specific cosmological solutions.Dynamical systems analysis has recently been employed to assess the stability of various scalar field dark energy models in comprehensive reviews [37][38][39][40][41], as well as in the context of f (R) gravity [42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57].A systematic and comprehensive exploration of dynamical systems within f (R) gravity theories has been carried out in a series of research endeavors documented in [58][59][60][61] and references therein.These studies aimed to identify and classify models that faithfully portray the correct cosmological evolution among a diverse set of f (R) gravity models.Consequently, they contributed to the construction and development of 'cosmologically viable' f (R) gravity models, characterized by trajectories in the dynamical system phase space that align with observed cosmic phenomena.These viable f (R) gravity models also conform to local gravity constraints, as extensively discussed with examples in references [42][43][44].What makes these models appealing is their ability to provide cosmic acceleration without requiring a cosmological constant or the introduction of dark energy as an additional field component in the universe.Instead, dark energy is 'curvature-driven' -the dynamics described by the modified field equations incorporate the function f (R) of Ricci scalar (R), which possesses the potential to drive cosmic acceleration, even when facing stringent constraints from both local gravity tests and various observational limitations.
In this study, we explore a scenario in which curvature-driven dark energy interacts with (dark) matter, all within the framework of viable f (R) gravity models.We investigate whether the inclusion of these interactions in viable f (R) models has a significant impact on achieving accurate cosmological acceleration when compared to f (R) models that do not include matter-curvature interactions.In f (R)-gravity theories which maintain a minimal coupling to matter [30][31][32][33][34], the dynamics originate from an action composed of two components: a modified 'geometric' part obtained by replacing R with f (R) in the Einstein-Hilbert action, and a 'matter' component arising from the conventional constituents (matter and radiation) of the universe.The involvement of the metric g µν in the matter component ensures minimal coupling with the geometry.The modified field equations can be restructured to assume the appearance of the standard field equations, expressed as G µν = 8πGT effectively serves as a representation of the stress-energy tensor, simulating a fluid-equivalent of the curvature function, while T µν is the adjusted stress-energy tensor for 'matter and radiation' in the presence of gravity modifications arising from f (R).In FLRW spacetime background, the conservation of T (tot) µν leads to a continuity equation that expresses the conservation in terms of the energy density and pressure for the comprehensive fluid, comprising radiation, matter, and curvature components.The continuity equation for each individual sector may allow for the inclusion of source terms, as long as the overall continuity equation for the comprehensive fluid remains valid without any source term.In this study, we considered radiation to be decoupled from matter and curvature and proposed that interactions between the matter and curvature components introduce a source term Q in their respective continuity equations, with opposite signs.The term Q serves as a measure of rate of energy exchange between the matter and curvature sectors owing to their interaction.In prior research [58], the expression for Q has been determined on the multiplicative function of matter density and Hubble parameter.In this subsequent analysis, we take a further step by formulating the interaction term as Q = κ 2 α 3H ρm ρ curv where, α signifies the strength of the coupling between matter and curvature with ρm and ρ curv as their respective energy densities.
The incorporation of interactions between matter and curvature thus introduces the additional parameter α in into the analytical framework.This coupling parameter, in conjunction with the parameters defining the function f (R), becomes intricately entwined within the evolution equations, playing a substantial role in shaping the evolutionary dynamics in the matter-curvature coupled scenario.We employed dynamical analysis technique as a method to explore this interacting curvature-matter scenario in the context of viable f (R)-gravity models.This process entailed setting up the representative autonomous equations in terms of appropriately defined dynamical variables describing the cosmic evolution in this context.The essential components of our study involved identifying the fixed points of the system and examining their stability using the technique of linear stability analysis.These are crucial for achieving a comprehensive understanding of critical aspects in cosmic evolution, particularly those influenced by the interactions between curvature and matter, as considered in the present framework of our study.In this investigation, we've opted for two specific f (R) models, namely, the generalised Λ-CDM model: Both models have the potential to drive cosmic acceleration while maintaining their cosmological viability.In general, the fixed points and characterisation of their stability are de-pendent on the parameters involved in the analysis, such as α, (b, c) or n.Additionally, values of the density parameters for various components and the equation of state for the overall fluid at these fixed points provide crucial insights into the distinctive critical phases associated with each point.We obtain the parameter limits within which a given fixed point may exhibit stable, unstable, or saddle-like behaviour.We also examine whether it can serve as an acceleration solution while staying within the boundaries of various cosmological constraints, as elaborated in subsequent sections.
Finally, we explore how curvature-driven dark energy influences various phases of cosmic evolution within the framework of matter-curvature interactions.To achieve this, we examined the evolution of matter-to-curvature energy density ratio (denoted as r mc in the text), as well as two additional cosmographic parameters: deceleration (q) and jerk (j).The evolution of r mc provides quantitative information about the dominance of curvature over matter at different stages of cosmic evolution.On the contrary, the deceleration parameter, along with the jerk, provides a detailed account of the kinematic aspects of evolution.These investigations shed light on the potential scenarios that may arise within the context of matter-curvature interactions.The observed evolution of the deceleration and jerk parameters suggests that all the mentioned models ultimately converge toward the Λ-CDM model in the distant future.
The paper is organised as follows: In Sec.[2], we establish a theoretical framework that addresses the interaction between the curvature and matter sectors within the context of a flat FLRW metric.Sec.[3] is dedicated to formulating the autonomous equations of the dynamical system, which incorporates the interaction between matter and curvature-driven dark energy within the framework of f (R) gravity.Within this section, we also discussed the two specific f (R) models considered for this study.We present the fixed points of the autonomous system for each model and discuss their stability characteristics.Within Sec.[4], we examine the influence of curvaturedriven dark energy on different phases of cosmic evolution.We summarize our conclusions in Sec.[5].
2 Theoretical framework of interacting curvature-matter scenario in f (R) gravity Within the framework of f (R) theory of gravity that incorporates minimal curvature-matter coupling, the action can be expressed as follows [30,31] where κ 2 = 8πG, f (R) is an arbitrary function of the Ricci scalar R and L M commonly referred to as the 'matter' Lagrangian, represents the Lagrangian density associated with both the radiation and matter (baryonic and CDM) components of the universe.The symbol g represents the determinant of the space-time metric tensor g µν .In metric formalism, the variation of this action with respect to the field g µν gives the corresponding modified field equation as Here, F (R) ≡ df /dR, and T (M) µν is the stress-energy tensor corresponding to radiation and matter components given by Eq. ( 2) may be reconfigured as [30] G where, G µν is the Einstein tensor and and In its reformulated form, the modified field equation (4) describes a FLRW universe featuring a fluid with a comprehensive energy-momentum tensor denoted as . The influence of the function f (R) is readily apparent in both of these distinct components T is exclusively determined by f (R) and its higher derivatives and it vanishes for f (R) = R.This component effectively acts as a portrayal of the stress-energy tensor, mimicking a fluid-equivalent of the curvature function f (R).

The component T (M)
µν is obtained by modifying T (M) µν using a functional multiplier of 1/F (which equals 1 when f (R) = R).It serves as the revised stress-energy tensor for 'matter and radiation' in presence of gravity modifications resulting from f (R).In the context of a flat FLRW spacetime background, when we regard the matter and radiation content as a perfect fluid, the (unmodified) stress-energy tensor T (M) µν becomes diagonal with entries being completely determined by the combined energy densities (ρ m + ρr ) of dark matter (ρ m ), radiation (ρ r ), and the radiation pressure ( Pr = ρr /3).The matter component is considered to be non-relativistic pressureless dust.Consequently, the modified stress energy tensor µν corresponds to a fluid with energy density (ρ m + ρ r ) and pressure P r , where (ρ r , ρ m , P r ) ≡ (ρ r /F, ρm /F, Pr /F ).Note that, ρ i , P r (i = r, m) are always positive because of the inherent positivity of ( ρi , Pr ) and positivity of F as well -a necessary condition for any cosmologically viable f (R) model (discussed later in Sec.[3.2]).However, under such considerations, the '00' and 'ii' components of the modified field equation (4) respectively yield the following modified versions of the Friedmann equations; where a is the FLRW scale factor, H ≡ ȧ/a (the Hubble parameter) and dot (•) denotes derivative with respect to cosmic time.The quantities ρ curv and P curv are defined through the following equations where ′ denotes derivative with respect to R. It follows from eq. ( 5) that, in a flat FLRW spacetime, the stress-energy tensor T curv µν representing curvature-fluid can be equated with that of an ideal fluid with energy density and pressure given by ρ curv and P curv as provided in eqs.(9) and (10).Also note that, we may express eq.(7) as where Ω m , Ω r and Ω curv are modified form of density parameters defined by following equations: It's also worth mentioning that the specific case of f (R) = R (F = 1, F ′ , F ′′ = 0) leads to ρ curv = 0, P curv = 0, (ρ i , P r ) = (ρ i , Pr ) implying reduction of modified Friedmann eqs.( 7) and (8) to conventional Friedmann equations.
Taking divergence of the both sides of eq. ( 4) and use of Bianchi's identity lead to the conservation equation of total stress-energy tensor T tot µν = T (M ) µν + T curv µν : In FLRW spacetime, eq. ( 13) takes the form of a continuity equation of the all-inclusive fluid with energy density ρ tot ≡ ρ r + ρ m + ρ curv and pressure P tot ≡ P r + P curv as ρtot + 3H(ρ tot + P tot ) = 0 or ρtot + 3Hρ tot (1 where  µν as its further sub-components) to be separately conserved.This feature provides scope for incorporating 'interactions' between the matter and curvature sectors.
Considering the negligible contribution of the radiation component during the late time stages of cosmic evolution, we assume absence of interactions between radiation and combined mattercurvature sectors while simultaneously preserving the possibilities of interactions between matter and curvature sectors.With such considerations, we posit the conservation eq.(13) as a system of equations: (I) ∇ µ T µν .In FLRW spacetime, (I) implies ρr +4Hρ r = 0 (ρ r ∼ a −4 ) and equation set (II) corresponds to non-conservation equations of the form: with a source term Q, which typically signifies a time-varying function characterizing the instantaneous rate of energy exchange between the curvature and matter sectors, reflecting their interactions.Interaction between dark matter and curvature-driven dark energy can either be introduced at the Lagrangian level, guided by field theoretic considerations, or at a phenomenological level by modelling the non-zero source term Q (in eqs.( 15),( 16)) in temrs of suitable cosmological quantities.We adopt the latter approach, chosing a specific form of Q given by Q = κ 2 α 3H ρm ρ curv , where the coupling parameter α is dimensionless maintaining the dimensional consistency in both eqs.(15) and (16).In this context, we may mention that, in a previous study [58], the impact of interaction was explored using a functional form of the source term that exclusively depended on the energy density of the matter sector, Hubble parameter, and coupling constant in a multiplicative form.Expanding upon this framework, our chosen form of the source term involves the multiplication of the energy densities of both sectors capturing the instantaneous effects of both sectors on the source term with a dimensionless coupling constant.Our analysis reveals that such an interacting term also has the potential to address the coincidence problem, as discussed in the concluding section of the paper.
We imposed the constraint ρ curv > 0 in attributing the sense of 'energy density' (of the fluidequivalent of the curvature) to ρ curv .This constraint restricts f (R) models to adhere to the condition RF − f > 6H ṘF ′ as evident from eq. ( 9).It also constraints the coupling parameter α as ρ curv is interconnected with it through the source term Q as described by eq. ( 15).Under such considerations, each of the modified density parameters Ω r , Ω m and Ω curv remain positive and collectively subject to the constraint given by Eq. ( 11).This allows us to navigate and explore the parameter space relevant to the scenario of curvature-matter interaction by varying the parameter Ω m within the range 0 ⩽ Ω m ⩽ 1, eventually leading to obtaining constraints on the parameters of the model.
The acceleration phase (ä > 0) corresponds to ω tot = P tot /ρ tot < −1/3 as indicated by eq. ( 8).F (R) and its derivatives are involved in ω tot through ρ curv and P curv which, in turn are linked to the coupling parameter α via Equations ( 15) and ( 16).The values of the parameter α and those involved in the parametrisation of F (R) that satisfy the condition ω tot < −1/3 during the latetime phase of cosmic evolution corresponds a 'Dark Energy' scenario driven by matter-curvature interactions.Non-phantom dark energy scenarios adhere to the constraint −1 < ω tot < −1/3.
3 Dynamical analysis of f (R) models with curvature-matter interactions

Autonomous equations of dynamical system
We employ dynamical analysis as a tool to examine the curvature-matter interaction scenario within the framework of cosmologically viable f (R) gravity models.In this subsection we furnish the set of autonomous equations in terms of suitably constructed basic dynamical variables which serve to characterize the cosmological evolution within the context under consideration.This sets the ground for the application of techniques of dynamical analysis in our current context.
The basic dynamical variables are [44] By introducing the dimensionless parameter denoted as N = ln a to address temporal variations and applying the given definitions of dynamical variables (X i 's), we can express eqs.( 9), (15), and ( 16) which capture the evolutionary dynamics within the curvature-matter interaction context, as the following set of 4 autonomous equations, forming a 4-D dynamical system: where and in this context another parameter r is introduced as These dimensionless parameters m and r both depend on R through f (R), allowing us to write m as a function of r i.e. m = m(r).Each specific functional relation m = m(r) corresponds to a distinct class of f (R) models.When the coupling parameter α is set to zero, the above set of equations reduces to the autonomous equations corresponding to the scenario that does not take into account interactions between the curvature and matter sectors.Also note that, if f (R) = R, then the parameter m (as well as X 1 ) goes to zero, rendering the construction of the dynamical system infeasible.
By virtue of eqs.( 11) and ( 12) the dynamical variables in eq. ( 17) are constrained by the equation where Ω m is subject to the inequality 0 ⩽ Ω m ⩽ 1 and we can additionally write Also, the 'grand' EoS parameter ω tot can be expressed as EoS parameter associated with the solo curvature sector [57] can be written as, The ratio of matter to curvature energy density, denoted as r mc ≡ Ω m /Ω curv , along with two additional cosmographic parameters -the deceleration (q ≡ −aä/ ȧ2 ) and the jerk (j ≡ ... a /aH 3 ) -serves as valuable indicators for evaluating various facets of cosmic dynamics.During late cosmic epochs, the ratio r mc , approximately reflects the ratio of 'dark matter' to 'dark energy' densities, as observations suggest baryonic contributions to be significantly smaller compared to dark matter contributions in the present-day universe.The parameter r mc thus holds the potential to address the issue of the apparent coincidence between dark matter and dark energy densities in present-day universe (the cosmic coincidence problem).The deceleration (q) and jerk (j) parameters directly describe the 'kinematic' characteristics of cosmic evolution.A negative q signifies an accelerating universe, which is a hallmark of dark energy dominance.The constant q = −1 value signifies presence of cosmological constant.A switchover from a positive to a negative value of q marks the transition point from a decelerating to an accelerating phase cosmic expansion.Furthermore, varying q values correspond to dynamical dark energy models, indicating the evolving nature of the universe's acceleration.The jerk parameter, j is a measure of rate of change of acceleration (i.e. of q) and thus captures even finer temporal details of cosmic acceleration.Both q and j encapsulate kinematic aspects of cosmic evolution which leave distinct imprints on the large-scale structure of the universe.Consequently, these parameters have a crucial role to play in our endeavours to understand the impact of dark energy across the different stages of cosmic evolution.Furthermore, these parameters are instrumental in facilitating comparisons and distinctions among different 'dark energy models' that utilize various mechanisms to trigger cosmic acceleration.
In the context of the matter-curvature interaction scenario considered in this study, the expressions for the three parameters (r mc , q, j) can be derived in terms of the dynamical variables (X i 's) as In Sec.3.2 we have computed the temporal evolution of these three parameters for a FLRW universe.
The fixed points in the 4D dynamical system described above correspond to solutions for X i 's of dX i /dN = 0, (i = 1, • • • , 4)).These stationary solutions and the assessment of their stability are crucial for a comprehensive grasp of critical aspects in cosmic evolution driven by curvaturematter interactions, as explored within the scope of this study.To investigate the system's behaviour around these critical points, we employ the linear stability analysis which involves performing a first-order Taylor expansion of the autonomous equations derived above, which have a generic expression of the form: s given by right hand side of the set of autonomous equations derived in eqs.(18,19,20,21).The Jacobian of the linear transformation ⃗ X → ⃗ f is the 4 × 4 matrix By examining the eigenvalues of J computed at the critical points we may deduce the nature of the stability of the fixed points and classify them accordingly.When all the eigenvalues have negative (positive) real parts, the fixed point is classified as asymptotically stable (unstable).On the other hand, if any pair of eigenvalues occur with a relative opposite sign in their real parts, the corresponding fixed point is a saddle point.However, if any of the eigenvalues approaches zero at the fixed point, the linear stability theory proves inadequate, necessitating the application of center manifold theory for a more insightful exploration of the characteristics of these critical points.However, in our analysis detailed in later sections, we did not find any non-hyperbolic fixed points within the scope of the various f (R) models examined in this study.This suffices for adhering to linear stability analysis for the models considered here without necessitating an extended exploration using methods beyond linear stability analysis [47].

On the choice of f (R)-models
The main objective of this study is to investigate and analyze the dynamics of f (R)-driven dark energy models taking into account interactions between the curvature and matter sectors.
To carry out this investigation, we have selected specific f (R) models that can induce cosmic acceleration while maintaining their cosmological viability.In the context of the metric formalism employed in this study, any f (R) function that is cosmologically viable must satisfy a set of stringent conditions which have been comprehensively discussed in [62,63].Below, we enumerate these conditions and provide the underlying rationale for each of them.
(iii) f (R) is close to (R − 2Λ) for R ≫ R 0 : Consistency with local gravity tests and presence of matter dominated era.
Here R 0 represents the value of the Ricci scalar at present epoch.The relationship m = m(r) between the dimensionless parameters m and r, as outlined in eqs.( 22) and ( 23), is determined by the choice of f (R).Conversely, a specific form of m(r) corresponds to a specific f (R).
We investigate the autonomous system furnished in Sec.3.1, focusing on two separate viable scenarios: The first scenario relates to a f (R) resulting in a constant m, while the second concerns a form of f (R) where m is a function of r, represented as m(r).We designate these scenarios as (A) and (B) and discuss their specifications below.
(A) We understand that a scenario with a constant m can be achieved through the modified gravity model given by . This model has been previously examined without taking into account matter-curvature interactions in [44,61].This was proposed with the aim of extending the Λ-CDM model to address local gravity constraints justifying it as a viable f (R) model.We will refer to this scenario as 'generalised Λ-CDM model' or sometimes simply as 'model-A' from here on.When constrained to c ⩾ 1 and bc ≈ 1, this model converges to the Λ-CDM model with m = 0, ensuring a viable cosmological evolution.The m and r parameters for the model are given by m and consequently, the parameter m can be expressed as m(r) = 1−c c r+b−1.A dynamical analysis (see [44,61] for details) of this model indicates that at stable points r = −1 − m, which results in m = −1 + bc, a constant value.
(B) This scenario involves consideration of m as a specific function of r given by m(r) = n(1+r) r which can be realised for the power law form f (R) = R − γR n with (γ > 0, 0 < n < 1).We refer to this scenario in subsequent texts as 'power law model' and sometimes for convenience as 'model-B'.As comprehensively discussed in [44,64], the power law model correctly describes cosmic evolution within framework of non-interacting curvaturematter scenarios and satisfies the condition for cosmological viability as well for the above mentioned range of γ and n.This form of f (R) gives m = γ(n−1)nR n R−γnR n and r = R−γnR n R−γR n , the elimination of γ from which leads to the form m(r) = n(1+r) r .
In the context of our study, choice of the above two cases: (A) m = −1 + bc = constant and (B) m = m(r) = n(1+r) r are intended to explore implications of the two above mentioned f (R)-models in presence of matter-curvature interactions.

Stability analysis of fixed points in generalised Λ-CDM and power law models
In the context of model-A and model-B, we identify the fixed points of the autonomous system described in Sec.[3.1] representing curvature-matter interactions within the framework of f (R) gravity.The fixed points are obtained by solving the equations dX k /dN = 0 (k = 1, 2, 3, 4), where dX k /dN 's are given by eqs.(18,19,20,21) for model-B.The fact that m is constant for model-A and r-dependent for model-B with r = X 3 X 2 i.e. expressible completely in terms dynamical variables X 2 and X 3 , enables us to formulate the system as a closed set of four autonomous equations involving four dynamical variables and the constant parameters (α, b, c, n) without the necessity for additional dynamical variables.
The coupling parameter α serves as a common parameter in the analytical framework across both models due to the consideration of matter-curvature coupling in our analysis.Moreover, the f (R) model includes distinct model parameters: (b, c) for model-A and (n) for model-B.We can also evaluate the values of density parameters (Ω r , Ω m , Ω curv ), grand EoS parameter ω tot and EoS parameter associated with curvature sector ω curv at the fixed points using eqs.( 17), ( 24), ( 25), ( 26), (27).As per our considerations mentioned in Sec.[2] each of the modified density parameters are positive and subject to the constraint given by eq. ( 11).Adherence to these constraints determines which fixed points hold cosmological significance.Therefore, providing a description of fixed points and their stability characteristics, which are dependent on model parameters, requires appropriate referencing of the relevant constraints within the model parameter space.Moreover, the density parameter values and the EoS parameter value at each fixed point provide crucial insights into characterizing the distinct critical phase associated with that point.If ω tot < −1/3 for any fixed point, it represents cosmic acceleration.Specifically, if the value falls within the range of −1 < ω tot < −1/3, the fixed point corresponds to a 'nonphantom' acceleration and when ω tot = −1, it corresponds to the de-Sitter acceleration solution.
The density parameter values at any given fixed point provide indications of the predominant component during the cosmic phase represented by that specific fixed point.
A total of 10 (real) fixed points are identified in model-A, while model-B comprises 7 fixed points.To distinguish them within the text, we label them as follows: For scenario (A) : {P 1 , P 2 , P 3 , P 4 , P 5A , P 6A , P 7A , P 8A , P 9A , P 10A } For scenario (B) : {P 1 , P 2 , P 3 , P 4 , P 5B , P 6B , P 7B } These fixed points, along with the corresponding values of Ω m and ω tot for both models, are displayed in Tabs.[1], [2], and [3].The arrangement of distributing the fixed points across three tables employs a classification scheme to facilitate a comparison of the characteristics attributed to the fixed points in both models.Fixed (X 1 , X 2 , X 3 , X 4 ) Ω m ω tot ω curv points P 1 (−4 , 5 , 0 , 0) 0 Common fixed points of model-A and model-B.The values of Ω m , ω tot and ω curv at these fixed points are also presented.
The 4 points (P 1 , P 2 , P 3 , P 4 ) presented in Tab.[1] serve as common fixed points in both model-A and model-B and they are also independent of the parameters of their respective f (R)-model parameters.Notably, the first two points P 1 and P 2 are even independent of the coupling parameter α implying that they would also persist as fixed points in the dynamical systems associated with non-interacting curvature-matter scenarios, regardless of any applied modifications to the gravity function f (R).Fixed points P 2 and P 4 correspond to ω tot = −1, indicating an accelerating de-Sitter solution.The α-dependence of the other two points (P 3 and P 4 ) arises directly from the consideration of matter-curvature interaction in the analysis, and they would not serve as (finite) fixed points in scenarios that do not account for such interactions.This can Fixed (X 1 , X 2 , X 3 , X 4 ) Ω m ω tot ω curv points (1 , 0 , 0 , 0 ) 0 P 7A (0 , 0 , 0 , 1 ) 0 undefine Table 2: Fixed Points of model-A that are independent of f (R) model parameters (b, c) but do not serve as fixed points of model-B.The values of Ω m , ω tot and ω curv at these fixed points are also presented.
be understood by observing that the values of X 2 associated with these two points diverge as the matter-curvature coupling (α) approaches zero.It is worth noting that for a critical value of α = −3, the points P 3 and P 4 converge to the points P 1 and P 2 respectively.
The 3 fixed points (P 5A , P 6A , P 7A ) in model-A, listed in Tab.[2], appear to be independent of the corresponding f (R) model parameters (b, c) but do not function as fixed points in model-B, which pertains to a different f (R) model.The emergence of these points can be attributed to the constancy of the value of m = −1 + bc within the dynamical equations for model-A, in contrast to the situation in model-B where m = n(X 2 +X 3 ) itself enters as a function of the dynamical variables, as discussed earlier.It is worth noting that, the fixed point P 7A corresponds to a state characterized by complete radiation dominance (Ω r = 1, ω tot = 1/3) and none of the three points listed in Tab.[2] are representatives of cosmic acceleration.
The remaining fixed points (P 8A , P 9A , P 10A ) in model-A, as well as (P 5B , P 6B , P 7B ) in model-B, depend on the parameter(s) of their respective f (R) models.They are collectively presented in Tab.[3].Since the expressions for ω tot at these fixed points depend on their respective model parameters, they would correspond to non-phantom accelerating solutions only if appropriate parameter values exist for which the condition −1 < ω tot < −1/3 is satisfied.
Regarding the fixed points identified in our analysis, it's worth noting that if we were to exclude matter-curvature interactions from our considerations by setting α = 0, the resulting dynamical system would consist of 8 fixed points model-A and 5 fixed points for model-B, as detailed below.
Table 3: Fixed points in model-A and model-B that are dependent on parameters of corresponding f (R) models.The values of Ω m , ω tot and ω curv at these fixed points are also presented.Displays a saddle-like behavior, lacking both stability and cosmic acceleration for any parameter values.In the limit as bc → 1, it is purely dominated by radiation.Otherwise, all three density components Ω r , Ω m , and Ω curv are present.

Table 4: Stability and cosmological features of the fixed points of model-A
Note that, here we have examined the stability properties of fixed points of the relevant autonomous system that occur with finite coordinate values.However, when dealing with scenarios involving non-compact phase spaces, it becomes crucial to explore the possibility of fixed points occurring at infinity and evaluate their stability criteria.The methodology for investigating features in the asymptotic regime is extensively discussed in [52,53].One may observe that all the dynamical variables of this system are bounded owing to the imposition of of the various energy density constraints (0 ≤ Ω m , Ω r ≤ 1), −1 < ω tot < 1/3 throughout the evolutionary era from radiation to dark energy domination and viable f (R) dark energy condition.This eliminates the necessity to explore fixed points at infinity.We examine the stability characteristics of the various fixed points using linear stability analysis, as outlined in Sec.[3.1].This entails an examination of the eigenvalues of the Jacobian matrix J (Eq. ( 31)).The expression for these eigenvalues typically involves the coupling parameter α and the parameters of the relevant f (R) model.Therefore, to explore the stability properties of the fixed points, it becomes essential to examine the eigenvalues of J across the parameter space that is pertinent to the particular model being studied.Due to the substantial and intricate analytical expressions of most eigenvalues for arbitrary values of model parameters, we opt not to include their detailed expressions in the article.We perform a thorough scanning of the parameter space allowing for a wide range of values for α, (b, c), while concurrently enforcing the constraint 0 < n < 1, which is essential to uphold the cosmological viability of the power law model.We obtain the (i) parameter space limitations under which a given fixed point demonstrates stable / unstable or saddle-like behaviour.We also identify (ii) the parameter range within which the fixed point represents an acceleration solution of a specific kind -phantom/ non-phantom or de-Sitter.Si-multaneously, we isolate (iii) the parameter constraints that ensure that Ω m remains within the range of 0 < Ω m < 1.
We perform a comprehensive scrutiny of the parameter constraints derived from the three schemes (i), (ii), and (iii) and offer a succinct overview of the stability characteristics and other cosmological aspects related to each of the fixed points.This summary, including description of relevant parameter ranges, is presented in Tab.[4] and [5] for model-A and model-B respectively.Moreover, for visualization purposes, we specifically selected three fixed points, namely, P 8A , P 9A , and P 6B , to illustrate the constraints within the parameter space.In fig.
[1], we present three distinct panels, each illustrating the permissible parameter ranges for individual points, wherein they separately (i) exhibit stability, (ii) demonstrate cosmic acceleration (−1 < ω tot < −1/3), and (iii) satisfy the condition of 0 < Ω m < 1.Note that, the condition 0 < Ω m < 1 is not applicable for the case of point P 8A as it corresponds to an Ω m value that is independent of the model parameters (Ω m = 0 for P 8A , see Tab. [3]).The areas of overlap between (i) and (ii) for point P 8A and among (i), (ii), and (iii) for point P 9A define the specific parameter domains for these two points within which 'stable cosmic acceleration is achieved through curvature-driven (non-phantom) dark energy'.No intersection of the three domains corresponding to (i), (ii) and (iii) are found for point P 6B .

Fixed Stability & characteristics points P 1
For any value of the coupling parameter α, this point exhibits saddle-type behavior within the pre-defined range of n.This point neither represents matter dominated epoch nor dark energy epoch.
This ponit is spiral stable if α > −3 and 0 < n < 32  25 .Within the viable rage of n this point signifies a stable accelerating de-Sitter solution, predominantly governed by curvature contribution.P 3 This point exhibits simililar behavior like P 1 .Saddle point behaviour with no interesting cosmological sequences regardless of the values of α.P 4 This point is also stable for α ⩽ −3 within the viable range of n. de-sitter like acceleration can be realised in presence of both matter and curvature components for any value of α except α = −3.For α = −3, this de-sitter attractor point has driven by the curvature only like P 2 .

P 5B
Stable for α > −3, but there is no common region in the parameter space of n and α where stability and accelerating condition simultaneously satisfy within the viable range of n.P 6B Stable for α < 0.8, but like P 5B this points represents stable non-accelerarting fixed points illustrated in fig.
[1](c).But in n → 1 limits this critical point embodies the conventional saddle matter dominated scenario.P 7B The system exhibits stable or saddle-like behavior, contingent on the value of α within the feasible range of n.In this viable range, ω tot is greater than 1 3 .However, in the limit as n approaches 1, it is solely dominated by radiation.Otherwise, there is no cosmologically interesting scenario has been observed here.

Table 5: Stability and cosmological features of the fixed points of model-B
It's important to note that, within the framework of f (R)-gravity scenarios featuring curvaturematter interactions, the fixed points P 2 and P 4 (which are common to both model-A and model-B), as well as fixed point P 8A and P 9A in model-A, are the only fixed points that represent 'stable acceleration solutions.'These solutions remain stable within the specified range of α values, as outlined in Tabs.[4] and [5].Specifically, P 2 and P 4 correspond to stable de-Sitter attractors, while P 8A and P 9A correspond to stable acceleration driven by non-phantom dark energy.It's important to emphasize that there are no other points that exhibit both stability and acceleration while adhering to the constraint 0 < Ω m < 1 for any combination of parameter values.
4 Evolutionary dynamics of curvature-matter coupled scenario in the context of cosmological and cosmographic variables In the context of f (R) gravity models featuring matter-curvature interactions, we analyzed cosmic evolution across the entire timeline -from the radiation epoch to the matter epoch and finally to the late-time cosmological era.To comprehend the complete dynamics of the intertwined curvature-matter scenario, we introduced three cosmological parameters, namely the density parameters associated with curvature and matter (Ω i ), the grand equation of state (EoS) parameter (ω tot.), the ratio of matter density to curvature density (r mc ), and two crucial cosmographic parameters -deceleration (q) and jerk (j).These cosmographic parameters, previously introduced in eqs.( 28), (29), and (30), play a pivotal role in providing a quantitative characterization of the evolutionary progression.
We numerically solve the system of autonomous eqs.(18,19,20,21) for the dynamical variables X i in the two cases -one with a constant m corresponding to model-A and the other with m(r) = n(1+r) r , where r = X 3 /X 2 corresponding to model-B.Due to the non-linearity of the coupled autonomous equations, the solutions are highly sensitive to the initial conditions necessary for solving this set of equations.However, to visually illustrate the potential evolutionary dynamics in f (R)-gravity models with matter-curvature interactions, we establish specific valid boundary conditions that allow us to obtain solutions for X i based on carefully selected parameter values.These obtained solutions are then employed in equations ( 24), (25), and (26) to compute the temporal profiles of density parameters (Ω m , Ω curv ) and the equation of state parameter (ω tot ).Additionally, we determine the temporal evolution of the triad of parameters (r mc , q, j) using eqs.( 28), (29), and (30).The depicted outcomes in fig.[2] model-A and fig.
[3] for model-B demonstrate the evolution of Ω m , Ω curv , ω tot , and r mc in the left panel, while showcasing the values of q and j in the right panel.This depiction covers a range of values for N = ln a, extending from 0.001 to 1000.The chosen parameter values (α, bc or n) for these representations are carefully selected to ensure that, they correspond to stable accelerating solutions and conform to the constraint 0 < Ω m < 1.The specific parameter values are provided in the figure captions for their respective cases.(scenario (B)).Left Panel: Plots of r mc , Ω m , Ω curv and ω tot vs N (= ln a) computed with n = 0.9, α = −5, Right Panel: Plots of cosmographic parameters q and j computed with n = 0.9, α = −2.
The graphs of Ω curv , Ω m , and ω tot suggest the possibility of early cosmic phases where curvature energy density dominates over (dark) matter, yielding a positive ω tot (with the value of 1  3 for the cases presented in figs.[2] and [3]) and exhibiting a radiation-like EoS.As time progresses, curvature energy density diminishes while matter-energy density grows, and their values tend to converge.At the matter-dominated epoch (where, ω tot.= 0), the values of the two densities become equal, triggering a shift in the nature of their temporal behaviour, resulting in the rise of curvature energy density and the decline of matter energy density.It's worth noting that curvature energy density consistently exceeds matter energy density throughout the depicted evolutionary phases in the left panel of figs.[2] and [3].During moments when their values come together, the equation of state parameter ω tot , which had remained constant and positive since the initial epochs experiences a rapid decline and swiftly falls below − 1 3 .This indicates the onset of late-time cosmic acceleration.For model (A), ω tot parameter can't cross the phantom barrier line, instead it saturates near the value of −0.8.As this phase unfolds, the values of Ω curv and Ω m once again separate, with curvature density regaining dominance over matter energy density.This dynamic facet of cosmic evolution is common to both scenarios (A) and (B), as can be inferred from figs.[2] and [3].
With a more in-depth comparison of figs.[2] and [3], one may also recognize the distinctive im- pacts that the corresponding f (R) modifications in model-A and model-B exert on the generic features of cosmic evolution mentioned above.In model-A, as cosmic acceleration initiates, the rapid decline of ω tot does not cross the value of −1, but rather gracefully saturates at a value slightly above −1, indicating a non-phantom nature of curvature-driven dark energy.In contrast, in model-B, the sharp fall of ω tot takes it below −1, but it subsequently rebounds and eventually settles at the value of −1, never exceeding it.This feature implies that within the framework of the power-law model with matter-curvature interactions, the curvature-induced dark energy responsible for cosmic acceleration can exhibit characteristics of phantom darkenergy.Nevertheless, it's important to note that such a scenario does not exhibit stability since ω tot eventually saturates -1 with the passage of time.Another notable distinction lies in the possibility for late-time cosmic acceleration to be purely curvature-dominated (Ω curv = 1) in model-A.On the contrary, in model-B, cosmic acceleration can occur with both matter and curvature energy densities present in approximately equal proportions (on the order of ∼ 1), with an edge towards curvature dominance over matter.
Furthermore in model-B, where f (R) follows a power law form indicated by m(r) = n(1 + r)/r, as the universe progresses into its late-time acceleration phase, the temporal profiles illustrating Ω curv and Ω m display intricate and non-trivial patterns, rather than a simple, monotonic progression, until they eventually stabilize in the distant future.In contrast, these traits are absent in the constant m scenario (model-A).This distinction can be attributed to the differences resulting from the specific forms of f (R) in each scenario.
The ratio of matter to curvature density r mc ≡ Ω m /Ω curv , serves as a gauge of the extent to which one component dominates over the other.The value r mc = 0 signifies complete curvature dominance, while r mc ≈ 1 suggests a coincidence of matter and curvature energy densities in the universe.Therefore, the evolution of the parameter r mc , as depicted in fig.
[3] (model-B), provides quantitative insights into the dominance of curvature over matter at different stages of cosmic evolution in the respective scenarios.The impact of the curvaturematter coupling strength on cosmic evolution becomes apparent when we compare the profiles of r mc computed using different values of the interaction parameter α.According to eqs. ( 15) and ( 16), a negative value of α indicates energy transfer from the matter sector to the curvature sector, while a positive value of α implies the opposite, with energy flowing from curvature to matter.In fig.[4], we display the evolution of the parameter r mc for different α-values in model-A (left panel) and model-B (right panel).These computations use specific values of f (R)-model parameters, as mentioned in the caption of fig.[4].Peaks in the r mc -profiles signify the moments when the universe attains its highest matter density proportion.The trend of the peak height approaching unity for α ∼ −5 or even more negative values implies the presence of epochs with almost perfect coincidence between curvature and matter energy densities for significantly negative α values.For such negative α values, the late-time and distant-future accelerated phase of the universe experiences the presence of both matter and curvature energy densities in comparable proportions.Scenarios within the framework of the power law model (right panel of fig.[4]) may emerge for significantly negative α values where the matter energy density share may gradually rise as the universe progresses towards distant future epochs.However, the r mc values consistently remain below one, indicating that curvature density consistently exceeds matter energy density, even if by a slightly smaller margin during epochs when matter density contribution is at its highest.At earlier epochs during the decelerated phase of the universe, the negative α values result in very low r mc ∼ 0 in both models, indicating curvature dominance over matter during this period.For the case of α = 0 representing the standard f (R) gravity framework without any additional matter-curvature interactions, the profile of r mc exhibits a nearly symmetrical pattern centered around the epoch of maximum matter-energy density contribution, with the matter energy density remaining exceedingly minimal (r mc ∼ 0) during earlier as well as during distant future epochs.A positive α, on the other hand, results in a non-zero r mc at early stages, which subsequently decreases over time, reaching a minimum before resuming its ascent to reach its highest point during the epoch when the universe maximizes its matter density share.
The right panels of figs.[2] and [3] depict the evolution of deceleration and jerk parameters (q and j) for some benchmark values of the parameters corresponding to scenarios in model-A and model-B.These plots capture the kinematic aspects of cosmic evolution.During the early stages of evolution, they show positive values of both q and j, implying a decelerated phase with an increasing rate of deceleration.However, at a specific point, the deceleration parameter (q) rapidly declines, crosses the zero line, and becomes negative, signifying the onset of cosmic acceleration.It then sharply drops below the −1 threshold, rebounds, and ultimately stabilizes at −1 in the distant future.This decline of q and its eventual stabilization at q = −1 in the distant future are further characterized by a significant increase in the jerk parameter j, reaching its peak.This signifies a growing rate of acceleration, culminating in its maximum value.Subsequently, the jerk parameter starts to decrease, indicating a reduction in the rate of acceleration, and eventually settles at a positive value of unity in both model-A and model-B.All these intricate features of cosmic evolution result from the complex interplay among the interacting components of matter and curvature within the framework of modified gravity.

Conclusion
In this work, we explored the possibility for accommodating interactions between curvature and matter within the framework of viable f (R) gravity models, all the while ensuring that these models retain their ability to achieve cosmic acceleration.At the action level, we consider minimal coupling between matter and curvature.The radiation component is considered to be decoupled from both matter and curvature, implying the conservation of its energy-momentum tensor.We incorporate interactions between the matter and curvature components by introducing a source term Q in their respective continuity equations with opposite signs.This ensures energy conservation within the combined matter-curvature sector, enabling the potential for a continuous exchange of energy between the two sectors at a rate determined by Q.A specific form of the interaction term Q is chosen to illustrate a dynamic connection between the curvature and matter sectors by incorporating their individual energy densities in a multiplicative manner within Q, along with a dimensionless parameter (α) signifying the strength of coupling between the two sectors.The sign of α dictates the direction of energy flow between the two sectors.
An essential component of this study is the establishment of a 4-D dynamical framework that portrays the dynamics of curvature-matter interacting scenario.The introduction of the interaction term Q leads to alterations in the autonomous equations that govern the dynamical system, distinguishing it from those resulting from f (R) theories that do not consider matter-curvature interactions.We identified the critical points of the dynamical system and analyzed their stability, separately in the context of two viable f (R) models viz. the generalized Λ-CDM model (model-A) and the Power-law model (model-B) .The interaction parameter α assumes a crucial role, significantly shaping the dynamics of the system.Its impact on the stability characteristics of various fixed points within the dynamical system, as well as on the evolutionary dynamics of the universe during different phases, have been thoroughly investigated.
The incorporation of the curvature-matter coupling (α ̸ = 0) within the framework of f (R) gravity leads to an increase in the number of critical points, introducing new stable de-Sitter attractors.In both models, we notice the appearance of two stable de-Sitter attractors P 2 and P 4 , in contrast to the scenario without interactions (α = 0) where P 2 exists as the only de Sitter point.An in-depth analysis of the stability traits of these fixed points across an extensive range of α-values reveals intriguing patterns within the stability landscape of the matter-curvature interaction scenario.For α = −3, the stable de-sitter attractors P 2 and P 4 represents the same fixed point in both models.In model-A, for α > −3, P 2 takes the role of stable de-Sitter attractor, whereas in the range of −18 < α ≤ −3, it is P 4 that assumes this stabilizing role, provided the f (R)-model parameters b,c satisfy 1 < bc < 1.65.On the other hand, in Model B, stable de-Sitter acceleration is achievable for α > −3, with P 2 serving as the attractor and for α ≤ −3, P 4 emerges as the stable de Sitter attractor, operating within the viable range of the f (R) model parameter 0 < n < 1.As evident from the above discussion, a significant contrast between the two models is that, in Model B, the entire span of α permits stability of de-Sitter acceleration solution, whereas in Model A, this range is notably more restricted.Shifting our focus to identifying fixed points that signify stable acceleration of non-phantom nature (−1 < ω tot < −1/3), we observe that model-A presents two stable options, namely P 8A and P 9A , with their stability demonstrated for α > −4 and α ⩽ −4 respectively.Notably, these conditions are dependent on the parameters (b,c) of model-A, as illustrated in figs.[1](a) and [1](b).In contrast, within the viable range of n in Model B, no stable fixed points with non-phantom acceleration have been found for any value of α.
The evolution of various cosmological and cosmographic parameters have been explored for a nuanced assessment of the influences of curvature-driven dark energy on various phases of cosmic evolution within the framework of matter-curvature interactions.Examining the evolution of modified energy densities and the EoS parameter, particularly the matter-to-curvature density ratio (r mc ), along with two cosmographic parameters -deceleration (q) and jerk (j), offers valuable insights into potential evolutionary scenarios within scenarios involving the interplay of matter and curvature.In model-A, we observed that in both the early and late stages, curvature energy density exhibits significant dominance over matter energy density.During the transitional phase between these two phases, curvature density decreases while matter density increases over time, eventually reaching its peak.Subsequently, there is a resurgence in curvature's contribution and a decline in matter's contribution, ultimately reestablishing the pronounced dominance of curvature in the later stages.During the transition from early to later phases, the EoS parameter swiftly drops below −1/3, marking the onset of late-time cosmic acceleration and stabilizes around ∼ −0.8 in distant future, never reaching below −1, affirming the non-phantom nature of curvature-driven dark energy.This stable acceleration phase in the distant future is achieved with prominent curvature dominance, indicated by r mc ∼ 0. In model-B, while curvature density plays a dominant role in the early universe, matter density also becomes significant during the later phases, resulting in an order of magnitude around 1 (but always < 1) for r mc in the later stages.The EoS parameter drops from its positive value to values below −1/3 as the universe enters its accelerated phase.The descent persists, surpassing the non-phantom threshold of −1 before rebounding and stabilizing at this value in the later stages signifing the emergence of an unstable phantom region in model-B, albeit for a brief period.The evolution of parameters (q, j) remains consistent across both models.The positive jerk (j) and negative deceleration (q) observed during the late-time phase in both models confirm the presence of stable late-time cosmic acceleration in interacting curvature-matter scenario.We examined and illustrated (fig.[4]) the progression of r mc at specific reference values of the coupling parameter α.The parameter r mc is a measure of the extent to which matter dominates over the curvature component.In the absence of any interaction, the r mc profile remains symmetric around the epoch of maximum energy density contribution of the matter component, with r mc values tending to zero at both earlier and later time epochs in both models.When interactions are present (α ̸ = 0), positive values of α result in earlier epochs showing a notable contribution from the matter component.Conversely, with negative α values, the later epochs display a substantial contribution from the (dark) matter component.Our study is able to address the issue of cosmological coincidence in the late-time epoch with a negative α.
The ratio of matter density to curvature density, denoted as r mc and defined as Ωm Ωcurv , serves as a measure of the dominance of one component over the other in the cosmic landscape.r mc = 0 indicates complete curvature dominance, while r mc ≈ 1 suggests a coincidental balance between matter and curvature energy densities in the universe.The evolution of the parameter r mc is illustrated in fig.[2] (model-A) and fig.[3] (model-B), providing quantitative insights into the dominance of curvature over matter at different stages of cosmic evolution in the respective scenarios.The impact of the curvature-matter coupling strength becomes evident when comparing r mc profiles computed with different values of the interaction parameter α.Within the framework of equations ( 15) and ( 16), a negative α implies energy transfer from matter to curvature, while a positive α signifies the opposite, with energy flowing from curvature to matter.fig.[4] displays the evolution of r mc for different α values in model-A (left panel) and model-B (right panel), using specific values of f (R)-model parameters as mentioned in the caption of fig.[4].Peaks in the r mc profiles represent moments when the universe attains its highest matter density proportion.The trend of peak height approaching unity for α ∼ −5 or more negative values indicates epochs with nearly perfect coincidence between curvature and matter energy densities.For significantly negative α values, the late-time and distant-future accelerated phase of the universe experiences both matter and curvature energy densities in comparable proportions.However, r mc values consistently remain below one, suggesting that curvature density consistently exceeds matter energy density, even during epochs with the highest matter density contribution.At earlier epochs during the decelerated phase of the universe, negative α values result in very low r mc , indicating curvature dominance over matter.In the case of α = 0, representing the standard f (R) gravity framework without additional matter-curvature interactions, the r mc profile exhibits a nearly symmetrical pattern centered around the epoch of maximum matter-energy density contribution.A positive α results in a non-zero r mc at early stages, decreasing over time and reaching a minimum before resuming its ascent to reach the highest point during the epoch when the universe maximizes its matter density share.
While numerous viable f (R) models exist to address the enigma of dark energy in cosmology, the qualitative insights from our study can form a foundation for more intricate examinations.These findings serve as supplementary constraints when assessing the viability of f (R) models, alongside the broader range of cosmological constraints applied to f (R) theories.Each f (R) gravity model leaves its distinct imprint on the evolution of cosmological perturbations.Specifically, the scrutiny of structure formation in the universe proves to be highly responsive to the nature of interaction between dark energy and dark matter within the framework of f (R) gravity.Notably, the interaction between curvature driven dark energy and dark matter in the early stages of the universe may potentially influence the epoch of matter-radiation equality.The manifestation of anisotropies resulting from this phenomenon can be quantitatively assessed in terms of the growth of structure formation within the curvature-matter interacting framework of f (R) gravity.Interestingly, the interaction between f (R)-dark energy and dark matter during the early phase of cosmological evolution may offer an explanation for the disparity in the Hubble tension parameter observed between local measurements and those derived from cosmic microwave background observations.Furthermore, the study of matter perturbations and the local gravity approximations of interacting dark energy in f (R) gravity holds promise for shedding light on the inconsistent predictions between non-interacting and interacting f (R) dark energy models.Additionally, exploring alternative models of interacting dark energy within the framework of f (R) gravity holds the potential to extract further physical implications and discern their cosmological consequences.
In summary, this research deepens our understanding of the intricate dynamics in modified f (R) gravity theories, particularly in the presence of interactions between curvature and matter.The emergence of stable de-Sitter attractors, nuanced stability characteristics, and distinctive evolutionary patterns in this scenario are instrumental in elucidating the various facets of the complex behavior of the system.The findings presented here not only enhance our understanding of the dynamics of these systems but also open up new avenues for further exploration in the field of cosmology, paving the way for gaining valuable insights into the fundamental dynamics of the universe.
µν , where G µν is the Einstein tensor and T µν is determined by f (R) and its derivatives, and vanishes for f (R) = R reducing to Einstein gravity.T (M) µν is a modification obtained by scaling the usual energy-momentum tensor Tµν from Einstein's equation by the factor [df (R)/dR] −1 .T (curv) µν and (II) ∇ µ T (m) µν = −∇ µ T (curv) µν ≡ −Qν ̸ = 0 which collectively, validate the overall conservation of the comprehensive tensor T (tot)

Figure 4 :
Figure 4: Evolution plot of the parameter r mc for different values of coupling strength α.Left Panel: For f (R) model with generalised Λ-CDM model with parameter bc = −3.Right Panel: For power law form of f (R) with value of model parameter n = 0.9.
tot ≡ P tot /ρ tot represents the 'grand' equation of state (EoS) parameter of the comprehensive (radiation+matter+curvature) fluid.Though the combined stress-energy tensor T