Preventing Curvature Singularities in $f(R)$ Dark Energy Models

TL;DR

This paper explores how adding higher-order curvature corrections to $f(R)$ gravity models can prevent curvature singularities, making these models more viable for describing dark energy across different cosmological epochs.

Contribution

It introduces specific higher-order corrections compatible with inflation to stabilize $f(R)$ dark energy models and prevent singularities.

Findings

Higher-order corrections can effectively prevent curvature singularities.

Models remain consistent from inflationary to late-time cosmology.

Corrections improve the viability of $f(R)$ dark energy models.

Abstract

The curvature singularity problem in $f(R)$ dark energy models poses a significant challenge to their viability as alternatives to the $\Lambda$CDM paradigm. In this work, we investigate the possibility of resolving this issue by incorporating higher-order corrections that are compatible with the inflationary phase. We analyze the effects of adding $R^2$, $R^{\frac{m+2}{m+1}}$ and $\alpha$-attractor representation in $f(R)$ gravity terms to types of dark energy $f(R)$ models, focusing on their ability to prevent singularities in high-density environments. Our results demonstrate that these corrections can effectively stabilize the models, ensuring their consistency across both inflationary and late-time cosmological scales.

Equations

$$S=\frac{1}{2{\kappa }^2}\int d^{4}x\sqrt{-g}f(R)+S_M,$$

$$f_R{R}_{\mu \nu }-\frac{f}{2}{g}_{\mu \nu }-\left({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square \right)f_R={\kappa }^2{T}_{\mu \nu }^{(M)},$$

$${\tilde{g}}_{\mu \nu }={\Omega }^2{g}_{\mu \nu },{\Omega }^2=f_R(R).$$

$$S_E=\int d^{4}x\sqrt{-\tilde{g}}\left(\frac{1}{2{\kappa }^2}\tilde{R}-\frac{1}{2}{\tilde{g}}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -V(\varphi )\right),$$

$$\varphi =\sqrt{\frac{3}{2}}\frac{1}{\kappa }\ln f_R(R).$$

$$V(\varphi )=\frac{R{f}_R(R)-f(R)}{2{\kappa }^2(f_R(R){)}^2}.$$

$$ϵ=\frac{1}{2{\kappa }^2}{\left(\frac{V^{\prime }(\varphi )}{V(\varphi )}\right)}^2,\eta =\frac{1}{{\kappa }^2}\frac{V^{\prime \prime }(\varphi )}{V(\varphi )},$$

$$N=\ln \left(\frac{a_f}{a_i}\right)\simeq {\kappa }^2{\int }_{{\varphi }_{\mathrm{end}}}^{{\varphi }^{\ast }}\frac{V(\varphi )}{V^{\prime }(\varphi )}d\varphi ,$$

$$n_s=1-6ϵ+2\eta ,r=16ϵ,A_s=\frac{{\kappa }^{4}V(\varphi )}{24{\pi }^2ϵ}.$$

$$f^p(R)$$

$$f^e(R)$$

$$f^p(R)$$

$$f^e(R)$$

$$R{f}_R-2f-R+3\square f_R={\kappa }^{2}T,$$

$$T(t)=-T_0\left(1+\frac{t}{t_{ch}}\right),$$

$$\square f_R=-{\partial }_t^2(f_R)=\frac{1}{3}[{\kappa }^{2}T-R{f}_R+2f+R].$$

$$y_p^{\prime \prime }+2n\frac{{y_p^{\prime }}^2}{y_p}+\frac{({\kappa }^2{T}_0){\gamma }_p^2{\beta }^{2n+1}}{6n\lambda (2n+1)}{y}_p^{-2n}\left[\left(1+\frac{\tau }{{\tau }_{\text{ch}}}\right)-\frac{1}{y_p}\right]=0,$$

$$\beta =\frac{{\kappa }^2{T}_0}{R_{\text{ch}}},{\gamma }_p^2=\frac{6n(2n+1)\lambda }{R_{\text{ch}}}{\beta }^{-(2n+2)}.$$

$$y_p^{\prime \prime }+2n\frac{{y_p^{\prime }}^2}{y_p}+y_p^{-2n}\left[\left(1+\frac{\tau }{{\tau }_{\text{ch}}}\right)-y_p^{-1}\right]=0,$$

$$y_e^{\prime \prime }+(\beta y_e^{-1}-2)\frac{{y_e^{\prime }}^2}{y_e}+y_e^2{e}^{\beta y_e^{-1}}\left[\left(1+\frac{\tau }{{\tau }_{\text{ch}}}\right)-y_e^{-1}\right]=0,$$

$${\gamma }_e^2=\frac{3\lambda }{R_{\text{ch}}}.$$

$$f(R)=R+\frac{R^2}{6M^2},$$

$$n_s\approx 1-\frac{2}{N},r\approx \frac{12}{N^2}.$$

$$n_s\approx 0.960-0.967,r\approx 0.003-0.004,$$

$$f(R)=R-\lambda R_{\text{ch}}\left[1-{\left(\frac{R_{\text{ch}}}{R}\right)}^{2n}\right]+\frac{R^2}{6m_R^2},$$

$$y_p^{\prime \prime }+2n\frac{{y_p^{\prime }}^2}{y_p}+g_{pI}{y}_p^{-(2n+2)}\left[y_p^{\prime \prime }-\frac{2{y_p^{\prime }}^2}{y_p}\right]+y_p^{-2n}\left[\left(1+\frac{\tau }{{\tau }_{\mathrm{ch}}}\right)-y_p^{-1}\right]=0,$$

$$g_{pI}=\frac{R_{\mathrm{ch}}}{m_R^{2}n\lambda (2n+1)}{\beta }^{2n+2}.$$

$$f(R)=R-\lambda R_{\mathrm{ch}}\left(1-e^{-R/R_{\mathrm{ch}}}\right)+\frac{R^2}{6m_R^2},$$

$$y_e^{\prime \prime }+(\beta y_e^{-1}-2)\frac{{y_e^{\prime }}^2}{y_e}+g_{eI}{e}^{\beta \left(\frac{1}{y_e}-1\right)}\left(y_e^{\prime \prime }-\frac{2{y_e^{\prime }}^2}{y_e}\right)+y_e^2{e}^{\frac{\beta }{y_e}}\left[\left(1+\frac{\tau }{{\tau }_{\mathrm{ch}}}\right)-y_e^{-1}\right]=0,$$

$$g_{eI}:=\frac{2R_{\mathrm{ch}}}{3m_R^2\lambda }{e}^{\beta }.$$

$$f(R\gg \gamma )\simeq \gamma \left(\frac{m+1}{m+2}\right){\left(\frac{1}{4(m+2)}\right)}^{\frac{1}{1+m}}{\left(\frac{R}{\gamma }\right)}^{\frac{m+2}{m+1}}.$$

Full text

11footnotetext: * Preprint of an article submitted for consideration in International Journal of Modern Physics A © 2025 World Scientific Publishing Company. https://www.worldscientific.com/worldscinet/ijmpa

1 Introduction

Recent cosmological observations, such as those from the Planck satellite and Type Ia supernovae, strongly support the conclusion that the late-time universe is undergoing accelerated expansion [1]. This phenomenon is commonly attributed to the presence of dark energy, a mysterious component that constitutes approximately 70% of the universe’s energy density. The simplest explanation for dark energy is the cosmological constant $\Lambda$, which fits well within the $\Lambda$CDM model. However, the cosmological constant faces theoretical challenges, such as the fine-tuning problem and the coincidence problem [2, 3]. Alternative explanations include dynamical dark energy models, modified gravity theories, and interactions between dark energy and dark matter [4, 5].

One promising alternative to the cosmological constant is $f(R)$ gravity, a modification of General Relativity (GR) where the gravitational action is generalized to an arbitrary function of the Ricci scalar $R$ [6, 7]. This theory naturally accommodates dark energy without the need for additional fields, as the modified dynamics of gravity can drive the late-time accelerated expansion. Several successful $f(R)$ models, such as the Hu-Sawicki model [8] and the Starobinsky model [9], have been proposed and tested against observational data. However, these models often suffer from curvature singularities, where the Ricci scalar $R$ diverges in high-density environments, leading to unphysical predictions [10, 11]. Previous attempts to resolve these singularities include introducing additional terms in the $f(R)$ function or modifying the theory to ensure stability under cosmological perturbations [12, 13, 10, 14, 15].

In this paper, we address the curvature singularity problem in $f(R)$ dark energy models by incorporating higher-order corrections that are compatible with the inflationary phase. Our goal is to demonstrate that such corrections can effectively cure the singularities while maintaining the viability of the models across different cosmological scales [16, 15, 17, 3, 18, 19]. We explore several modifications, including the addition of $R^{2}$, $R^{\frac{m+2}{m+1}}$ terms, and the $\alpha$-attractor polynomial extension, and analyze their impact on the stability and dynamics of the models [20]. By doing so, we aim to identify a viable $f(R)$ gravity model that can consistently describe both the early inflationary epoch and the late-time accelerated expansion, while satisfying all cosmological constraints. The paper is organized as follows: in Section 2, we review the theoretical framework of $f(R)$ gravity and its application to dark energy. In Section 3, we discuss the curvature singularity problem in detail. Sections 4, 5 and 6 present our proposed solutions, and Section 7 concludes with a discussion of the implications of our findings.

2 Viable $f(R)$ Gravity Theories

2.1 Overview of $f(R)$ Gravity Models

In modern cosmology, $f(R)$ gravity theories have gained prominence as a compelling extension of General Relativity (GR), achieved by generalizing the Einstein-Hilbert action. Rather than adopting the Ricci scalar $R$ as in GR, these theories introduce a general function $f(R)$ of the scalar curvature. The action for $f(R)$ gravity is written as

[TABLE]

where $\kappa^{2}=8\pi G$, $G$ denotes Newton’s gravitational constant, $g$ is the determinant of the metric tensor $g_{\mu\nu}$, and $S_{M}$ represents the matter action, which includes both relativistic and non-relativistic components. Varying this action with respect to the metric $g_{\mu\nu}$ yields the modified field equations,

[TABLE]

where $f_{R}\equiv df(R)/dR$, $\nabla_{\mu}$ denotes the covariant derivative, $\Box\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the d’Alembertian operator, and $T^{(M)}_{\mu\nu}$ is the energy-momentum tensor for matter. Unlike the second-order Einstein field equations in GR, these equations are fourth-order in derivatives of the metric, leading to a richer but more complex dynamical structure [7].

2.2 Inflation in $f(R)$ Gravity

Inflation refers to a brief epoch of accelerated cosmic expansion in the early universe and is now a cornerstone of the standard cosmological paradigm. It resolves several major issues in the original Big Bang framework, including the horizon, flatness, and monopole problems [21, 22]. While the conventional picture relies on a scalar field (the inflaton) within GR to drive this expansion, $f(R)$ gravity offers a purely geometric alternative for realizing inflationary dynamics [23, 24]. In this framework, the inflationary behavior emerges from the additional scalar degree of freedom inherently present in the gravitational sector via the $f(R)$ function [25, 26].

To facilitate the analysis, it is convenient to work in the Einstein frame, obtained through a conformal transformation of the metric:

[TABLE]

Under this transformation, the action becomes

[TABLE]

where $\phi$ is the scalar field (often called the scalaron) defined by

[TABLE]

The corresponding scalaron potential is given by

[TABLE]

This formulation establishes the equivalence of $f(R)$ gravity with a scalar-tensor theory in the Einstein frame, where the scalaron $\phi$ effectively plays the role of the inflaton [27, 28].

Inflation takes place when the scalaron’s potential energy $V(\phi)$ dominates over its kinetic energy, resulting in a nearly exponential expansion of the scale factor. This regime is achieved under the so-called slow-roll conditions, where the evolution of the scalar field is sufficiently gradual. The slow-roll parameters are defined as

[TABLE]

where $V^{\prime}(\phi)=dV/d\phi$ and $V^{\prime\prime}(\phi)=d^{2}V/d\phi^{2}$. The slow-roll conditions $\epsilon\ll 1$ and $|\eta|\ll 1$ ensure a sustained period of quasi-de Sitter expansion [29, 30].

During inflation, the scale factor evolves approximately as $a\approx\exp(N)$, where the number of e-folds $N$ quantifies the total amount of expansion. It is computed as

[TABLE]

with $\phi^{*}$ and $\phi_{\mathrm{end}}$ denoting the values of the scalaron field at horizon exit and at the end of inflation, respectively. For a successful inflationary scenario consistent with cosmological observations, the number of e-folds is typically required to lie within the range $55\leq N\leq 65$.

The inflationary model’s predictions can be characterized by key observables, including the spectral index $n_{s}$, the tensor-to-scalar ratio $r$, and the amplitude of scalar perturbations $A_{s}$ at horizon crossing. These quantities are expressed in terms of the slow-roll parameters as

[TABLE]

To be viable, these predictions must align with current cosmological constraints, such as those reported by the Planck 2018 results: $n_{s}=0.9649\pm 0.0042$, $r<0.036$, and $A_{s}\approx 2.1\times 10^{-9}$ [1]. New measurements from the Atacama Cosmology Telescope (ACT) indicate a higher scalar spectral index ($n_{s}$) compared to previous Planck results. The combined analysis of ACT and Planck data yields $n_{s}=0.9709\pm 0.0038$. This value further increases to $n_{s}=0.9743\pm 0.0034$ when additional cosmological datasets are included, incorporating CMB observations, baryon acoustic oscillations (BAO), and Dark Energy Spectroscopic Instrument (DESI) measurements [31, 32].

2.3 Viable Dark energy Models

For $f(R)$ gravity to serve as a viable framework for explaining late-time cosmic acceleration and the nature of dark energy, it must satisfy a set of theoretical and observational viability conditions. These constraints ensure the physical consistency of the theory and its compatibility with cosmological and astrophysical data. The principal requirements are as follows [33, 5]. First, the effective gravitational coupling must remain positive, which requires $f_{R}>0$, thereby guaranteeing that gravity remains attractive. Second, the stability of scalar perturbations demands that $f_{RR}>0$, avoiding ghost-like instabilities and ensuring the stability of the scalar degree of freedom. Third, the theory must asymptotically approach the $\Lambda$CDM model at high curvatures, such that $f(R)\rightarrow R-2\Lambda$, in order to reproduce the correct cosmological behavior at early times. Fourth, the existence of a stable de Sitter solution at low curvature is necessary to account for the present accelerated expansion. Finally, a viable $f(R)$ model must incorporate a chameleon mechanism to suppress deviations from General Relativity in high-density environments, thereby satisfying local gravity constraints such as those tested in the Solar System.

Models that fulfill these criteria are referred to as viable $f(R)$ gravity models. Based on their asymptotic behavior, these models can generally be categorized into two classes: power-law type and exponential type. Representative forms of each class are given by [9, 4, 34, 16, 35, 36]:

[TABLE]

where $\lambda_{s}$, $\lambda_{e}$, and $n$ are dimensionless parameters, and $R_{ch}$ denotes a characteristic curvature scale. In the high-curvature limit, both models approach $f(R)\rightarrow R-\lambda R_{ch}$, effectively mimicking the $\Lambda$CDM model. In this limit, $\lambda R_{ch}$ plays the role of an effective cosmological constant, and $R_{ch}$ is typically of the order $R_{ch}\sim\Lambda_{obs}/\lambda$, where $\Lambda_{obs}=\kappa^{2}\rho_{DE}^{(0)}$ is defined by the present dark energy density. Consequently, a smaller $R_{ch}$ is required for larger values of $\lambda$ to match observational constraints, leading to a larger ratio $R/R_{ch}$. In the limit where $\lambda\gg 1$ and $R_{ch}\rightarrow 0$, $f(R)$ gravity becomes indistinguishable from the $\Lambda$CDM model.

In this context, the Hu-Sawicki model [37] and the Starobinsky model [9] are examples of the power-law class, while the Tsujikawa model [4], the exponential gravity models [34, 16, 35, 5, 38], and the Appleby-Battye (AB) model [39, 12] fall into the exponential-type category.

The viability of these $f(R)$ gravity models has far-reaching implications for cosmology. They provide an alternative to the cosmological constant for explaining the accelerated expansion of the universe, relying instead on a dynamical mechanism that naturally emerges from the modified gravitational action. Furthermore, the presence of the chameleon mechanism allows $f(R)$ theories to satisfy stringent local gravitational constraints, enhancing their applicability across both cosmological and astrophysical scales. Future observations, particularly those involving gravitational waves and large-scale structure surveys, are expected to offer critical tests for distinguishing $f(R)$ gravity from other modified gravity frameworks [7].

In the high-curvature regime, where the scalar curvature $R$ is much larger than the characteristic scale $R_{ch}$, the functional forms of the power-law and exponential models simplify as follows:

[TABLE]

These simplified forms will be utilized in the subsequent analysis to study the emergence of curvature singularities in dense astrophysical environments within the context of $f(R)$ dark energy models.

3 Curvature singularity in $f(R)$ gravity dark energy model

Starting from the field equations given in Eq. 2, one can derive the trace equation by contracting with the metric tensor. This yields:

[TABLE]

where $T$ is the trace of the matter energy-momentum tensor $T_{\mu\nu}^{(M)}$. This equation reduces to the standard trace equation in General Relativity in the limit $f(R)=0$.

In a localized region of spacetime that can be approximated as a perfect, homogeneous, and isotropic fluid, and where the non-relativistic matter density $\rho_{m}$ is significantly larger than the background density $\rho_{ch}$, the trace of the energy-momentum tensor can be assumed to evolve homogeneously with time. This behavior is modeled as [10, 14]:

[TABLE]

where $T_{0}=\rho_{m}^{(0)}-3P_{m}^{(0)}\simeq\rho_{m}^{(0)}$ is the initial local energy density and $t_{ch}$ represents the characteristic timescale over which this density evolves.

Assuming that the gravitational field of the object is sufficiently weak, the background spacetime can be approximated by a flat Minkowski metric. Under this approximation, the d’Alembert operator in Eq. 14 reduces to a second-order time derivative, thereby simplifying the trace equation to a more tractable form:

[TABLE]

To further facilitate the analysis, we introduce a change of variables from $R$ and $t$ to dimensionless variables $y=\beta(R_{\text{ch}}/R)$ and $\tau=\gamma^{-1}t$, where $\beta$ and $\gamma$ are scaling parameters. This transformation allows the time evolution of the system to be expressed in terms of $\tau$, with derivatives denoted by a prime (′). It is important to emphasize that the functional forms of $y$ and $\gamma$ depend on the specific $f(R)$ model considered. To reflect this model dependence, subscripts are used accordingly.

For the power-law type $f(R)$ dark energy model, the dimensionless variable is defined as $y_{p}=\beta(R_{\text{ch}}/R)$. The corresponding equation of motion takes the following form:

[TABLE]

where the scaling parameters are given by

[TABLE]

Substituting these into Eq. (17), the equation simplifies to:

[TABLE]

where $\tau_{\text{ch}}=\gamma_{p}^{-1}t_{\text{ch}}$ represents the rescaled characteristic timescale.

The parameter $\beta$ is chosen to satisfy $\beta=\kappa^{2}T_{0}/R_{\text{ch}}\simeq\rho_{m}/\rho_{c}\gg 1$, ensuring that $y$ remains finite and physically meaningful in high-density environments. For instance, for a matter density of $\rho_{m}\sim 10^{-24}\,\text{g/cm}^{3}$ and a characteristic curvature scale $R_{\text{ch}}\propto R_{0}\sim 10^{-35}\,\text{s}^{-2}$, we find $\beta\sim 1.67\times 10^{5}$. Assuming $\lambda\simeq 1$, the scaling parameter $\gamma_{p}$ evaluates to approximately $\gamma_{p}\sim 519\,\text{s}$ for $n=2$ and $\gamma_{p}\sim 0.004\,\text{s}$ for $n=3$.

Equation (19) can be solved numerically using the initial conditions $y_{p}(0)=\kappa^{2}T_{0}/R\simeq 1$ and $y_{p}^{\prime}(0)=0$. These conditions model the initial configuration of the high-density environment, such as during gravitational collapse. The differential equation is nonlinear due to the terms $2ny_{p}^{\prime 2}/y_{p}$ and $y_{p}^{-2n}$, precluding closed-form solutions. Nevertheless, numerical solutions reveal that the system exhibits oscillatory behavior. The equilibrium trajectory of the system can be determined by setting $y_{p}^{\prime\prime}=0$ and $y_{p}^{\prime}=0$. Under these conditions, the equilibrium line is found to be $y(\tau)=1/(1+\tau/\tau_{\text{ch}})$, which is independent of $n$ and asymptotically approaches zero as $\tau\to\infty$. This equilibrium curve is included in the graphical analysis. The oscillatory dynamics can be interpreted by examining the roles of the nonlinear terms. The term $2ny_{p}^{\prime 2}/y_{p}$ is always positive for $y_{p}\geq 0$ and $n>0$, thus acting as a damping component that drives $y_{p}^{\prime\prime}$ negative. Meanwhile, the term $y_{p}^{-2n}\left[\left(1+\frac{\tau}{\tau_{\text{ch}}}\right)-y_{p}^{-1}\right]$ acts as a time-dependent restoring force, driving the system toward the equilibrium configuration whenever deviations occur.

Numerical results for Eq. 19 are illustrated in Fig. 1, showing that the oscillations of $y_{p}$ are centered around the equilibrium curve $y=1/(1+\tau/\tau_{\text{ch}})$. As the value of $n$ increases, the frequency of oscillation decreases. Since $y_{p}\propto 1/R$, a curvature singularity emerges as $y_{p}\to 0$. The figure indicates that the power-law dark energy model leads to curvature singularities at finite values of $\tau/\tau_{\text{ch}}$, depending on the value of $\tau_{\text{ch}}$. These singularities occur at later times for larger $\tau_{\text{ch}}$. For example, when $\tau_{\text{ch}}=20$, the singularity occurs at $\tau/\tau_{\text{ch}}\approx 0.7434$ for $n=2$ and $\tau/\tau_{\text{ch}}\approx 0.3853$ for $n=3$. In contrast, for $\tau_{\text{ch}}=80$, the singularities appear at $\tau/\tau_{\text{ch}}\approx 1.4927$ for $n=2$ and $\tau/\tau_{\text{ch}}\approx 0.7722$ for $n=3$. These results suggest that in high-density environments, curvature singularities can develop within a relatively short timescale, given by $t_{\text{sing}}=\gamma_{p}\tau=\alpha\gamma_{p}\tau_{\text{ch}}$, where $\alpha$ is a constant of order unity.

For the exponential-type $f(R)$ dark energy model, Eq. 13 is substituted into Eq. 16. To distinguish this model from the power-law case, we define the dimensionless variable $y_{e}=\beta(R_{\text{ch}}/R)$ and rescale the time variable using $\gamma_{e}$, yielding $\tau=\gamma_{e}^{-1}t$. With these definitions, the governing differential equation becomes

[TABLE]

where

[TABLE]

To analyze this equation, we impose initial conditions that represent the physical state at the onset of the high-density regime. Specifically, we set $y_{e}(0)=\kappa^{2}T_{0}/R$, normalized to unity, and $y_{e}^{\prime}(0)=0$. These conditions ensure that the system begins in equilibrium prior to the emergence of curvature oscillations. Numerical solutions of Eq. 20 are shown in Fig. 2. The dynamics are governed by the interplay between a damping term, $(\beta y_{e}^{-1}-2)y_{e}^{\prime 2}/y_{e}$, and a restoring term, $y_{e}^{2}e^{\beta y_{e}^{-1}}\left[\left(1+\tau/\tau_{\text{ch}}\right)-y_{e}^{-1}\right]$. The damping term, which scales with $\beta$, gradually suppresses the oscillations. For values of $\beta$ in the range of 10–20, the damping is moderate, allowing oscillations to persist for some time before decay sets in.

The equilibrium configuration of the system is defined by the trajectory $y_{e}(\tau)=1/(1+\tau/\tau_{\text{ch}})$. As $\tau\to\infty$, this trajectory asymptotically approaches zero, implying that the Ricci scalar $R$ diverges and a curvature singularity emerges. However, the numerical results show that the system exhibits damped oscillations around this equilibrium path before reaching the singular state. The amplitude and frequency of these oscillations depend on both $\beta$ and the characteristic timescale $\tau_{\text{ch}}$.

The behavior of the exponential-type model stands in contrast to the power-law $f(R)$ model, where oscillations are generally more pronounced and the singularity tends to occur later for larger values of $\tau_{\text{ch}}$. This distinction arises due to the presence of the exponential term $e^{\beta y_{e}^{-1}}$ in the restoring force, which significantly enhances the system’s sensitivity to variations in $R$. For moderate values of $\beta$ (in the range 10–20), representing a curvature regime that is not excessively high, the exponential model exhibits stronger damping and earlier singularity onset compared to the power-law model.

The parameter $\gamma_{e}$ plays a central role in setting the temporal scale of the system’s evolution. Given its definition, $\gamma_{e}^{2}=3\lambda/R_{\text{ch}}$, and using representative values $\lambda\sim 1$ and $R_{\text{ch}}\sim 10^{-35}\,\text{s}^{-2}$, we estimate $\gamma_{e}\sim\sqrt{3\times 10^{35}}\,\text{s}\approx 1.73\times 10^{17}\,\text{s}$. The singularity in the exponential-type model arises as $y_{e}\to 0$, corresponding to $R\to\infty$. As shown in Fig. 2, this occurs at finite values of $\tau/\tau_{\text{ch}}$, with the precise timing dependent on the parameters $\beta$ and $\tau_{\text{ch}}$. For example, when $\tau_{\text{ch}}=0.1$, the singularity occurs at $\tau/\tau_{\text{ch}}\approx 0.2048$ for $\beta=12$. For $\tau_{\text{ch}}=0.5$, it occurs at $\tau/\tau_{\text{ch}}\approx 0.4146$ for $\beta=17$. These findings imply that the exponential-type $f(R)$ dark energy model predicts curvature singularities within a finite, though relatively long, timescale given by $t_{\text{sing}}=\gamma_{e}\tau=\alpha\gamma_{e}\tau_{\text{ch}}$, where $\alpha\in\mathcal{O}(10^{-1})$. Based on the numerical simulations, we estimate the singularity time to lie within the range $t_{\text{sing}}=[10^{15}-10^{16}]\,\text{s}$, which remains below the current age of the universe, $t_{U}\simeq 4\times 10^{17}\,\text{s}$.

4 Solving singularity by adding $R^{2}$ term

Inflationary dynamics within the framework of $F(R)$ gravity have been the subject of extensive investigation over several decades [22, 40, 41, 42, 43, 13, 44, 45]. Among the various models proposed, the Starobinsky model, characterized by the functional form [22]:

[TABLE]

stands out as one of the most successful and widely studied. This model, despite its simplicity, exhibits exceptional predictive accuracy in describing inflationary phenomena. The Starobinsky model provides well-defined predictions for key inflationary observables, such as the scalar spectral index $n_{s}$ and the tensor-to-scalar ratio $r$, which are expressed as:

[TABLE]

For typical values of the e-folding number $N\approx 50-60$, these predictions yield:

[TABLE]

which are in excellent agreement with observational data from the Planck satellite. The inflaton mass scale $M\approx 1.13\times 10^{-5}M_{\text{Pl}}$ is determined by the amplitude of scalar perturbations. Remarkably, the predictions for $n_{s}$ and $r$ remain robust against variations in $M$, demonstrating the model’s stability across different parameter regimes.

Recent joint analyses incorporating data from Planck, ACT, and DESI yield a scalar spectral index of $n_{s}=0.974\pm 0.003$ [31, 32]. Notably, this value positions the Starobinsky inflationary model at the boundary of the $2\sigma$ confidence interval derived from the combined observational constraints, under the assumption of approximately 60 e-folds of inflation[32]. This proximity has sparked renewed interest and debate within the cosmology community, prompting critical reassessments of canonical inflationary models and motivating the exploration of alternative theoretical frameworks in light of the high-precision measurements from ACT [46, 47, 48, 49]. Nevertheless, a more refined treatment of the Starobinsky scenario, as presented in [50], demonstrates that the model remains compatible with current observational data at the $2\sigma$ level for e-folding numbers $N\geq 58$.

Regarding curvature singularities, studies [10, 11] have shown they can be avoided by introducing additional terms $R^{m}/M^{2(m-1)}$ with $1<m\leq 2$ in $f(R)$ dark energy models. The Starobinsky model itself has been proven to prevent such singularities [10, 12].

In the high-curvature regime $R\gg R_{\rm ch}$, the power-law $f(R)$ dark energy model with Starobinsky correction becomes:

[TABLE]

where $m_{R}$ represents an adjustable mass scale. The corresponding differential equation derived from Eq. 14 is:

[TABLE]

with the coefficient:

[TABLE]

Here we maintain the rescaling parameters $\beta$ and $\gamma_{p}$ used in deriving Eq. 17.

The inclusion of the $R^{2}$ correction term modifies the evolution of the curvature scalar $R$ in dense environments. Equation 26 governs the dynamics of the rescaled curvature variable $y_{p}=\beta(R_{\mathrm{ch}}/R)$ and incorporates the higher-order contribution from the Starobinsky term. Relative to Equation 17, the primary distinction lies in the additional term $g_{pI}y_{p}^{-(2n+2)}\left[y_{p}^{\prime\prime}-\frac{2y_{p}^{\prime 2}}{y_{p}}\right]$, which directly originates from the $R^{2}$ component.

To understand its physical impact, recall that Equation 17 describes an oscillatory system driven by a damping term $2n\frac{y_{p}^{\prime 2}}{y_{p}}$ and a restoring force $y_{p}^{-2n}\left[\left(1+\frac{\tau}{\tau_{\mathrm{ch}}}\right)-y_{p}^{-1}\right]$. This system oscillates around the attractor solution $y_{p}(\tau)=1/(1+\tau/\tau_{\mathrm{ch}})$, with singularities appearing when $y_{p}\to 0$ (i.e., $R\to\infty$). The $R^{2}$ correction modifies this evolution by introducing a feedback mechanism that couples the second derivative $y_{p}^{\prime\prime}$ and the nonlinear term $y_{p}^{\prime 2}/y_{p}$ via a factor $g_{pI}y_{p}^{-(2n+2)}$. Since the exponent is negative, this term becomes increasingly dominant as $y_{p}$ decreases—precisely in the regime where singularities are most threatening. Physically, the correction acts as a stabilizer: when $y_{p}^{\prime\prime}$ is positive, the opposing term $-2y_{p}^{\prime 2}/y_{p}$ can reduce the rate of divergence, dampening the trajectory toward $y_{p}=0$. Conversely, when $y_{p}^{\prime\prime}$ is negative, the term reinforces the restoring effect, countering further deviation. As a result, the $R^{2}$ term serves as a regulator that may prevent the onset of curvature singularities.

The magnitude of this regulation depends critically on the parameter $g_{pI}$. When $g_{pI}$ is small, the dynamics remain governed by the original oscillatory terms, allowing singularities to develop. As $g_{pI}$ increases, the influence of the stabilizing term grows, reducing the oscillation amplitude and potentially preventing $y_{p}$ from reaching zero. This reveals a clear threshold effect: for sufficiently large $g_{pI}$, the curvature evolution is regularized.

Numerical solutions to Equation 26 are displayed in Figure 3, where $n=2$ and $\tau_{\mathrm{ch}}=80$, enabling direct comparison with Figure 1. The figure highlights the impact of varying $g_{pI}$. For $g_{pI}\ll 1$, the curvature scalar undergoes pronounced oscillations, and the singularity time $\tau_{\mathrm{sing}}$ is highly sensitive to $g_{pI}$. Specifically, $\tau_{\mathrm{sing}}\simeq 505.76\tau_{\mathrm{ch}}$ for $g_{pI}=10^{-4}$, $\tau_{\mathrm{sing}}\simeq 55.76\tau_{\mathrm{ch}}$ for $g_{pI}=10^{-8}$, and $\tau_{\mathrm{sing}}\simeq 4.33\tau_{\mathrm{ch}}$ for $g_{pI}=10^{-12}$. As $g_{pI}$ increases, the oscillations are progressively damped, and $y_{p}(\tau)$ approaches the equilibrium trajectory $y=1/(1+\tau/\tau_{\mathrm{ch}})$. When $g_{pI}>1$, singularities are fully suppressed.

These results confirm that the $R^{2}$ correction term can effectively resolve curvature singularities in power-law $f(R)$ dark energy models, provided that $g_{pI}>1$. For $\lambda=1$, $R_{\mathrm{ch}}\simeq 4\Lambda\simeq 2.024\times 10^{-83}M_{\mathrm{Pl}}^{2}$, and $m_{R}=M$, one finds $g_{pI}\simeq 10^{-44}$ for $n=2$, $g_{pI}\simeq 10^{-34}$ for $n=3$, and $g_{pI}\simeq 10^{-24}$ for $n=4$. These values indicate that the Starobinsky inflationary scale, $M\approx 1.13\times 10^{-5}M_{\mathrm{Pl}}$, is too large to suppress singularities in this context. Achieving $g_{pI}>1$ requires lower mass scales: $m_{R}\simeq 1.28\times 10^{-27}M_{\mathrm{Pl}}$ for $n=2$, $m_{R}\simeq 1.26\times 10^{-22}M_{\mathrm{Pl}}$ for $n=3$, and $m_{R}\simeq 1.37\times 10^{-17}M_{\mathrm{Pl}}$ for $n=4$. These findings underscore the crucial role of the $R^{2}$ term in stabilizing curvature dynamics and preventing singularity formation in modified gravity theories.

In the regime where $R\geq R_{\mathrm{ch}}$, the exponential-law model takes the form

[TABLE]

where the final term represents the $R^{2}$ correction. The corresponding differential equation derived from Equation 14 is given by

[TABLE]

where the coefficient $g_{eI}$ is defined as

[TABLE]

Here, $\beta$ and $\gamma_{e}$ retain their definitions as introduced in the uncorrected model.

The third term in Equation 29, given by $g_{eI}e^{\beta(1/y_{e}-1)}\left(y_{e}^{\prime\prime}-\frac{2y_{e}^{\prime 2}}{y_{e}}\right)$, originates from the $R^{2}$ correction applied to the exponential-type $f(R)$ dark energy model. This term plays a crucial role in regulating the curvature dynamics and mitigating singular behavior through several interrelated mechanisms. The exponential prefactor $e^{\beta(1/y_{e}-1)}$ renders this term highly sensitive to the curvature regime: as $y_{e}\to 0$ (corresponding to $R\to\infty$), the term increases rapidly, thereby dominating the evolution near potential singularities. The differential structure $y_{e}^{\prime\prime}-2y_{e}^{\prime 2}/y_{e}$ acts as a feedback mechanism. When $y_{e}^{\prime\prime}>0$, indicating an accelerating collapse toward a singularity, the negative contribution $-2y_{e}^{\prime 2}/y_{e}$ can counterbalance this acceleration. Conversely, when $y_{e}^{\prime\prime}<0$, the term can reinforce the restoring force, inhibiting divergent behavior. This response contrasts with the polynomial prefactor in the power-law model of Equation 26, where the exponential term here enables stronger and more localized regulation of curvature growth, even for relatively small values of the coupling constant $g_{eI}$.

The numerical solutions illustrated in Figure 4 elucidate the impact of the $R^{2}$ correction on the system’s behavior. For small values of $g_{eI}$ ($g_{eI}\ll 1$), the evolution of $y_{e}(\tau)$ exhibits damped oscillations around the equilibrium trajectory $y_{e}=1/(1+\tau/\tau_{\rm ch})$. Compared to the uncorrected scenario, the oscillation amplitudes are significantly suppressed by the exponential correction. As $g_{eI}$ exceeds unity, the singular behavior is entirely suppressed: the function $y_{e}(\tau)$ asymptotically approaches the equilibrium solution without ever reaching zero, indicating complete stabilization of the curvature scalar.

This parameter dependence yields important physical insights. In particular, achieving $g_{eI}\geq 1$ necessitates the mass scale $m_{R}$ of the $R^{2}$ correction to be on the order of $10^{-38}M_{\rm Pl}$, assuming typical values of $\lambda\sim 1$ and $R_{\rm ch}\sim 4\Lambda$. This required scale is dramatically lower than the inflationary mass scale associated with the Starobinsky model, $M\sim 10^{-5}M_{\rm Pl}$. Thus, these findings suggest that effective singularity resolution via $R^{2}$ corrections in the exponential $f(R)$ model operates at energy scales far below those characteristic of early-universe inflation.

5 Solving singularity by adding $R^{\frac{m+2}{m+1}}$ term

We examine whether $f(R)$ models capable of driving inflation can simultaneously avoid curvature singularities in the context of $f(R)$ dark energy. Our analysis focuses on a generalized $f(R)$ gravity framework that extends the Starobinsky scalar potential [13], specifically considering the case where the potential takes the form $V(\phi)\simeq\gamma\exp(-m\phi)$ with $m>0$. In the high-curvature regime ($R\gg\gamma$), this yields:

[TABLE]

The corresponding primordial power spectrum amplitude, spectral indices, and tensor-to-scalar ratio are given by:

[TABLE]

These predictions remain consistent with Planck data [1] for $m=0.1,0.2$, with $\gamma\simeq 10^{-11}M_{\rm Pl}^{2}$, indicating that Starobinsky-like $R^{2}$ gravity models can successfully drive this inflationary scenario.

For the power-law $f(R)$ dark energy model with $R^{\frac{m+2}{m+1}}$ correction in the limit $R\gg R_{\rm ch}$, we obtain:

[TABLE]

where

[TABLE]

Here, $m$ represents an arbitrary positive real number, not necessarily integer-valued. The correction term $R^{\frac{m+2}{m+1}}$ approaches $R^{2}$ for small $m$ and converges to $R$ for large $m$. From this, we derive the differential equation:

[TABLE]

with

[TABLE]

Notably, Eq. (35) reduces to Eq. (25), while Eq. (36) converges to Eq. (27) as $m\rightarrow 0$.

Following our established methodology, we numerically solve the differential equation for $y_{p}(\tau)$ using initial conditions $y_{p}(0)=1$ and $y_{p}^{\prime}(0)=0$. Figure 5 presents these solutions for various $m$ and $g_{pII}$ values in the high-density regime ($\beta\simeq 10^{5}$), with fixed parameters $n=2$ and $\tau_{ch}=80$ to enable direct comparison with the $R^{2}$ correction case shown in Fig. 3.

The results reveal several key features. First, the power-law parameter $m$ exhibits relatively weak influence on the solutions, with minimal differences between small ($m=0.1$) and large ($m=10$) values. The $m=0.1$ case closely resembles the $R^{2}$ correction results, confirming that the $R^{2}$ scenario represents a special case of our more general $R^{\frac{m+2}{m+1}}$ framework.

While $m$ appears primarily relevant for inflationary-scale corrections, we find that curvature singularities can be avoided even for large $m$ values. However, larger $m$ values do enhance both the amplitude and frequency of $y_{p}(\tau)$ oscillations. The parameter $g_{pII}$ emerges as the dominant factor controlling singularity formation - when $g_{pII}\geq 1$, both oscillation amplitude and frequency decrease significantly, with $y_{p}(\tau)$ asymptotically approaching the equilibrium solution $y=1/(1+\tau/\tau_{ch})$. This behavior robustly prevents curvature singularities at finite times.

The exponential-law $f(R)$ dark energy model incorporating an $R^{\frac{m+2}{m+1}}$ correction term takes the following form:

[TABLE]

which yields the corresponding differential equation derived from Eq. 14:

[TABLE]

where the parameter $g_{eII}$ is defined as:

[TABLE]

In the limit $m\rightarrow 0$, Eqs. (5) and (39) reduce to Eqs. (29) and (30), respectively, recovering the $R^{2}$ correction case.

Figure 6 presents numerical solutions to the differential equation (5) for model parameters $\beta=17$ and $\tau_{ch}=0.05$, chosen to enable direct comparison with the $R^{2}$ correction results shown in Figure 4. Following the same convention as in the power-law dark energy analysis, we display results for small $m$ values ($m=0.1$) in the left panels (a, c, e) and large $m$ values ($m=10$) in the right panels (b, d, f).

For the small $m$ case ($m=0.1$), the numerical solutions show excellent agreement with the $R^{2}$ correction results from Figure 4. The differential term containing the $g_{eII}$ factor proves particularly effective at regularizing curvature singularities that would otherwise occur when $y_{e}(\tau)\rightarrow 0$. Furthermore, when $g_{eII}\geq 1$, we observe significant suppression of both the oscillation frequency and amplitude of $y_{e}(\tau)$.

In the large-$m$ regime ($m=10$), where the term $R^{\frac{m+2}{m+1}}$ asymptotically approaches $R^{n}$ with $n\sim 1$, we observe qualitatively similar behavior. The $g_{eII}$ term continues to effectively regularize curvature singularities while suppressing both the amplitude and frequency of oscillations when $g_{eII}\geq 1$. In contrast, for $g_{eII}<1$, the system exhibits amplified oscillatory behavior, with increased amplitude and frequency, potentially resulting in finite-$\tau$ singularities that arise more rapidly than in the small-$m$ case. These findings demonstrate that $R^{n}$ corrections within the range $1<n\leq 2$ can successfully mitigate curvature singularities in exponential-law dark energy models, provided the condition $g_{eII}\geq 1$ is satisfied.

6 Solving singularity by adding $\alpha$-attractor term

Among inflationary models proposed to address the initial singularity problem, the $\alpha$-attractor framework stands out for its capacity to unify diverse inflationary scenarios, including the well-established Starobinsky model [51, 52, 53, 54]. This class of models comprises two principal variants: (1) the T-model with potential $V_{T}\propto\tanh^{2}(\varphi/\sqrt{6\alpha})$, and (2) the E-model characterized by $V_{E}\propto(1-e^{-\sqrt{2/3\alpha}\varphi})^{2}$. Remarkably, despite their differing functional forms, both variants yield nearly identical predictions for inflationary observables.

The $\alpha$-attractor models generate universal predictions for the scalar spectral index $n_{s}$ and tensor-to-scalar ratio $r$:

[TABLE]

where $\alpha>0$ represents a dimensionless parameter controlling model flexibility, and $N$ denotes the required e-folding number. When $\alpha=1$, these predictions align precisely with those of the Starobinsky and Higgs inflation models, while the $\alpha\to\infty$ limit recovers chaotic quadratic inflation.

Within $f(R)$ gravity, the $\alpha$-attractor model corresponds to a polynomial modification of the Starobinsky model in the high-curvature regime [55]:

[TABLE]

where

[TABLE]

Analysis of $n_{s}$ and $r$ as functions of $\zeta$ reveals excellent consistency with Planck observational constraints for $-0.02<\zeta\lesssim 0.8$ (corresponding to $0.96<\alpha\lesssim 25$). The upper limit $\zeta\approx 0.8$ corresponds to an $R^{1.8}$ correction to the Starobinsky potential, demonstrating the model’s capacity to encompass diverse inflationary dynamics. A recent analysis in [56] investigated the impact of adding an $R^{3/2}$ correction to the Starobinsky model. The study found that this extension improves agreement with ACT data, without affecting the model’s consistency with other observational results. Interestingly, this type of modification is closely related to the $\alpha$-attractor class, suggesting that such models are worth further theoretical and phenomenological study as possible refinements of standard inflationary theory.

The power-law dark energy model incorporating $\alpha$-attractor corrections takes the form:

[TABLE]

The associated differential equation derived from Eq. (14) becomes:

[TABLE]

with coefficients defined as:

[TABLE]

This differential equation reveals how $\alpha$-attractor corrections modify curvature oscillations through terms containing $g_{pIIIa}$ when $\zeta\neq 1$.

Numerical solutions of Eq. (6) with initial conditions $y_{p}(0)=1$ and $y^{\prime}_{p}(0)=0$ are presented in Figure 7 for various parameter combinations ($n,\zeta,\tau_{ch},\beta,g_{pIIIa}$). Throughout our analysis, we maintain $0<g_{pIIIa}\leq 1$ and $\beta\gg 1$ to emphasize the $\alpha$-attractor correction’s influence. Parameter dependencies are organized as: $\zeta$ in panels (a) and (d), $g_{pIIIa}$ in (b) and (e), and $g_{pIIIb}\simeq 1$ in (c) and (f).

Panels (a) and (d) reveal that significant curvature fluctuations occur primarily for $0<\zeta\ll 1$, with amplitude decreasing sharply as $\zeta$ increases. This suggests that $\alpha$-attractor terms with $\zeta\gtrsim 0.3$ (equivalent to Starobinsky inflation with $R^{1.3}$ corrections) effectively suppress both curvature fluctuations and singularities in the power-law dark energy model. Panels (b) and (e) demonstrate that smaller $g_{pIIIa}$ values reduce fluctuation amplitudes while increasing their frequency. When $0<g_{pIIIa}\ll 1$ (yielding $g_{pIIIb}\simeq 1$), the solution $y_{p}(\tau)$ closely follows the equilibrium curve $y=1/(1+\tau/\tau_{ch})$ for all $\beta\gg 1$ cases. The condition $g_{pIIIb}\simeq 1$ emerges as crucial for avoiding curvature singularities at arbitrary $\beta\gg 1$, allowing us to express the mass parameter as $M_{\alpha}^{2}/R_{c}=f(\beta,\zeta,n,\lambda)$. This dependence makes $M_{\alpha}^{2}$ particularly sensitive to the local density parameter $\beta$, resulting in values distinct from those constrained purely by inflationary considerations.

The exponential-law dark energy model with $\alpha$-attractor corrections is given by:

[TABLE]

with corresponding curvature fluctuation equation:

[TABLE]

where

[TABLE]

Similar to the power-law dark energy case with $\alpha$-attractor corrections, the exponential-law dark energy model exhibits modifications to curvature oscillations through the $g_{eIIIa}$ term when $\zeta\neq 1$. Figure 8 displays the numerical solutions for $y_{e}(\tau)$ obtained from Eq. (6) across various parameter combinations of $\beta$, $\zeta$, $\tau_{ch}$, $\lambda$, and $g_{eIIIa}$. The left panels (a) and (d) illustrate the dependence on $\zeta$, while panels (b) and (e) demonstrate the effects of varying $g_{eIIIa}$ within the range $0<g_{eIIIa}<1$. The right panels (c) and (f) present results for different $\beta$ values under the extreme condition $g_{eIIIa}=10^{-320}$, which corresponds to $g_{eIIIb}$ values between $10^{2}$ and $10^{4}$. Notably, when $g_{eIIIa}$ approaches such extremely small values, the solutions begin to develop curvature singularities.

The numerical results reveal several important features of the model. First, the amplitude of curvature oscillations in the exponential-law dark energy model decreases systematically with increasing $\zeta$ within the range $0<\zeta<1$. Remarkably, the solution $y_{e}(\tau)$ asymptotically approaches the equilibrium form $1/(1+\tau/\tau_{ch})$ even for very small $\zeta$ values ($\zeta=10^{-2}$). The solution profile further improves with larger values of the density parameter $\beta$ and smaller values of $g_{eIIIa}$, the latter leading to $g_{eIIIb}>10^{2}$. This behavior indicates that the $R^{1+\zeta}$ correction originating from the Starobinsky model - which emerges as a consequence of the $\alpha$-attractor effect - simultaneously enhances both the oscillation profile and the curvature singularity behavior in the exponential-law dark energy model. This improvement persists even for $\zeta\ll 1$, provided the condition $g_{eIIIa}\ll 1$ maintains $g_{eIIIb}>10^{2}$.

Further numerical investigation shows that solution profiles qualitatively similar to panel (d) can be obtained for parameter combinations such as $\beta=30$ with $\zeta=0.4$ and $g_{eIIIa}\simeq 1.0$. However, when $\beta$ is reduced to values like $\beta=17$ while maintaining $\zeta=0.4$ and $g_{eIIIa}\simeq 1.0$, significant curvature fluctuations emerge, producing profiles resembling panel (c). This demonstrates that when $g_{eIIIa}\simeq 1.0$, proper selection of $\zeta$ within $0<\zeta<1$ can effectively suppress both singularities and curvature fluctuations in the exponential-law dark energy model.

7 Conclusion

Dark energy models formulated within $f(R)$ gravity frameworks inherently face the challenge of curvature singularities developing in high-density regimes. This pathological behavior emerges when examining the trace equation derived from $f(R)$ field equations, which governs the dynamics of curvature fluctuations. The singularity becomes apparent as the scalar curvature $R$ diverges ($R\to\infty$), mathematically equivalent to the condition where the rescaled variable $y\equiv\beta(R_{\text{ch}}/R)$ approaches zero. Numerical solutions of this system reveal that $y$ exhibits characteristic oscillations around the equilibrium trajectory $y=1/(1+\tau/\tau_{\text{ch}})$, where $\tau$ represents the rescaled time coordinate and $\tau_{\text{ch}}$ the rescaled characteristic timescale. These singularities manifest at finite times in both power-law and exponential-law $f(R)$ models when the local density parameter satisfies $\beta\gg 1$, with the physical singularity timescale $t_{\text{sing}}=\gamma_{p}\tau_{\text{sing}}$ falling within cosmologically relevant ranges of $10^{15}-10^{16}$ seconds.

To address this fundamental limitation, we implement curvature correction terms inspired by successful inflationary $f(R)$ models, specifically the Starobinsky $R^{2}$ term, its generalized extension $R^{\frac{m+2}{m+1}}$, and the $\alpha$-attractor polynomial extension. These corrections serve the dual purpose of resolving curvature singularities while preserving compatibility with inflationary dynamics, thereby offering a unified framework for cosmic acceleration across different epochs. The stabilization mechanism operates through distinct dimensionless parameters: the $R^{2}$ correction is governed by $g_{pI}$ for power-law models and $g_{eI}$ for exponential models, while the generalized extension introduces $g_{pII}$ and $g_{eII}$ parameters, and the $\alpha$-attractor modification involves $g_{pIIIa,b}$ for power-law models and $g_{eIIIa,b}$ for exponential models.

Comprehensive numerical simulations demonstrate that curvature singularities can be completely avoided when specific threshold conditions are met: $g_{pI,eI}\geq 1$ for the Starobinsky case, $g_{pII,eII}\geq 1$ for the generalized extension, and $g_{pIIIb}\approx 1$ and $g_{eIIIa}\approx 1$ for the $\alpha$-attractor model. These stability conditions necessitate inflationary parameters that substantially differ from conventional inflation-constrained values. Notably, the effective mass scale $m_{R}$ required for successful $R^{2}$ corrections must be approximately 33 orders of magnitude lower than the standard Starobinsky scale. The exponent $m$ in generalized extensions shows optimal performance when approaching the Starobinsky limit ($m\to 0$), while the $\alpha$-attractor parameter $\zeta=1-\frac{1}{\sqrt{\alpha}}$ demonstrates maximum efficacy within the range $0<\zeta<1$.

The analysis yields three principal insights regarding model construction. First, the Starobinsky $R^{2}$ term emerges as the most effective singularity resolution mechanism, particularly when other extensions asymptotically approach its functional form. Second, successful singularity avoidance requires the introduction of new mass scales distinct from traditional inflationary parameters. Third, polynomial extensions of the form $\propto R^{\zeta+1}+R^{2}$ demonstrate promising capacity to simultaneously describe early-universe inflation and late-time dark energy while maintaining mathematical regularity. These findings collectively establish that incorporating inflation-motivated corrections provides a viable pathway to singularity-free $f(R)$ dark energy models, with the necessary parameter adjustments preserving the models’ ability to unify the description of cosmic acceleration across different cosmological epochs.