1 INTRODUCTION↩︎

The universe is currently experiencing an epoch of accelerated expansion, a phenomenon that has been documented over the past century through various observational efforts. Among the recent key sources of evidence are the measurements of luminosity distance and redshift of Type Ia supernovae [1]–[3], which have played a central role in this discovery. Other important observations, such as the cosmic microwave background radiation [4], baryon acoustic oscillations [5], [6], and the determination of the Hubble constant [7], have further contributed to our understanding of this accelerated phase. Due to sheer lack of any information about what’s driving this late-time expansion, cosmologists have attributed it to the mysterious dark energy which exerts a large negative pressure to counteract gravitational collapse.

Although many theoretical models have been proposed, the true nature of dark energy is still unknown. One of the first and most widely studied models is the \(\Lambda\)CDM model, which includes the cosmological constant (\(\Lambda\)) and cold dark matter. This model has been popular, especially among the particle physics community. However, the severe fine-tuning associated with the observed value of \(\Lambda\), along with the coincidence problem [8], [9], suggests that a dynamical mechanism may be responsible for cosmic acceleration. This has driven interest in exploring dynamical dark energy models based on scalar fields, where a scalar field is responsible for generating negative pressure, leading to the accelerated expansion of the universe at late times. Depending on the form of the Lagrangian, there are mainly two types of scalar field models studied so far: the quintessence models [10]–[25] and the \(k\)-essence models [26]–[34].

Another promising direction involves modifying the gravitational sector itself. Among such theories, \(f(R)\) gravity, which generalizes the Einstein-Hilbert action by promoting the Ricci scalar \(R\) to a function \(f(R)\), has attracted considerable interest [35]–[40], sotiriou2010?, defelice2010?. These models can naturally drive late-time acceleration without introducing exotic matter components. However, pure \(f(R)\) models have their limitations. In particular, they often face challenges in fitting both cosmological and local gravitational tests simultaneously, and they may lack the flexibility needed to reproduce a wide range of cosmological behaviors [41].

To overcome these shortcomings, an extended class of theories has been proposed in which a scalar field \(\phi\) is explicitly coupled to the curvature, leading to models of the form \(f(R, \phi)\) [42]–[46]. These models merge the advantages of both scalar field dynamics and modified gravity, offering a richer framework to explore dark energy and late-time acceleration. More recently, certain classes of \(f(R,\phi)\) inflation models have also been shown to produce a slightly enhanced scalar spectral index corresponding to the recent results from Atacama Cosmology Telescope (ACT). The latest ACT data release has disfavored long-standing inflationary paradigms such as Higgs inflation and Starobinsky inflation, though there have been efforts to show that modified dynamics may be able to restore these models to within the acceptable range [47]–[55]. More importantly, it has been highlighted by some authors that the non-minimal coupling \(\phi R\) is well-suited to produce the required enhancement in the scalar spectral index.

An important subtlety in studying \(f(R,\phi)\) models lies in the choice of variational principle. In the metric formalism, the action is varied with respect to the metric, and the connection is taken to be Levi Civita. In contrast, the Palatini formalism treats the metric and affine connection as independent variables [56]–[58]. The resulting equations of motion are second order in derivatives, thereby avoiding the instabilities that typically emerge in higher derivative metric gravity theories. Furthermore, the Palatini approach often leads to new scalar tensor equivalences that are not present in the metric formulation, making it a compelling framework for exploring modified gravity.

A standard technique for revealing the effective degrees of freedom is the application of a Weyl transformation which maps the Jordan frame action into its equivalent Einstein frame counterpart (wherein gravitational sector resembles GR and the modifications are absorbed into scalar field terms with non standard kinetic and potential structures) [59], [60]. In the case of \(f(R,\phi)\) models in Palatini formalism, the transformed action exhibits a particularly novel structure: the scalar sector acquires a higher order kinetic interaction term and an effective potential [61]. The appearance of higher order kinetic terms places the model within the broad class of \(k\)-essence theories [28], [62]. These models are known to produce rich cosmological behavior, including dynamical dark energy with varying equation of state, attractor solutions, and even phantom like regimes without introducing ghosts [63]. In this paper, we investigate the cosmological behavior of a class of Palatini \(f(R, \phi)\) gravity models in the Einstein frame and examine their viability as an explanation for the observed acceleration of the universe. Our goal is to determine whether such models can provide an improved or complementary description to \(\Lambda\)CDM while remaining consistent with current observational data. We perform a dynamical system analysis of cosmological evolution which provides a systematic way to classify possible asymptotic states of the Universe by converting the cosmological equations into an autonomous system of differential equations [19], [64]. Our main goal is to identify conditions under which such models admit stable accelerated attractors at late times and thus providing a framework for dynamical dark energy. In addition, we highlight the role of the induced quartic kinetic interaction term in driving acceleration.

The paper is organized as follows: In Section 2, we present the theoretical framework for Jordan frame \(f(R,\phi)\) gravity theories such that they could effectively describe a late-time accelerated epoch. In Section 3, we review the corresponding Einstein frame action in metric formalism and provide a generalized canonicalization approach without fixing specific forms of potentials and coupling functions. Next, in Section 4, we briefly look at the Einstein frame action in Palatini formalism and choose a generalized scalar-Starobinsky model (based on consistency conditions). Its dynamical stability is analyzed in Section 5 where we list fixed points that could provide an efficient candidate for dark energy at late-times. Constraints and other bounds for all the relevant parameters are found with respect to late-time observational data in Section 6. Finally, Section 7 summarizes our conclusions and outlines possible avenues for future research.

2 Jordan Frame \(f(R,\phi)\) theories: Early and Late-time Acceleration↩︎

We consider the following action in the Jordan frame, \[\label{eq:jordan-action} S=\int d^{4}x\sqrt{-g}\left[\frac{1}{2\kappa}f(R,\phi)-\frac{1}{2}\nabla{_\mu}\phi\nabla{^\mu}\phi-V(\phi)\right].\tag{1}\] Here \(g\) is the metric determinant, \(R\) is the Ricci scalar for our given spacetime, and \(\kappa=M_{\rm Pl}^{-2}\) where \(M_{\rm Pl}\) is the reduced Planck mass. Deriving the corresponding Einstein field equations and the equation of motion for the scalar field \(\phi\), we obtain, \[\begin{align} &FR_{\mu\nu}-\frac{1}{2}fg_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}F+g_{\mu\nu}\square F=\kappa\left[(\nabla_{\mu}\phi\nabla_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\nabla^{\alpha}\phi\nabla_{\alpha}\phi)-Vg_{\mu\nu}\right],\tag{2}\\ &\square\phi+\frac{1}{2}\left(\frac{f_{,\phi}}{\kappa^2}-2V_{,\phi}\right)=0,\tag{3} \end{align}\] where \(F=\frac{\partial f}{\partial R}\). Now, for the background spacetime, we consider a spatially flat FLRW metric given by the line element, \[\label{eq:flrw-metric} ds^2=-dt^2+a^2(t)(dx^2+dy^2+dz^2),\tag{4}\] where \(a(t)\) denotes the scale factor of expansion. Assuming homogeneity in the background, i.e. \(\phi\equiv\phi(t)\), 00 and \(ii\) components of the Einstein field equations 2 can be expressed respectively as, \[\begin{align} 6F(\dot{H}+H^2)+(2\kappa V-f)-6\dot{F}H+\kappa{\dot{\phi}}^2&=0,\tag{5}\\ -\kappa{\dot{\phi}}^2+2F(\dot{H}+H^2)+4FH^2+(2\kappa V-f)-4H\dot{F}-2\ddot{F}&=0.\tag{6} \end{align}\] Further, the equation of motion of \(\phi\) given by 3 now becomes: \[\label{eq:frphi-equation-phi} \Ddot{\phi}+3H\dot{\phi}-\frac{1}{2}\left(\frac{f_{,\phi}}{\kappa}-2V_{,\phi}\right)=0.\tag{7}\] These models have been studied extensively in reference to the inflationary paradigm [45], [58], [65]–[68]. However, as discussed earlier, the allure of modified gravity and scalar-tensor theories (STTs) lies in their ability to explain a variety of phenomena ranging from inflation to late-time acceleration. To demonstrate this, we shall try to find conditions for a general \(f(R,\phi)\) theory of the form 1 that could provide a stable dark energy candidate. In general \(f(R,\phi)\) theories, both \(\phi\) and \(F\) are important scalar degrees of freedom (DOF) and either can be responsible for late-time expansion. We can find conditions for \(\phi\) to support late-time acceleration by rewriting 7 as, \[\Box\phi=V_{,\phi}-\frac{f_{,\phi}}{2\kappa}\equiv\frac{\partial U_{\rm eff}}{\partial\phi},\] Now, we can find the condition of existence of an extremum as, \[\begin{align} \frac{\partial U_{\rm eff}}{\partial\phi}&=0,\nonumber\\ f_{,\phi}&=2\kappa V_{,\phi}. \end{align}\] And for this extremum to be a stable minimum, we find, \[\begin{align} \frac{\partial^2 U_{\rm eff}}{\partial\phi^2}&>0,\nonumber\\ 2\kappa V_{,\phi\phi}&>f_{,\phi\phi}. \end{align}\] It is clear that we can go no further without assuming explicit forms for \(f(R,\phi)\) and \(V(\phi)\). But as this goes against the foundation of this work, we stop here and shift our focus on the other viable scalar DOF: \(F\). The trace equation corresponding to 2 , without the assumption of homogeneous background, is given by, \[3\Box{F}+FR-2f=\kappa T=-\kappa[\partial_{\mu}{\phi}\partial^{\mu}{\phi}+4V(\phi)],\] where \(T\) represents trace of the total stress-energy tensor that includes the scalar field \(\phi\) and other non-gravitational components (that are no longer decoupled from the system at late times). We can write rewrite the trace equation as [69], \[\Box{F}=\frac{1}{3}(\kappa T+2f-FR)\equiv\frac{\partial{V_{\rm eff}}}{\partial{F}}.\] Now, for the existence of a local extremum, we should have, \[\begin{align} \frac{\partial{V_{\rm eff}}}{\partial{F}}=&\frac{1}{3}(\kappa T+2f-FR)=0,\nonumber\\ \implies\;\;\;\;\;R=&\frac{1}{F}(2f+\kappa T). \end{align}\] Assuming a de Sitter expansion at late times, \(R\approx 12H^2>0\). Thus we obtain, \[2f+\kappa T>0.\] For this extremum to be a local minimum, i.e. for a stable late-time cosmic expansion, we require \[\begin{align} \label{eq:condition-stable-latetime-expansion} &\frac{\partial^2{V_{\rm eff}}}{\partial{F^2}}=\frac{1}{3}\left(2\frac{\partial{f}}{\partial{F}}-R-F\frac{\partial R}{\partial F}\right)>0,\nonumber\\ \implies &2\frac{\partial{f}}{\partial{\phi}}\left(\frac{\partial{F}}{\partial{\phi}}\right)^{-1}+F\left(\frac{\partial{F}}{\partial{R}}\right)^{-1}-R>0, \end{align}\tag{8}\] which is true only when \(\frac{\partial{F}}{\partial{\phi}},\frac{\partial{f}}{\partial{R}}\neq0\). Similar condition for \(f(R)\) gravity has been derived for a stable inflationary era in [69]. This condition can be further simplified assuming that \(f(R,\phi)\) is factorizable as \(f=p(R)q(\phi)\), for some functions \(p(R)\) and \(q(\phi)\). To do so, we list below a few properties of \(f(R,\phi)\) and their consequences for \(p(R)\) and \(q(\phi)\):

• \(f(R,\phi)\) has 2 mass dimensions (see 1 ). Consequently, we choose \(p(R)\) to have 2 mass dimensions while \(q(\phi)\) remains dimensionless.

• \(f(R,\phi)\) must be analytic in both \(R\) and \(\phi\). This is required because, as \(\phi\to0\), \(f(R,\phi)\) must become some function of \(R\). \(R\to0\) must also be a well defined limit in order to explain transitions between various epochs of the universe. In factorizable scenarios, analyticity of \(p(R)\) and \(q(\phi)\) implies that as \(\phi\to0\), \(q(\phi)\) returns to some finite value, i.e. there exists a non-zero minimum of \(q\).

• In weak gravity regimes, \(f(R,\phi)\) must reduce to the Einstein-Hilbert action coupled with some function of \(\phi\), i.e. \(p(R)\to R\).

The condition 8 can then be expressed as, \[\label{eq:condition-stable-latetime-expansion-pq} \frac{2p}{p_R} + \frac{p_{R}}{p_{RR}}-R>0\tag{9}\] where the subscript \(R\) represents a differentiation w.r.t. \(R\). By assuming factorizability of \(f\), we arrive at the weaker condition 9 which is valid as long as \(p_R,\, p_{RR}\neq0\). Note that this condition is completely independent of the scalar coupling \(q(\phi)\). For \(p=R+\alpha R^2\), one can verify that the condition is satisfied for \(\alpha>0\) in both weak gravity \(\alpha R\gg1\) and strong gravity \(\alpha R\ll1\) regimes. Performing cosmological perturbation theory as a next step in this analysis can also help probe whether the DOFs \(F\) and \(\phi\) support an early-universe inflationary era by deriving observables and comparing them against exiting datasets [41], [66].

3 Einstein Frame: Metric Formalism↩︎

While Jordan frame can be quite informative in terms of modifications in the gravity sector, cosmological analyses are typically performed by going to the equivalent Einstein frame via Weyl transformations [58]. This is done by first introducing to the action 1 the term \(\frac{F}{2\kappa}(R-\chi^2)\) (where \(\chi\) is an auxiliary field) and then performing a Weyl transformation given by \(g_{\mu\nu}\rightarrow Fg_{\mu\nu}\) to obtain, \[S_E=\int d^{4}x\sqrt{-g}\left[\frac{R}{2\kappa}-\frac{3}{2\kappa F^2}\partial_{\mu}{F}\partial^{\mu}{F}-\frac{1}{2F}\partial{_\mu}\phi\partial{^\mu}\phi-\frac{1}{F^2}\left(\frac{1}{2\kappa}(\chi^2 F-f)+V\right)\right],\] where \(F\equiv F(\chi,\phi)\). The corresponding equations of motion are obtained as, \[\frac{1}{2\kappa}{G_{\mu\nu}}=\frac{3}{\kappa F^2}\partial_{\mu}{F}\partial_{\nu}{F}+\frac{1}{F}\partial{_\mu}\phi\partial{_\nu}\phi-{g_{\mu\nu}}\left[\frac{3}{2\kappa F^2}\partial_{\alpha}{F}\partial^{\alpha}{F}+\frac{1}{2F}\partial_{\alpha}{\phi}\partial^{\alpha}{\phi}+\frac{1}{F^2}\left(\frac{1}{2\kappa}(\chi^2 F-f)+V\right)\right],\] For this tensor equation, the 00 component is given as, \[\frac{-3}{F^2}\dot{F}^2-\frac{\kappa}{F}\dot{\phi}^2-\frac{2\kappa}{F^2}\left(\frac{1}{2\kappa}(\chi^2 F-f)+V\right)=-3{H^2},\] and \(ii\) components are, \[\frac{3}{F^2}\dot{F}^2+\frac{\kappa}{F}\dot{\phi}^2-\frac{2\kappa}{F^2}\left(\frac{1}{2\kappa}(\chi^2 F-f)+V\right)=-2{\dot{H}}-3{H^2}.\] These equations can be used to perform perturbative and late-time analyses similar to the Jordan frame. But instead of these cumbersome equations, we can instead work with their canonicalized fields which greatly simplify calculations and localize modifications compared to standard STTs to certain terms, leaving others intact. Canonicalization is also a necessary step before proceeding with any quantum field theoretical analysis, such as finding unitarity violation scales or checking renormalizability of the theory [70]. Proceeding with canonicalization, we can express, \[\partial_{\mu}{F}=\frac{\partial{F}}{\partial{\phi}}\partial_{\mu}{\phi}+\frac{\partial{F}}{\partial{\chi}}\partial_{\mu}{\chi},\] which gives, \[\begin{align} \partial_{\mu}{\phi}&=\left(\partial_{\mu}{F}-\frac{\partial{F}}{\partial{\chi}}\partial_{\mu}{\chi}\right)\left(\frac{\partial{F}}{\partial{\phi}}\right)^{-1},\\ \frac{3}{\kappa F^2}\partial_{\mu}{F}\partial^{\mu}{F}&=\frac{3}{\kappa F^2}\left[\left(\frac{\partial{F}}{\partial{\phi}}\right)^2{\partial{_\mu}\phi\partial^{\mu}\phi} +\left(\frac{\partial{F}}{\partial{\chi}}\right)^2{\partial{_\mu}\chi\partial^{\mu}\chi} +2\left( {\frac{\partial{F}}{\partial{\phi}}}\cdot{\frac{\partial{F}}{\partial{\chi}}}\right)\partial{_\mu}\phi\partial^{\mu}\chi\right]. \end{align}\] Now, \[\begin{align} \partial_{\mu}{\ln{F}}&=\frac{F_{\phi}}{F}\partial_{\mu}{\phi}+\frac{F_{\chi}}{F}\partial_{\mu}{\chi},\\ \partial_{\mu}{(\ln{F}})\partial^{\mu}(\ln{F})&=\left(\frac{F_{\phi}}{F}\right)^2\partial{_\mu}\phi\partial^{\mu}\phi+\left(\frac{F_{\chi}}{F}\right)^2\partial{_\mu}\chi\partial^{\mu}\chi+2\left(\frac{F_{\chi}F_{\phi}}{F^2}\right)\partial{_\mu}\phi\partial^{\mu}\chi. \end{align}\] Let \(\sqrt{\frac{2}{\kappa}}\ln{F}=\theta\). Then, we have, \[\frac{1}{F^2}\partial_{\mu}{\phi}\partial^{\mu}{\phi}=\exp\left({-\sqrt{\frac{\kappa}{3}}}\theta\right)\partial_{\mu}{\phi}\partial^{\mu}{\phi}=\partial_{\mu}{\phi}\partial^{\mu}{\phi}-\sqrt{\frac{\kappa}{3}}\theta\partial_{\mu}{\phi}\partial^{\mu}{\phi}+\frac{\kappa}{6}\theta^2\partial_{\mu}{\phi}\partial^{\mu}{\phi}+...,\] which gives an infinite series of terms that can be truncated based on the scale of \(F\) and the higher order kinetic interactions can be used to check renormalizability and unitarity violation scales. But this method leaves the fate of the term \((\chi^2F-f)\) uncertain. The canonicalization problem is better resolved assuming \(f(\chi,\phi)\) is factorizable, i.e. \(f(\chi,\phi)=p(\chi)q(\phi)\). In that case, \[\frac{3}{\kappa F^2}\partial_{\mu}{F}\partial^{\mu}{F}=\frac{3}{\kappa} \left[ \left(\frac{q_{\phi}}{q}\right)^2 {\partial{_\mu}\phi\partial^{\mu}\phi} +\left(\frac{p_{\chi^2\chi}}{p_{\chi^2}}\right)^2 {\partial{_\mu}\chi\partial^{\mu}\chi}+2 \left(\frac{p_{\chi^2\chi}}{p_{\chi^2}}\right)\left(\frac{q_{\phi}}{q}\right) {\partial{_\mu}\chi\partial^{\mu}\phi} \right],\] where the subscripts represent variable with respect to which differentiation is being done. Now, considering, \[\sqrt{\frac{3}{\kappa}}\frac{p_{\chi^2\chi}}{p_{\chi^2}}\partial_{\mu}{\chi}=\partial_{\mu}{\Omega},\] or \[\label{eq:p-chi-canonicalization-metric} \partial_{\mu}\left[{\frac{3}{\kappa}}\ln(p_{\chi^2})\right]=\partial_{\mu}{\Omega}\implies p_{\chi^2}=\exp{\left(\sqrt{\frac{\kappa}{3}} \Omega \right)}.\tag{10}\] Similarly, \[\frac{1}{p_{\chi^2} q}\partial_{\mu}{\phi}\partial^{\mu}{\phi}=\exp{\left(-\sqrt{\frac{\kappa}{3}} \Omega \right)}\frac{1}{q}\partial_{\mu}{\phi}\partial^{\mu}{\phi}=\frac{1}{q}\partial_{\mu}{\phi}\partial^{\mu}{\phi}-\sqrt{\frac{\kappa}{3}}\frac{1}{q}\Omega\partial_{\mu}{\phi}\partial^{\mu}{\phi}+\frac{\kappa}{3}\frac{\Omega^2}{q} \partial_{\mu}{\phi}\partial^{\mu}{\phi}.\] Then, the total kinetic term of \(\phi\) is, \[\left(\sqrt{\frac{3}{\kappa}}\frac{q_{\phi}}{q}+\frac{1}{\sqrt{q}}\right)\partial_{\mu}{\phi}=\partial_{\mu}{\theta}.\] Canonicalization is now complete. But the term \((\chi^2F-f)\) still needs to be expressed in terms of canonicalized fields \(\Omega\) and \(\theta\). In factorized form, \[\chi^2F-f=(\chi^2 p_{\chi^2}-p)q.\] Now, if \((\chi^2 p_{\chi^2}-p)\) can be expressed as an explicit function of only \(p_{\chi^2}\), then using canonicalization transformations, \((\chi^2F-f)\) can be recast as an interaction between \(\Omega\) and \(\theta\) fields. While this treatment in Einstein frame is straightforward, the situation is a little more complicated when we follow the Palatini formalism because the auxiliary field is non-dynamical, as we shall see in the next section.

4 Einstein Frame: Palatini Formalism↩︎

Proceeding with the same Weyl transformation of the Jordan frame action 1 but in the Palatini formalism, we arrive at the following Einstein frame action: \[\label{eq:palatini-einstein-action} S=\int d^{4}x\sqrt{-{g}}\left[\frac{{R}}{2\kappa}-\frac{1}{2F}\partial{_\mu}\phi\partial{^\mu}\phi-\frac{1}{F^2}\left(\frac{1}{2\kappa}(\chi^2 F-f)+V\right)\right].\tag{11}\] It is clear that \(\chi\) in this case is a non-dynamical field [41]. Since it is a superfluous DOF in the system, varying the action with respect to \(\chi^2\) can help us obtain a constraint equation for the system, \[\frac{1}{2F^2}\left(\frac{\partial F}{\partial \chi^2}\right) \partial_{\mu}{\phi}\partial^{\mu}{\phi}+\frac{2}{F^3}\left(\frac{\partial F}{\partial \chi^2}\right)\left(\frac{1}{2\kappa}(\chi^2 F-f)+V\right)-\frac{1}{2\kappa F^2}\left(F+\chi^2 \frac{\partial F}{\partial \chi^2} -F\right)=0,\] which on simplification yields, \[\label{eq:palatini-constaint-chi} \kappa F \partial_{\mu}{\phi}\partial^{\mu}\phi + F \chi^2 -2f +4\kappa V =0.\tag{12}\] This equation can be used to eliminate the non-dynamical DOF \(\chi\) from 11 , which is necessary to proceed with any physical analysis in this system. We also require this constraint in order to canonicalize \(\phi\) given in 11 where the field redefinition would obviously depend of \(F(\chi,\phi)\). We will first try to find valid forms of \(p\) for which the the constraint equation can be expressed in a form that could eliminate \(F(\chi,\phi)\). It would be easier to find conditions for which the non-dynamical DOF can be eliminated from the system if we can rewrite the constraint equation 12 by factorizing \(f(\chi,\phi)\) as \(f(\chi,\phi)=p(\chi)q(\phi)\) and rewriting \(p\equiv p(R)\) for typographical ease. Now, we have, \[\begin{align} &\kappa p_{R} q \partial_{\mu}{\phi}\partial^{\mu}{\phi} +\chi^2 p_{R} q -2pq + 4 \kappa V=0,\nonumber\\ \implies &\chi^2 p_{R} - 2p=-\frac{4\kappa V}{q}-\kappa p_{R} \partial_{\mu}{\phi}\partial^{\mu}{\phi}.\label{eq:palatini-constraint-simplified-terms} \end{align}\tag{13}\] Assuming \(p(R)=R\;u(R)\) for the same reasons as earlier, we find, \[p_{R}=u(R)+R\frac{du}{dR}.\] This equation has the same complimentary solution as in metric formalism. However, the terms we need to manage are slightly different. From the left hand side (LHS) of 13 , we have, \[\label{eq:palatini-differential-equation-2} R p_{R}-2p=R^2\frac{du}{dR}-uR.\tag{14}\] Since, we have assumed analyticity of \(p\), we restrict ourselves only to certain forms of \(p\). In theory, many varied forms of this function may exist, each with their signature phenomenological features and high-energy behavior. We also do not consider any forms that may involve weak gravity cut-offs such as [71], [72]. The two simplest cases we analyze are listed below as examples:

• Case 1: When \(p_{R}\) is polynomial in \(R\) as \(p_{R}=a+bR+cR^2+..\), the complete solution has the form, \[u(R)=\frac{m}{R} + Ap_{R} +B.\] Substituting in 14 , we get, \[\begin{align} &Rp_{R} - 2p = R^2 \left[-\frac{m}{R^2} +A(b+2cR+...) \right]-m-AR(a+bR+cR^2+...)-BR,\\ \implies &Rp_{R} - 2p=-2m-(aA+B)R+... \end{align}\] By substituting the RHS from this expression in 13 and solving the resulting polynomial equation in \(\chi^2\), we can completely eliminate \(\chi^2\) using the constraint equation. As an example, consider \(p=R+\alpha R^2\), i.e. a Starobinsky-like gravity with \(a=1\), \(b=2\alpha\), \(A=2\), \(B=-1\), and the other constants set to 0. We, then, find, \[\begin{align} &\chi^2=\frac{4\kappa V+\kappa q\partial_{\mu}{\phi}\partial^{\mu}{\phi}}{q(1-2\kappa \alpha \partial_{\mu}{\phi}\partial^{\mu}{\phi})},\\ \implies &p_{R}q=\frac{q+8\kappa \alpha V }{1-2\kappa \alpha \partial_{\mu}{\phi}\partial^{\mu}{\phi}}. \end{align}\] For higher-order polynomials \(p(R)\), the expressions become extremely cumbersome and may even require additional information about the relative magnitudes of the potentials and kinetic terms to avoid complex solutions of \(\chi\). Fixing these relative magnitudes at this stage would defeat the purpose of the dynamical analysis we shall performed in the next section.

• Case 2: When \(p_R\) is an exponential function of \(R\), \(p_{R}=ae^{\lambda R}\) for which, \[u=\frac{m}{R} + A e^{\lambda R}.\] Following the same procedure as the previous case, we arrive at the following equation: \[-2m+\left[1+ \ln\left(\frac{p_{R}}{a} \right) \right]\frac{A}{a} \ln \left(\frac{p_{R}}{a} \right) p_{R} + \kappa p_{R} \partial_{\mu}{\phi}\partial^{\mu}{\phi} = -\frac{4\kappa V}{q},\] which gives no straightforward solution except the trivial case where \(\lambda=0\), i.e. \(p_{R}=a\) or \(p(R)=aR\). Clearly, eliminating \(\chi\) in exponential models where \(p_{R}=ae^{\lambda R}\) is extremely difficult.

It appears that among the polynomial expansions of \(p(R)\), Starobinsky-like model is a special case for which one can safely remove the non-dynamical DOF \(\chi\) in the Einstein frame without imposing additional constraints on the scalar potentials. We shall, therefore, proceed with this choice for the rest of the analysis. Rewriting the action accordingly,

\[\label{eq:palatini-action} S=\int d^{4}x\sqrt{-{g}} \left[\frac{{R}}{2\kappa} - \frac{1}{2(q+8\kappa \alpha V)} \partial_{\mu}{\phi}\partial^{\mu}{\phi} + \frac{\kappa \alpha }{2(q+8\kappa \alpha V)} (\partial{_\mu}\phi\partial{^\mu}\phi) (\partial{_\nu}\phi\partial{^\nu}\phi) -U \right],\tag{15}\] where, \[U\equiv\frac{V}{q (q+8\kappa \alpha V)},\] is the effective potential. Compared to metric formalism, where a wider variety of \(f(R,\phi)\) models could give favorable conditions, it appears that Palatini formalism demands a smaller subset of analytic functions \(p(R)\) in order for the non-dynamical scalar to be eliminated. Now that we have isolated these conditions solely from consistency requirements, we shall now verify whether 15 in particular can support a proper cosmological evolution from early to late-times.

5 Dynamical Stability Analysis of Einstein Frame Palatini \(f(R,\phi)\) model↩︎

So far, we have left the potential \(V\) and the factorized coupling function \(q\) as arbitrary. We shall continue to do so in this section as well and check if we can find some conditions for both \(V\) and \(q\) that ensure proper late-time acceleration. We will, however, make one well-posed simplification in 15 to make the forthcoming analysis easier. Starobinsky inflation in metric formalism can be obtained from a general \(f(R,\phi)\) theory by assuming that \(q(\phi)=1\) and \(p(R)=R+\alpha_mR^2\), where \(\alpha_m\approx10^9\kappa\) from Planck 2018 [73] constraints on the total energy density of the universe during the inflationary epoch. In terms of \(\alpha_m\), the total energy density during a slow-roll inflation in metric Starobinsky model can be expressed as \(\approx (8\kappa\alpha_m)^{-1}\) [74]. Since observational data is independent of the choice of model or formalism, if we assume that 15 drives early-universe inflation instead, then during inflation, \[\label{eq:inflation-potential-planck} \frac{V}{q (q+8\kappa \alpha V)}\approx \frac{1}{8\kappa\alpha_m},\tag{16}\] where we have assumed that the slow-roll inflationary phase is driven by potential energy density domination. Now, considering the case where \(q\to1\) in Palatini Einstein frame action, we can make the identification \(\alpha\to\alpha_m\) if and only if, \[8\kappa \alpha V\gg q,\] i.e. the scalar potential \(V(\phi)\) (\(\neq0\)) doesn’t contribute and inflation is driven solely by the \(R^2\) term. Thus, Starobinsky inflation in Palatini formalism is expected to show similar behavior to metric formalism even in the absence of a dynamical scalaron mode [75]. Now, if we move away from a pure Starobinsky inflation such that \(\phi\) provides a non-trivial contribution to inflationary dynamics and \(q(\phi)>1\), the dimensionless quantity \(\kappa\alpha V\) must undergo careful calibration to ensure that the condition 16 remains true.

Post-inflation, however, \(U\) must decrease significantly to allow further cosmological evolution. Keeping \(\alpha\) constant during post-inflationary phases, we can claim without loss of generality that, \[\frac{8\kappa\alpha V}{q}\ll1,\] which is especially true during the present epoch where \(U\leq 3H_0^2/\kappa\), and \(H_0\) is the observed Hubble parameter at present time (the inequality suggests that the total energy density in the present epoch may not be dominated by the potential \(U\) alone). We can, then, rewrite the action 15 as, \[\label{eq:palatini-canonicalized-action-simplified} S=\int d^{4}x\sqrt{-{g}} \left[\frac{{R}}{2\kappa} - \frac{1}{2q} \partial_{\mu}{\phi}\partial^{\mu}{\phi} + \frac{\kappa \alpha }{2q} (\partial{_\mu}\phi\partial{^\mu}\phi) (\partial{_\nu}\phi\partial{^\nu}\phi) -U \right]+S_m,\tag{17}\] where \(U= V/q^2\) and \(S_m\) represents non-relativistic dark matter. Note that the only remnant of the Jordan frame \(R^2\) term is now the quartic-order kinetic coupling term. The Friedmann equations for this model containing both scalar field and a perfect fluid in FLRW background 4 are: \[\begin{align} H^2&=\frac{\kappa}{3}\left(\frac{\dot{\phi}^2}{2q}+\frac{3\alpha \kappa \dot{\phi}^4}{2q}+U+\rho_{m} \right),\tag{18}\\ \dot{H}&=-\frac{\kappa}{2} \left(\frac{\dot{\phi}^2}{q}+\frac{2 \alpha \kappa \dot{\phi}^4}{q}+p_m+\rho _m\right). \tag{19} \end{align}\] Here, \(\rho_{m}\) represents the energy density for dark matter (DM) taken as pressureless (\(p_{m} = 0\)) dust with equation of state parameter \(\omega_{m} = 0\). Now, by taking the variation of the action given in 17 with respect to the scalar field \(\phi\), we obtain the Klein-Gordon equation as, \[\ddot{\phi}+3 H \dot{\phi} \frac{(1+2 \alpha \kappa \dot{\phi}^2)}{(1+6 \alpha \kappa \dot{\phi}^2)}+\frac{V_{\phi}}{q(1+6 \alpha \kappa \dot{\phi}^2)}-\frac{q_{\phi} \dot{\phi}^2}{2q} \frac{(1+3 \alpha \kappa \dot{\phi}^2)}{(1+6 \alpha \kappa \dot{\phi}^2)}-\frac{2 q_{\phi} }{q^2} \frac{V}{(1+6 \alpha \kappa \dot{\phi}^2)}=0.\label{eq:phiddot}\tag{20}\] We define the following dimensionless variables, \[\label{eq:dynvar} x\equiv\frac{\kappa\dot{\phi}^2}{6 q H^2},\qquad y\equiv\frac{\alpha \kappa^2 \dot{\phi}^4}{2 q H^2},\qquad z\equiv\frac{\kappa U}{3H^2}, \qquad \lambda\equiv\frac{q_{\phi}}{q} \frac{\dot{\phi}}{H}.\tag{21}\] Using these, from 18 , we obtain the constraint equation: \[x+y+z+\Omega_{m}=1,\] where the relative energy density of matter \(\Omega_{m}=\frac{\kappa \rho_{m} }{3 H^2 }\) is taken to be positive, such that, \[\Omega _{m}=1-x-y-z \leq 1.\] We also find expressions for energy density \(\rho_{\phi}\) and pressure \(p_{\phi}\) respectively as: \[\begin{align} \rho_{\phi}&=\frac{\dot{\phi}^2}{2q}+\frac{3 \alpha \kappa \dot{\phi}^4}{2q}+U,\\ p_{\phi}&=\frac{\dot{\phi}^2}{2q}+\frac{ \alpha \kappa \dot{\phi}^4}{2q}-U. \end{align}\] Now, we can obtain the relevant cosmological parameters using the dynamical variables we defined in 21 . For the scalar field (DE), the density parameter and the equation of state can be respectively expressed as, \[\begin{align} \Omega _{\phi }&=\frac{\kappa \rho _{\phi }}{3 H(t)^2}=(x+y+z), \\ \omega _{\phi }&=\frac{p_{\phi }}{\rho _{\phi }}=\frac{3 x+y-3 z}{3 (x+y+z)}. \end{align}\] In addition, we can also express the effective equation of state parameter \(\omega_{\rm eff}\) as, \[\omega_{\rm eff}=\frac{p_{\phi}+p_{m}}{\rho_{\phi}+\rho_{m}}=x+\frac{y}{3}-z,\] Now, by using the Friedmann equations given in eqs. 18 , 19 and the Klein-Gordon equation 20 , we can obtain the following dynamical system for our analysis, \[\begin{align} x'&=2 I x+J x-\lambda x, \tag{22}\\ y'&=4 I y+J y-\lambda y, \tag{23}\\ z'&=J z-2 \lambda z+z \sigma, \tag{24}\\ \lambda'&=I \lambda +\frac{J \lambda }{2}+\lambda ^2 \rho -\lambda ^2, \tag{25} \end{align}\] where, \[\begin{align} I &=\frac{-3 x-2 y}{x+2 y}+ \frac{\lambda x+\lambda y+2 \lambda z-z \sigma }{2 (x+2 y)},\tag{26}\\ J &= 3 x+y-3 z+3,\tag{27}\\ \sigma& \equiv\frac{V_{\phi}}{V} \frac{\dot{\phi}}{H},\tag{28}\\ \rho &\equiv \frac{q\, q_{\phi\phi}}{q_\phi^2}.\tag{29} \end{align}\] Here \(\sigma\), \(\rho\) are parameters, and the \('\) over the quantity signifies a corresponding derivative with respect to the number of e-folds \(N\), which can be related to cosmic time using the relation \(dN=Hdt\). We will now solve the autonomous system (22 25 ) to find the fixed points and their corresponding eigenvalues, and also identify the conditions under which the model yields late-time accelerated attractors. The critical points along with their respective eigenvalues are given in Tables 1 and 2.

• Critical Point \(P_1\): This point represents a de Sitter vacuum solution dominated by the scalar field (\(\Omega_{\phi}=1\)) with an effective equation of state \(\omega_{\rm eff}=-1\). Its stability is determined by the parameter \(\sigma\); specifically, it possesses a zero eigenvalue, making it non-hyperbolic. However, for \(\sigma<0\), the centre manifold dynamics allow it to function as a stable late-time attractor, physically corresponding to a cosmological constant (this is also demonstrated in the dynamics of the reduced system; see Appendix). Alternatively, if \(\sigma>0\), the point becomes a saddle implying that the current acceleration era may be transient and may lead to a different epoch.

• Critical Point \(P_2\): The point \(P_2\) is a saddle point as it has mixed eigen values and with equation of state parameter \(\omega_{\rm eff}=\frac{1}{3}\), this point cannot serve as a late-time attractor and represents only a transient epoch in the cosmic history.

• Critical Point \(P_3\): This point represents a phantom dark energy attractor (\(\omega_{\rm eff}<-1\)) whenever \(\sigma>0\). In this parameter region, all eigenvalues are negative, implying that \(P_3\) a fully stable and offers a mathematically robust model for phantom acceleration as the final state of the universe. However, for \(\sigma \le 0\), it becomes unstable, thus losing its viability.

• Critical Point \(P_4\): This is the most dynamically complex point in the system, with an effective equation of state and stability dependent on both \(\rho\) and \(\sigma\). It can behave as a stable attractor, a repeller, or a saddle depending on parameters \(\rho\) and \(\sigma\). We illustrate the stability region of the fixed point \(P_4\) in the \((\sigma,\rho)\) parameter space in Fig. 1, where, \[\begin{align} \Delta(\rho,\sigma)=&\;256 \rho^{4} (5\sigma+16)^{2}-1024 \rho^{3}\!\left(16\sigma^{2}+85\sigma+112\right)+32\rho^{2}\!\left(479\sigma^{2}+2056\sigma+2208\right)\nonumber\\ &-64\rho\!\left(96\sigma^{2}+323\sigma+280\right) +873\sigma^{2}+2256\sigma+1600. \end{align}\]

• Critical Point \(P_5\): \(P_5\) can act as saddle or unstable depending on parameter \(\sigma\). Physically, \(P_5\) can behave as an early-time repeller if \(\sigma > -4\) (all positive eigenvalues), and ensures the universe evolves away from this point very quickly. For other choices of \(\sigma\), the point will act as saddle that would represent a transient phase.

• Critical Point \(P_6\): This point corresponds to stiff fluid phase where \(\omega_{\rm eff}=1\). The eigenvalue spectrum includes a positive value, indicating that \(P_6\) is unstable. While such stiff-matter solutions often appear as transient phases in the early universe, this point cannot support late-time acceleration. It functions dynamically as a saddle, bridging earlier cosmological epochs but never serving as a final attractor.

• Critical Point \(P_7\): It offers a phenomenologically promising alternative to the phantom point \(P_3\). It possesses the same effective equation of state form but distinct stability properties. For \(-3<\sigma<0\), it acts as a stable, accelerating attractor with quintessence-like behavior (\(-1<\omega_{\rm eff}<-1/3\)). This makes \(P_7\) another strong candidate for describing observed cosmic acceleration.

Figures 2 and 3 can be shown to represent two distinct scenarios corresponding to points \(P_1\) and \(P_3\) respectively depending on the parameters \(\sigma\) and \(\rho\). Please note that even though some points appear non-hyperbolic (with one or more zero eigenvalues), their stability has been verified in the reduced parameter space (see Appendix). Also note that even though \(\rho\) and \(\sigma\) may appear to be dynamical, their treatment here is as constant parameters. This is because if they are considered dynamical, then \(\rho'\) and \(\sigma'\) would contain higher derivatives of potentials \(V(\phi)\) and \(q(\phi)\) with respect to \(\phi\) which would require us to define variables for successive derivatives of these potentials to study their evolution. However, \(\rho\) and \(\sigma\) being positive or negative near a fixed point could have profound implications on the behavior of \(V(\phi)\) and \(q(\phi)\).

The analysis performed so far clearly establishes the theoretical consistency of the model and demonstrates the existence of a stable late-time accelerating solution. However, this analysis alone is qualitative and does not determine whether the predicted expansion history agrees with observations. To test the observational viability of the model and constrain its parameters, we now confront it with late-time cosmological data in the next section.

6 Data Analysis↩︎

Our subsequent objective requires performing data analysis to constrain the parameter space we have in the model. For achieving this, deriving an analytical expression for the Hubble parameter as a function of redshift is necessary. While this process is quite straightforward for the \(\Lambda\)CDM model, it poses challenges in our \(k\)-essence model owing to the inability to obtain an analytical solution to the Friedmann equation. Therefore, we opt for a dynamical systems approach employing the dynamical variables obtained in section 5. Solving these coupled differential equations, we capture the evolution of the Hubble parameter. In our analysis, we use the following publicly available late time datasets,

1. DESI:- We use the 13 DESI-BAO DR2 ([76], [77]) measurements across the redshift range \(0.1 < \boldsymbol{z} < 4.2\) obtained from observations of about 14 million galaxies and quasars which include bright galaxy sample (BGS), luminous red galaxies (LRG), emission line galaxies (ELG), quasars (QSO), and Lyman-\(\alpha\) tracers. These measurements are given in terms of the volume averaged distance \(D_{\rm V}(\boldsymbol{z}) / r_{\rm d}\), angular diameter distance \(D_{\rm M}(\boldsymbol{z}) / r_{\rm d}\) and comoving Hubble distance \(D_{\rm H}(\boldsymbol{z}) / r_{\rm d}\), where \(r_{\rm d}\) is the sound horizon at the drag era.

2. Cosmic Chronometers:- The Cosmic chronometers (CC) approach allows us to obtain observational values of the Hubble function at different redshifts \(\boldsymbol{z}\leq2\) directly. Since these measurements are independent of any cosmological model and Cepheid distance scale, they can be used to place better constraints on it. In the present analysis, we measure \(H(\boldsymbol{z})\) using the CC covariance matrix [78]–[80]

3. SNeIa:- We use the PantheonPlus (PP) dataset [81], which contains 1701 light curves for 1550 spectroscopically confirmed Type Ia supernovae (SNeIa) covering the redshift range \(0.001 < \boldsymbol{z} < 2.26\). Additionally, we incorporate the Union3 compilation comprising 2087 SNe [82]. We also use the DESY5 sample comprising 1635 photometrically classified SNe from the released part of the full 5 year data of the Dark Energy Survey collaboration (with redshifts in the range 0.1 \(< \;\boldsymbol{z}\;<\) 1.3), complemented by 194 low-redshift SNe from the CfA3 [83], CfA4 [84], CSP [85], and Foundation [86] samples (with redshifts in the range 0.025 \(< \;\boldsymbol{z}\;<\) 0.1), for a total of 1829 SNe [87].

To constrain the free and derived parameters in our model, we perform Markov-Chain-Monte Carlo (MCMC) simulations using the publicly available tool COBAYA [88], [89]. We will be analysing the general structure of the dynamical equations obtained in eqs. (22 25 ). Further, we rewrite them, by replacing the dynamical variables \(x\) and \(y\) with cosmologically relevant quantities such as \(\Omega_\phi\) and \(\omega_\phi\). Employing the numerical integration routines in Scipy [90], we solve the dynamical system starting from \(\ln a = -4\), and evolve it to the present time. Thus our parameter space consists of the initial values of the four dynamical variables \(\Omega_\phi^i\), \(\omega_\phi^i\), \(z_i\) and \(\lambda_i\), and the parameters \(\rho\) and \(\sigma\) that appear in their dynamical equations. Further, we are taking \(r_d\), the comoving sound horizon at the drag epoch, as a free parameter for analysing the BAO data from DESI. We constrain the parameters with uniform priors: \(\Omega_\phi^i \in [0, 0.3]\), \(\omega_\phi^i \in [-1,1]\), \(z_i \in [10^{-6}, 10^{-3}]\), \(\lambda_i \in [10^{-8},1]\), \(\rho \in [-10,2]\), \(\sigma \in [-10,2]\) and \(r_d \in [100, 200]\). The convergence of chains is ensured by having the Gelman-Rubin criterion \(|R-1| \leq 0.01\). We utilize GetDist [91] and BOBYQA [92], [93] to analyze and minimize the chains. The marginalised 1-D and 2-D posterior distributions of the dynamical variables and other parameters in our model for some of the data set combinations are shown in Fig. (5). In Table 3, we report the marginalized parameter constraints along with their \(1-\boldsymbol{\sigma}\) errors.

We can see from Table 3 that the combination of DESI with CC can put only a bound on the model parameter \(\rho\). However, with the inclusion of a supernova dataset to this combination, \(\rho\) can be better constrained in our framework. The same can be identified from the corner plot given in Fig. (5). Constraints can be obtained on the other model parameters as well. However, the parameters \(z_i\) and \(\lambda_i\) do not show any constraints and so we haven’t reported them in our results. The evolution of the cosmological parameters in our \(k\)-essence model, such as the fractional energy density and the effective equation of state, in accordance with the constraints obtained from various dataset combinations, is shown in Fig. 4. In this plot, the evolution is depicted against the number of e-folds \(N=\ln a\). The results demonstrate that the model exhibits a matter-dominated phase (i.e., \(\Omega_m>\Omega_\phi\)), for all the datasets, before transitioning to an accelerating epoch. During this matter dominated phase, the effective equation of state is close to zero, indicating attractive scaling behavior. As the system enters the late-time phase, the model shows accelerated expansion with \(\omega_{\rm eff}<-1/3\). Further, \(\omega_{\rm eff}\) is seen to approach to a value of -1 asymptotically in the future.

7 Discussion↩︎

\(f(R,\phi)\) gravity offers rich phenomenology and can even explain early- and late-time cosmic evolution under one unified umbrella. In most literature, \(f(R,\phi)\) theories are studied by assuming specific form factors. In this paper, however, we attempted to constrain Palatini \(f(R,\phi)\) theories in the Einstein frame based on consistency requirements, dynamical stability at late-time, as well as observational data to find compatibility with various parameters. We have used late time data sets to capture the behaviour of the cosmological and other parameters we have in our model. We have used a general potential in our analysis where it appear as a dynamical variable while its derivatives remain as parameters. With the observational data, we could not obtain robust constraints on all the parameters. However, we were able to constrain some of them, while obtaining bounds for the others. Particularly, the estimation of the parameters \(\rho\) and \(\sigma\) are important since they determine the stability of the fixed points in the dynamical system, thereby controlling the evolution.

We found that among the fixed points listed in Table 1, only \(P_1\) and \(P_3\) provide a stable late-time acceleration era with behavior summarized in Table 2 and plotted explicitly in Figs 2 and 3 respectively. Now, as mentioned earlier, \(\rho\) and \(\sigma\) are treated as constants throughout the fixed point analysis and their values are later fixed at \(N=-4\) using observational data in Table 3. Within 68% CL, we see that \(\rho<0\). Based on its definition and knowing that \(q\to1\) is the minimum allowed value to maintain consistency with general relativity, we find that \(q_{\phi\phi}<0\). Also, even though observations seem to favor a negative \(\sigma\) (i.e. \(V_\phi<0\)), small positive values may still be allowed.

If we consider \(\sigma<0\), we rule out \(P_3\) as a suitable fixed point since it demands \(\sigma>0\). Below we list the various parameter values for \(P_1\) as a stable fixed point and their corresponding implications: \[\begin{align} x=-2,\;y=3& \qquad\implies\qquad \dot{\phi}^2\to-(24\alpha \kappa)^{-1},\tag{30}\\ z=0& \qquad \implies\qquad U= Vq^{-2}\to0,\tag{31}\\ \lambda=0&\qquad\implies\qquad q_\phi\to0,\tag{32}\\ \sigma<0&\qquad\implies\qquad V_\phi<0,\tag{33}\\ \rho<0&\qquad\implies\qquad q_{\phi\phi}<0.\tag{34} \end{align}\] Moving step by step, conditions 32 and 34 imply at the coupling potential \(q(\phi)\) approaches its maximum value at \(P_1\). The scalar potential \(V(\phi)\) is found to decrease in magnitude as \(\phi\) increases, as shown in 33 and it approaches zero as per 31 . Also, the combination \(x<0\) and \(y>0\), given their definitions in 21 implies that \(\dot{\phi}^2<0\) which is characteristic of phantom fields (also implied by the effective equation of state parameter in Fig. 2). One could argue that we could consider that \(\dot{\phi}^2\) is positive and \(\alpha<0\) instead to eliminate the phantom menace. This, however, would mean that for \(x=-2\) and \(y=3\), \(q(\phi)\) would have to be negative. In Einstein frame, this would imply that the signs of the kinetic term and quartic kinetic coupling in the scalar part of action 15 are flipped, leading to \(\phi\) becoming a ghostly DOF⁵. Since, instabilities are, therefore, present in the system nonetheless, \(\dot{\phi}^2<0\) is a safer choice.

If, on the other hand, we consider \(\sigma\gtrsim0\) (which is allowed within 68% CL), the viable late-time attractor would be \(P_3\). Below we list its corresponding parameter values along with their implications: \[\begin{align} x=0,\;y=0& \qquad\implies\qquad \dot{\phi}^2\to0\qquad\implies\qquad\dot{\phi}\to0^{\pm},\tag{35}\\ z=1& \qquad \implies\qquad U= Vq^{-2}\to 3M_{\rm Pl}^2H^2,\tag{36}\\ \lambda=0&\qquad\implies\qquad q_\phi\to0,\tag{37}\\ \sigma\gtrsim0&\qquad\implies\qquad V_\phi>0\;\text{and}\;\dot{\phi}\gtrsim0,\quad \text{or}\quad V_\phi<0\;\text{and}\;\dot{\phi}\lesssim0,\tag{38}\\ \rho<0&\qquad\implies\qquad q_{\phi\phi}<0.\tag{39} \end{align}\] Here, conditions 37 and 39 again imply that \(q\) is approaching its maximum value, as for \(\sigma<0\). From condition 38 , we find two completely different cases: either \(\phi\) slowly increases in magnitude with time and \(V\) increases in magnitude with \(\phi\) or vice versa. In both cases, we can at least guarantee that \(V\) is directly proportional to \(\phi\). Since condition 36 implies that the scalar potential \(V(\phi)\) is approaching its maximum value \(3q^2M_{\rm Pl}^2H^2\) at late times, the only viable combination is \(V_\phi>0\) and \(\dot{\phi}\gtrsim0\).

An important thing to note is that the scale of \(N\) in Fig. 3 is much larger than the scale of \(N\) in Fig. 2. From Table 2, one can see that \(P_1\) could behave as a saddle point for \(\sigma>0\). As such, Figs 2 and 3 could imply that \(P_1\) is simply a transient stage before the universe settles into \(P_3\) at late-times. We reiterate that throughout the analysis we have treated \(\sigma\) as a constant parameter due to restrictions in finite-dimensional dynamical systems analysis. If \(\sigma\) is treated as a dynamical variable, it may eventually change behavior from \(\sigma<0\) (as seen from data analysis performed in Section 6 for \(N=-4\)) to \(\sigma>0\) much later than the present epoch. We cannot rule out this conclusion with absolute certainty since it requires us to define an infinite-dimensional dynamical system for general potentials \(V\) and \(q\). This is beyond the scope of the present work, but it could be pursued in the future as a follow-up.

Alternatively, if one assumes specific forms of the potentials \(V(\phi)\) and \(q(\phi)\), such that they conform to the behavior exhibited by points \(P_1\) and \(P_3\), the number of variables could become finite and the system could be solved using the methods presented in this work. But since we sought to study the behavior of the general \(f(R,\phi)\) theories with arbitrary potentials in this paper, we shall leave that analysis as a future work as well.

Appendix↩︎

Here, we present dynamical evolution of the system described in Section 5 under various constraints. This analysis is intended to demonstrate that the points that appear non-hyperbolic in Table 2 can become hyperbolic in the constrained parameter space.

7.1 \(V(\phi)=0\) and \(\lambda\) as Constant↩︎

The autonomous system (22 26 ) admits three critical points in the \((x,y)\) phase space, which are given in Table-4 and corresponding eigenvalues and stability conditions are given in Table-5. In the following, we analyse the stability and physical implications of each point separately.

• Critical Point \(A_{1}\): The fixed point \(A_{1}\) corresponds to a scalar field dominated state with \(\Omega_{\phi}=1\) with effective equation of state parameter \(\omega_{\phi}=\omega_{\rm eff}=1/3\). Its eigenvalues \((1, 2-\lambda/2)\) show that \(A_{1}\) is an unstable node for \(\lambda<4\) and a saddle for \(\lambda>4\). Thus, it cannot act as a late time attractor, but it can act as early time repeller for \(\lambda<4\) or a transient early time phase for \(\lambda>4\).

• Critical Point \(A_{2}\): Point \(A_{2}\) also satisfies \(\Omega_{\phi}=1\) and features an effective equation of state \(\omega_{\rm eff}=\lambda/3-1\), covering de Sitter (\(\lambda=0\)), quintessence-like (\(\lambda<2\)), and phantom-like (\(\lambda<0\)) regimes. The eigenvalues \((\lambda-3)\) and \((\lambda-4)(\lambda-6)/(\lambda-8)\) imply stability for \(\lambda<3\), and for \(\lambda<2\) the point becomes a stable accelerating attractor. These properties make \(A_{2}\) a strong candidate for the universe’s late-time evolution.

• Critical Point \(A_{3}\): The point \(A_{3}\) yields \(\Omega_{\phi}=1\) and a stiff-fluid equation of state \(\omega_{\rm eff}=1\), indicating strong deceleration. Its eigenvalues \((3, \lambda-6)\) ensure at least one positive eigenvalue for all \(\lambda\), rendering it a saddle for \(\lambda<6\) and an unstable node for \(\lambda>6\). Consequently, \(A_{3}\) cannot serve as a late-time attractor and is relevant only as a possible early stage phase.

Fig.6 illustrates the numerical evolution of the dynamical variables \(\Omega_\phi\), \(\Omega_m\), \(\omega_{\rm eff}\), \(x\), and \(y\) with respect to the number of e-folds \(N\). The plots correspond to the initial conditions \(\Omega_\phi(0)=0.68\), \(\omega_\phi(0)=-0.99\), and \(\lambda = 0.01\). As shown in the left panel, the matter density parameter \(\Omega_m\) decreases while \(\Omega_\phi\) increases and asymptotically approaches unity, causing the effective equation of state \(\omega_{\rm eff}\) to evolve towards \(-1\), signalling the onset of a late-time accelerating phase. The right panel shows the evolution of the variables \(x\) and \(y\), which settle to constant values associated with the stable attractor \(A_2\). Overall, the figure demonstrates the transition from a matter-dominated epoch to a scalar-field dominated accelerating regime for the chosen initial conditions.

7.2 \(V(\phi)=0\) and \(\lambda\) as variable↩︎

In this subsection we analyze the stability and cosmological interpretation of the fixed points of the three–dimensional autonomous system in the variables \((x,y,\lambda)\) for the case \(V=0\).

• Critical Point \(B_1\): This is corresponding to a de Sitter phase as \(\omega_{eff}=-1\). The eigenvalues \(\{-3,-3,0\}\) indicate two stable directions and one marginal direction. Hence \(B_1\) behaves as a late time attractor in the \((x,y)\) plane and is the only fixed point providing accelerated expansion. It is, therefore, the unique viable dark energy solution.

• Critical Point \(B_2\): The point \(B_2\) has all eigen values to be positive with effective equation of state parameter \(\omega_{eff}=\frac{1}{3}\), which makes it highly unstable point, it can act as early time repeller.

• Critical Point \(B_3\): This points has same behaviour as \(P_2\). The eigenvalues \(\left\{\frac{4-8\rho}{3-4\rho}, -1, 1 \right\}\) always include both positive and negative modes for all \(\rho \neq \tfrac{3}{4}\); hence \(B_3\) is a saddle point. It can act as a transient phase.

• Critical Point \(B_4\): This is corresponding to a stiff fluid \(\omega_{eff}=1\) with mixed eigenvalues \(\{-6,3,0\}\) classifying \(B_4\) as a saddle point. It may appear only as a transient early epoch and does not describe the present Universe.

Figure 7 shows how the quantities \(\Omega_\phi\), \(\Omega_m\), \(\omega_{\rm eff}\), \(x\), \(y\) and \(\lambda\) change with the number of e-folds \(N\). The plots are generated using the initial conditions \(\Omega_\phi(0)=0.68\), \(\omega_\phi(0)=-0.99\), and \(\rho = 0.1\). From the left panel, we see that the matter density \(\Omega_m\) gradually decreases, while the scalar-field density \(\Omega_\phi\) increases and approaches unity at late times. As this happens, the effective equation of state \(\omega_{\rm eff}\) moves towards \(-1\), indicating that the system enters an accelerated expansion phase. The right panel shows the evolution of the variables \(x\) and \(y\). Both of them smoothly settle to constant values, which correspond to the late-time attractor \(B_1\). Overall, the figure shows that, for these initial conditions, the system naturally evolves from a matter-dominated stage to a scalar-field dominated accelerating phase.

7.3 \(V(\phi) \neq 0\) and \(\lambda\) as Constant↩︎

In this case we are considering x, y and z to be dynamical variables while keeping \(\lambda\) to be a constant parameter. The critical points extracted from the set of equations (22 24 ) along with their respective eigenvalues are given in Table 8 and Table 9.

• Critical Point \(C_1\): The eigenvalues of \(C_1\) are \(\left\{ 1,\; \frac{4-\lambda}{2},\; 4-2\lambda+\sigma \right\}\), where the positive eigenvalue \(1\) immediately indicates that the point is unstable/saddle. Additionally, at \(C_1\), the effective equation of state \(w_{\text{eff}} = \frac{1}{3}\) corresponds to a decelerating phase. Therefore, this point is not suitable for modeling late-time acceleration. Instead, it can act as an early time repeller when \(\lambda < 4\) and \(\sigma > 2\lambda - 4\), while for other parameter choices it corresponds to a transient saddle phase.

• Critical Point \(C_2\): The critical point \(C_2\) has eigenvalues \(\left\{ \frac{\lambda-\sigma}{2},\, -4+2\lambda-\sigma,\, -3+2\lambda-\sigma \right\}\), and stability requires all of them to be negative, which leads to the condition \(\sigma>\max\{\lambda,\,2\lambda-3\}\). The effective equation of state is \(w_{\text{eff}}=\frac{1}{3}(2\lambda-\sigma-3)\), and accelerated expansion requires \(\sigma>2\lambda-2\). Therefore, critical point \(C_2\) corresponds to a stable late-time accelerating solution for \(\sigma>\max\{\lambda,\,2\lambda-2\}\).

• Critical Point \(C_3\): The critical point \(C_3\) has eigenvalues \(\{\lambda-3,\; \frac{(\lambda-4)(\lambda-6)}{\lambda-8},\; -\lambda+\sigma\}\) and an effective equation of state \(w_{\mathrm{eff}}=\frac{\lambda}{3}-1\). \(C_3\) behaves as a stable late-time attractor if \(\lambda<3\) and \(\sigma>\lambda\) satisfied. In addition, the solution corresponds to an accelerated expansion phase when \(\lambda<2\). Therefore, the critical point \(C_3\) represents a late time stable accelerating solution for \(\lambda<2\) and \(\sigma>\lambda\).

• Critical Point \(C_4\): The eigenvalues of \(C_4\) are \(\{3, \, -6 + \lambda, \, 6 - 2\lambda + \sigma \}\), where the positive eigenvalue \(3\) immediately indicates that the point is unstable/saddle. Additionally, at \(C_4\), the effective equation of state \(w_{\text{eff}} = 1\) corresponds to a decelerating phase, incompatible with the late-time accelerated expansion seen in the universe. Therefore, this point is not suitable for modeling late-time acceleration. Instead, it can act as an early time repeller when \(\lambda > 6\) and \(\sigma < 2\lambda -6\), while for other parameter choices it corresponds to a transient saddle phase.

• Critical Point \(C_5\): The critical point \(C_5\) is defined by the eigenvalues \((-\lambda+\sigma)\), \((-6+2\lambda-\sigma)\), and \((-3+2\lambda-\sigma)\), with an effective equation of state given by \(w_{\mathrm{eff}}=-1+\frac{2\lambda-\sigma}{3}\). This point is stable when all eigenvalues are negative, which occurs for \(2\lambda-3<\sigma<\lambda\), implying that the system naturally evolves toward \(C_5\) at late times. Moreover, the solution corresponds to an accelerating universe when \(w_{\mathrm{eff}}<-1/3\), or equivalently when \(\sigma>2\lambda-2\). As a result, for parameter values satisfying \(2\lambda-2<\sigma<\lambda\), the critical point describes the late time accelerated phase.

Figure 8 shows how the system evolves for the chosen parameters \(\lambda = 0.01\) and \(\sigma = -3\). The plots clearly indicate that the universe moves from an initial transient stage towards a stable final state. The scalar-field density \(\Omega_{\phi}\) slowly increases and becomes dominant, while the matter density \(\Omega_{m}\) decreases, and the effective equation of state approaches a value very close to \(-1\), showing that accelerated expansion is reached. The variables \(x\) and \(y\) also settle to constant values after some time, which means that the system is approaching a fixed point. The variable \(z\) quickly decreases to zero, matching the prediction from the stability analysis. Taken together, all three figures confirm that the system evolves towards the only stable critical point for these parameter values, \(C_{3}\), which represents a late-time accelerating solution.

Figure 9 illustrates the dynamical evolution of the system for the initial conditions \(\Omega_{\phi}(0)=0.68\), \(\omega_{\phi}(0)=-0.99\) and \(z(0)=10^{-4}\), with the parameter values \(\sigma=0.02\) and \(\lambda=0.01\). Figure 4(a) shows that the scalar-field energy density \(\Omega_{\phi}\) approaches unity, while the matter density \(\Omega_{m}\) decays to zero, causing the effective equation of state to settle at \(\omega_{\rm eff}=-1\), which is indicative of a accelerated expansion. Figure 4(b) displays the evolution of the dynamical variables \(x\) and \(y\), both of which approach towards zero as the system goes in far future. Figure 4(c) shows that the variable \(z\) grows monotonically from its initially negligible value and asymptotically approaches \(z\simeq 1\). Altogether, these behaviors confirm that the system evolves toward the stable late-time attractor \(C_{2}\), corresponding to a scalar-field dominated Universe with effective cosmological-constant behavior.

References↩︎

---

1. rahul19@iiserb.ac.in↩︎

2. abhijith.ajith.421997@gmail.com↩︎

3. sukanta@iiserb.ac.in↩︎

4. architmedes@gmail.com↩︎

5. A negative \(q(\phi)\) in the Jordan frame action would imply that gravity becomes repelling at late-times.↩︎