1. Introduction

The accelerating expansion of the universe recognized continuously in the recent past is supported by high precision observational and experimental data such as Supernovae Ia [1], cosmic microwave background (CMB) temperature anisotropies observed with the help of WMAP [2] together with high redshift supernovae surveys [3], observations of large-scale structure (LSS) [4], and baryon acoustic oscillations [5]. To account for this mystical cosmic outlook, two different classes of models are approached. The first kind of model assumes the presence of some sort of exotic energy component having negative pressure yet possess positive density, dubbed as dark energy (DE), which impelled this accelerated expansion. However, finding the origin of such DE models was a major problem. Among many others, the simplest DE candidate is the cosmological constant, but it also encounters difficulties like the fine-tuning problem and the coincidence problem [6, 7].

On the other hand, modified theories of gravity are introduced, claiming to explain the accelerated expansion while dealing with spacetime geometry. The approach of such gravities is the modification of general theory of relativity (GR) at low curvatures or at larger scales by reforming gravity part of the Einstein–Hilbert action. A broad spectrum of such modified gravity models is investigated including scalar–tensor theory [8], f(R) gravity [9], Gauss–Bonnet theory [10], f(R, T) theory [11], f(R, G) theory [12], braneworld models [13], TeVeS (tensor–vector–scalar) [14], and Galileon gravity [15]. The f(R) theory of gravity is the most elegant and propitious modified gravity among others, where the Ricci scalar R in the curvature part of Hilbert’s action is replaced by differentiable term f(R) without addition of any new kind of matter.

Initially, f(R) actions were carried out by Weyl and Eddington [16, 17]. Models based on f(R) gravity are expected to successfully develop early time inflation as well as late time acceleration. Nojiri and Odintsov [18] utilized positive as well as negative powers of curvature term R to yield the early inflationary and late time acceleration phases. Amendola et al. [19] deduced prerequisites for the f(R) models to be cosmologically viable.

To derive Einstein’s field equations, variation of the action takes place either through metric or the Palatini formalism [20, 21]. For the former, the variation is performed considering the metric to be independent and for later both metric and connection as independent variables. The field equations for both the formalisms are alike regarding the linear Lagrangian in R; however, the scenario is different for nonlinear f(R). The concept of these formalisms is comprehensively explained in [22, 23].

Starobinsky [24] introduced very first f(R) model of inflation, i.e., f(R) = R + αR², where the term αR² is responsible for accelerated expansion of universe. This inflationary model is compatible with the observed temperature anisotropies in CMB and a strong variant to scalar field models of inflation. On the other hand, the f(R) models with negative power of R, for example, f(R) = R–αR⁻ⁿ, (where α > 0, n > 0) were suggested to explain the late-time acceleration in metric formalism. The condition, f_RR < 0, gives rise to instability and as a result local gravity constraint is not satisfied. Also, the models of the form 1/R are found to be unstable in this formalism [25]. Due to the large coupling between DE as well as dark matter, such kinds of models do not possess matter-dominated epoch. The Palatini formalism seems attractive since the field equations are of second-order partial differential equations (PDEs), unlike the metric approach which yields fourth-order PDEs. Moreover, it not only satisfies the solar system test but also provides correct Newtonian limit. It is also observed that the restrictions like f_RR > 0 are not required in Palatini formalism to maintain stability [26–28]. In Palatini formalism, the model of the form f(R) = R − μR⁻ⁿ is free from instability in the presence of matter [29].

The field equations obtained for modified theories of gravity comprise huge nonlinear terms which complicate their analytic or numerical solutions, ultimately making their comparison with observational data troublesome. Therefore, various methods are devised to tackle such equations among which one technique is dynamical system analysis. Various applications of this analysis in context of cosmology are discussed by Wainwright and Ellis [30] and Coley [31].

In cosmology, the scheme of this analysis is to determine numerical solutions of a system which in turns help to figure out the qualitative behavior of the respective cosmological model or a system. The first step in this technique is to introduce dimensionless variables that transform nonlinear field equations into the ordinary differential equations (ODEs) and the second is to analyze the critical points from these ODEs. Afterward, the eigenvalues are determined using the Jacobian matrix for the corresponding fixed points. In Palatini formalism, solving the system of two autonomous equations for eigenvalues is very handy as it deals with a 2 × 2 Jacobian matrix unlike the metric formalism where difficulty arises while solving fourth-order matrix. Subsequently, the stability analysis is carried out regarding these eigenvalues.

In a two-dimensional system, if both the eigenvalues are unequal, real and positive, then the concerned fixed point is attractor and the trajectories are unstable. For distinct real and negative eigenvalues, the fixed point is a repeller and trajectories followed are stable, whereas if both eigenvalues are of opposite sign, then the fixed point is saddle, showing that some of the trajectories are repelled and others attracted. Considering the case where one of the two eigenvalues is zero and the other is positive, the fixed points are unstable, whereas if the other eigenvalue is negative, then the stability is inconclusive or fails to determine the stability. This methodology of dynamical system analysis is widely used, especially in framework of f(R) gravity. Samart, Silasan, and Channuie [32] examined the dynamics of interacting dark matter and DE given by Q = 3αHρ_m, in context of f(R) gravity using dynamical system analysis. They discovered that the presence of the interacting terms as well as the pertinent model parameters altered the dynamical profiles of the universe in the f(R) DE models. Lobato et al. [33] introduced new variables to study the dynamics of exponential gravity f(R) model using Palatini approach and obtained a fixed point which behaves like attractor. They also found multiple values of effective equation of state (EoS) for the different values of α. Shah, Samanta, and Capozziello [34] straightforwardly used Einstein gravity to analyze the cosmological dynamics for different mixtures of cosmic fluids. Roy [35] interpreted the dynamical analysis of different DE models including quintessence and chameleon models in detail. Recently, Shabani and Farhoudi [36] investigated the dynamics of f(R, T) cosmological models considering three different Lagrangian couplings and found that among them, the pure minimal coupling lacks matter epoch.

This paper is outlined as follows. Section 2 briefly describes the f(R) field equations in Palatini formalism and autonomous equations are derived using dimensionless variables. In this section, the general formalism for finding fixed points and eigenvalues is explained briefly. Section 3 comprises the AGPT f(R) gravity model to discuss the cosmological evolution as well as different dynamical phases of universe. Section 4 contains the concluding remarks concerning the obtained results.

2. The Autonomous Equations for Palatini f(R) Formalism

The cosmological background of universe in reference to the Einstein’s modified equation of motion could be derived from f(R) theory. The modification of the Einstein–Hilbert action in case of Palatini f(R) gravity is formulated as

()

Bear in mind that L (satisfying the conservation equations) is the Lagrangian density, where the subscripts m and r elaborate the dust-like matter as well as the radiation. Note that κ = 8πG, G simply being the gravitational constant. As mentioned above, Palatini formalism deals with metric and affine connection independently while varying the action. So, variation of equation (1) with respect to the metric tensor g_μν gives

()

where F = (∂f/∂R), R_μν is Ricci tensor, and T_μν is the energy momentum tensor regarding the matter field defined as

()

The variation regarding the connection is given as

()

where the covariant derivative  is defined with respect to the independent connection . Variation of equation (1) with respect to the independent connection along with the contractions over µ and ν results in

()

To comprehend the cosmological dynamics, we derive the field equations considering the universe to be flat Friedmann–Lema  tre–Robertson–Walker (FLRW), defined by the spacetime

()

where the scale factor a(t) is dependent upon time. Taking into consideration nonrelativistic matter and radiation, the generalized Friedmann equations are given as

()

()

where  is the Hubble parameter, dot represents differentiation with respect to time t, and ρ_m and ρ_r are energy densities satisfying continuity equations

()

()

Note that equation (8) is obtained by contracting equation (2). By connecting equations (7) and (8) with the continuity equations (9) and (10), we have

()

()

where prime denotes derivatives with respect to R and

()

We introduce dimensionless variables in order to derive autonomous equations to study the dynamical aspects of the universe under f(R) theories

()

Such a normalized set of variables is chosen to obtain a bounded dynamical system keeping in mind that the normalized variables are well behaved and possess direct physical consequence. For these dimensionless variables y₁ and y₂, constraint equation (12), which transforms to standard Friedmann equation H² = κ(ρ_m + ρ_r)/3 for f(R) = R, implies

()

Taking derivative of equation (12) and utilizing equations (7)–(15) yields

()

Using the above two equations will lead to the following evolution equations for variables y₁ and y₂:

()

()

where

()

Also, the following constraint equation holds, depicting that R as well as C(R) could be written in terms of variables y₁ and y₂.

()

Equations (17) and (18) clearly depict the dependence of variables y₁ and y₂ upon C(R) since y₁ and y₂ cease to attain equilibrium whenever C(R) diverges. For the discussion of the dynamics of cosmological systems, equilibrium, fixed points, and stabilities play a key role, so it is important to suppose that C(R) is well-behaved. The curves (y₁ (N), y₂(N)) describe the system solutions on the phase space of (y₁, y₂). We put equations (17) and (18) equal to zero to study the fixed point, i.e., (dy₁/ dN) = 0, (dy₂/ dN) = 0. Solutions closer to the critical points may be attracted or repulsed depending on their nature.

The effective EoS relates the total pressure P_tot to the total energy density ρ_tot of the cosmic fluids. In Palatini f(R) gravity, the modification to the gravitational field equations can have an impact on the evolution of the cosmic fluid such as radiation, DE, or dark matter, and consequently, the effective EoS parameter ω_eff may evolve differently compared to GR. Since the Hubble parameter is related to ω_eff as , equation (16) redefines the effective EoS parameter as

()

The specific behavior of ω_eff in Palatini f(R) gravity depends on the chosen function f(R), the cosmological scenario, and the dynamics of the universe. The above obtained fixed points will be helpful to easily evaluate the corresponding EoS using expression in equation (21).

3. Dynamical Analysis for the Cosmological Model

Dynamical system analysis is a qualitative method which offers to transform Einstein field equations into autonomous system of first ODEs as far as spatially homogenous cosmology is concerned. This section is devoted to carrying out the stability analysis for our model of interest used frequently in the literature utilizing the above formalism. We consider the f(R) model named as Amendola Gannouji Polarski Tsujikawa (AGPT) model [37], given as

()

with μ > 0, R_c > 0, being a small positive constant and 0 < p < 1. This specific form of f(R) gravity model avoids singularities since at lower curvatures, the correction term – will be negligible. Due to the negative sign in the correction term – of AGPT model, it mimics DE component which is responsible for accelerated expansion of universe. Moreover, for the chosen form of f(R) gravity model, since the value of p ranges between 0 and 1, the correction term  grows slower than R and thus mild deviations from GR are observed as compared to the polynomial forms of f(R) gravity. Prior to this study, the cosmological dynamics of this model is studied using the metric approach where the absence/presence of interactions between components of universe, e.g., dark matter, DE, and radiation, was assumed [38]. Examining the AGPT model using the Palatini formalism offers valuable understanding of the larger category of f(R) gravity theories. It offers a platform for investigating and understanding how alterations to general relativity can impact cosmic dynamics. Since p is restricted between 0 and 1, this choice helps in achieving the desired acceleration without the need for excessive fine-tuning. Equation (22) helps in explaining all the phases of universe, especially the early inflation (followed by radiation and matter domination) as well as late time acceleration considering no interaction between nonrelativistic dark matter (dust) and the radiation.

For the given model, equation (19) takes the following form:

()

Ultimately, constraint equation (20) gives

()

showing R in terms of dimensionless variables y₁ and y₂. This empowers C(R) in equation (23) to be function of y₁ and y₂, given as

()

A graphical representation of all the above eigenvalues is given in Figure 1 and relative properties are highlighted in Table 1.

Table 1. Critical points, eigenvalues, and nature of trajectories for .

( y1,  y2)	(λ1, λ2)	ωeff	Trajectory type	Ωr	Ωm	Ωgeo
Pm = (0, 0)	(3(1 − p), −1)	0	Saddle	0	1	0
Pr1 = (0, 1)	(4 − 3p, 1)	1/3	Unstable	1	0	0
Pr2 = (0, 1)	(1, 1)	−2/3 + 1/p	Unstable	1	0	0
Pd = (1, 0)	(−3, −4)	−1	Stable	0	0	1
P = (−(p − 1/p − 3), 0)	(0, 1)	1 − 1/p	—	0	(2(p − 1)/p − 3)	−(p − 1/p − 3)

The above stability analysis can also be complemented through the phase space portrait. Figure 1 provides the trajectory graph for system equations (17) and (18) in order to visually comprehend the behavior of trajectories in the vicinity to given fixed points. It is acceptable to use a graphical range between −1 and 1, as all the critical points are between 0 and 1. The fixed point (0, 0) is a saddle point because it is clearly observed that the trajectories that are drawn to it are being repulsed before they reach it. Eventually, every trajectory is drawn to (1, 0) and moves away from (0, 1), indicating that the critical point (1, 0) is stable and (0, 1) is the unstable one. Hence, the result obtained through phase portrait complements those that are obtained through eigenvalue analysis. It is to be mentioned here that the Mathematica and Maple software programs are used for plotting eigenvalues and trajectories.

3.1. Evolution of the Universe

This section is devoted to studying the dynamical evolution of the universe. From constraint equation (25), we have already seen that C(R) and R can be written in terms of dimensionless variables y₁ and y₂. The same concept can be utilized for transforming the effective EoS into a rather dimensionless expression by writing equation (21) in terms of y₁ and y₂. This mechanism allows to efficiently analyze the behavior of our dynamical system for multiple scales and under different conditions. It is favorable to separately evaluate the terms of equation (21). Firstly, as far as the systematic analysis for evolution of ω_eff is concerned, for our model, the third term of equation (21) becomes

()

From equation (13), variable ξ takes the form

()

Using this value of ξ in last two expressions of equation (21), deduce the following, respectively,

()

and

()

where the factor  appearing in above equations attains the value

()

Combining the above determined terms with equation (21), we have

()

Therefore, equation (31) gets the form

()

Figure 2 shows the evolution of dimensionless entities y₁ and  y₂ as well as the EoS parameter ω_eff where p = 0.5. The initial conditions are chosen as y₁ = 10⁻³⁰, y₂ = 1 − 3 × 10⁻⁶. The initial conditions are chosen due to the fact that smaller the ratio (y₁(0)/y₂(0)), longer will be the matter-dominated phase.

It can be clearly observed that the radiation domination ω_eff = 1/3, matter domination ω_eff = 0, and de-Sitter acceleration ω_eff = −1 epochs occur in a sequence at P_r1, P_m, and P_d, respectively. Hence, three out of the four eras of universe are possible for our dynamical system. Moreover, in Palatini gravity, different fractional energy density components are defined as

()

where Ω_m, Ω_r, and Ω_geo represent, respectively, the matter, radiation, and geometric components of these densities. Dynamical analysis encourages transformation of these fractional energy densities to dimensionless entities. In terms of y₁ and y₂, equation (33) gives

()

Figure 3, depicts that at the initial stage, universe is filled with radiation and the respective density is maximum as well at that point. As ω_eff started changing from 1/3 toward 0, the radiation density reduces and matter density prevails. The universe remained dominated with matter as long as ω_eff = 0. When the geometric component starts to dominate, accelerated expansion started as ω_eff went negative and abruptly reached −1, showing the de-Sitter phase. Hence, for the AGPT model, except the inflationary phase, the sequence of all the other three phases is achievable.

4. Concluding Remarks

Even though GR has been incredibly successful at explaining gravitational interactions at astronomical scales, it is non-renormalizable and hence cannot be quantized [39, 40]; therefore, the search for a more thorough theory of gravity has persisted. One such alternative theory that has attracted interest in the field of theoretical physics is f(R) gravity, which naturally serves to explain the accelerated expansion of the universe, unlike GR, where modeling of DE component is required for explaining this acceleration. Various models of f(R) theories provide unification of different phases of universe like the inflation and late time acceleration [41]. Such kind of unification was not possible in GR. These aspects motivate to study f(R) theories as well as other modified theories over GR for studying the universe.

Since the field equations in the Palatini formalisms are second-order PDEs rather than the fourth-order PDEs produced by the metric approach, they appear more appealing among the two formalisms of f(R) gravity. Additionally, it passes the solar system test, provides the correct Newtonian limit, and does not require constraints like f_RR > 0 to guarantee stability.

Prior to our study, a fixed-point analysis is performed in [42, 43] for the Starobinsky model, Hu–Sawicki model, and Gogoi–Goswami model in metric formalism; however, these models only explained the transitions from matter-dominated to DE-dominated universe (accelerated universe). Previously, the dynamics of models of the form f(R) = R + αR^m − βR⁻ⁿ, f(R) = R + αR² and f(R) = R − βR⁻ⁿ were studied in [44]. The sequence of three different phases, i.e, radiation dominated, matter dominated, and late time acceleration, was obtained for R − βR⁻ⁿ. However, f(R) = R + αR² was not able to depict such a sequence. Although the inflationary phase was visible in f(R) = R + αR^m − βR⁻ⁿ, it started with a nonstandard radiation phase.

In the present work, we have discussed the evolution and stability properties of specific dynamical system in Palatini f(R) gravity. For this purpose, we have established a qualitative technique dubbed as dynamical stability analysis in the background of FLRW universe. This technique followed the transformation of nonlinear DE to autonomous equations (equations (17) and (18)) by introducing dimensionless variables y₁ and y₂. With the aid of constraint equation (20), C(R) is declared to be function of dimensionless variables y₁ and y₂. For the AGPT f(R) model, i.e., , with μ > 0, R_c > 0, and 0  <  p  <  1 [31], four different critical points P_m  :  (0, 0),  P_r : (0, 1),  P_d  :  (1, 0), and P : ((p − 1/p − 3), 0) were obtained. The respective eigenvalues characterized P_m to be a saddle point in the range 0  <  p  <  1, P_r to be a repeller, and P_d an attractor.

The phase space plot in Figure 4 also shows that the point (0, 0) owns a saddle nature, since trajectories that are attracted toward origin are repelled before reaching it. Moreover, the unstable nature of the trajectories at point (0, 1) designates it to be a repeller. It can also be observed from Figure 4 that all the trajectories are attracted toward the fixed point (1, 0), pointing it to be the stable node. We have found that the restriction on p caused the critical point P to not lie in the asymptotic regime. The negative value of fixed point, y₁ < 0 contradicts the physical range of dimensionless variable (y₁, y₂) ∈ [0, 1]. Following the stability analysis, we have analyzed the evolution of the universe. In this context, we have utilized the EoS parameter, which plays a pivotal role in modeling different phases of accelerating universe. Its effectiveness depends on the appropriate selection of parameters for the specified model. Owing to the fact that R shows dependence on dimensionless variables, the EoS parameter ω_eff can also be featured as dimensionless quantity in Palatini gravity. Figure 2 represents the evolution of ω_eff as well as y₁ and y₂, where the initial conditions for y₁(0) and y₂(0) are chosen in such a way that their ratio (y₁(0)/y₂(0)) describes how long the matter-dominated phase could last.

Figure 2 clearly depicts the possibility of the radiation-dominated phase ω_eff = 1/3 at P_r = (0, 1), matter-dominated phase (ω_eff = 0) at P_m  :  (0, 0), and de-Sitter epoch (ω_eff = −1) at P_d  :  (1, 0).

We have observed that the EoS parameter varied directly with y₂ until the matter-dominated phase approached. Moreover, when variable y₁ varies from 0 to 1, EoS behaves inversely with respect to y₁ and attained the de-Sitter phase at (1, 0). We have found that these results correspond to f(R) = R − βR⁻ⁿ [44]. Thus, for the AGPT model, three among the four eras of universe evolution can be obtained unlike metric formalism where such a sequence of phases was not obtained. Also, our results differ from those obtained for f(R) = R + αR² and f(R) = R + αR^m − βR⁻ⁿ. From the analysis of Figure 3, we have seen that energy densities as well as the EoS evolution for our f(R) model are in accordance with those of Λ CDM model as discussed in [35]. Moreover, a sequence of these three phases is obtained for Hu–Sawicki f(R) gravity for Bianchi I and V cosmologies [45].

It is concluded that while the Palatini formalism provides more convenient and robust theoretical framework which aligns well with the observations as compared to metric formalism, it also encounters some challenges as far as dynamical analysis is concerned. For complex forms of f(R) gravity models, the equations of motions obtained are highly nonlinear and complex, which makes normalization as well as conducting a thorough dynamical system analysis difficult. Also, the equivalence of Palatini formalism to a scalar–tensor theory with an additional scalar field introduces additional degrees of freedom that need to be tackled while deriving equations of motions [27, 46]. However, these limitations do not invalidate the efficacy of Palatini formalism. Incorporating the results with other formalisms and gravity theories can strengthen the importance of this theory. From the above discussion, we affirm that the current study may play a useful role in probing the Palatini f(R) gravity and its extended versions.