Spherically Symmetric Potentials in Quadratic $f(R)$ Gravity

The quest to understand the fundamental nature of gravity has led to the development of numerous alternative theories, among which $f(R)$ gravity stands out for its mathematical elegance and phenomenological richness [1, 2, 3, 4, 5, 6, 7]. In these models, the gravitational action is generalized to depend on an arbitrary function of the Ricci scalar, $f(R)$, allowing for modifications to Einstein’s General Relativity (GR) that can address outstanding problems in cosmology and astrophysics [8, 9, 10, 11], such as the nature of dark energy [11, 12, 13, 14, 15, 16] and the dynamics of galactic rotation curves [17, 18, 19, 20]. Quadratic models represent a minimal extension with significant theoretical and observational consequences [21, 22, 23]. While the cosmological implications of quadratic $f(R)$ gravity have been extensively studied [24, 25, 26, 27, 28], much less attention has been paid to its impact on local gravitational potentials generated by realistic matter distributions.
In particular, static spherically symmetric sources—relevant for galaxies, dark matter halos, and compact objects—serve as ideal testbeds to explore the structure of the modified potential and its sensitivity to matter distributions. These systems provide a controlled environment where the interplay between gravity and the underlying mass profile can be studied in detail, offering insights into both fundamental theory and astrophysical phenomena. Standard mass profiles such as Plummer [29], Hernquist [30, 31], and Navarro–Frenk–White (NFW) [32, 33] have been widely used in astrophysics due to their empirical success in modeling luminous and dark matter components [34, 35], capturing the observed dynamics of stellar systems and the rotation curves of galaxies with remarkable accuracy [36, 37, 38]. Each of these profiles exhibits distinct inner and outer density behaviors, which in turn influence the resulting gravitational potential and its modifications under alternative gravity theories.
Despite their widespread application in classical gravity, a systematic analytical treatment of their gravitational potentials within modified gravity frameworks—such as quadratic $f(R)$ models or other extensions of general relativity—remains largely unexplored. Addressing this gap is crucial, as modifications to gravity can introduce new features in the potential, such as Yukawa-like corrections or altered asymptotic behavior, which may have observable consequences in galactic dynamics, gravitational lensing, and the structure of compact objects. A comprehensive analytical approach not only would clarify the theoretical implications of these modifications but also provide essential tools for interpreting astrophysical data and testing gravity beyond the Newtonian regime.

This paper is organized as follows. Section 2 focuses on deriving the exact solutions of the linearized fourth-order field equations in a quadratic $f(R)$ gravity model, under the assumption of static spherical symmetry. Section 3 provides an analytical solution to this equation, with integration constants determined by imposing appropriate physical conditions on the gravitational potential. In Section 4, the modified potential is computed and analyzed for various spherical mass distributions, including both classical profiles and newly proposed models. Finally, Section 5 discusses the main conclusions of this work.

We begin by considering $f(R)$ gravity as a generalization of the Einstein-Hilbert action in which the Ricci scalar $R$ is replaced by a function $f(R)$. This action takes the form

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

(1)

where $\kappa=8\pi$ in geometrized units, $g$ is the determinant of the metric tensor $g_{\mu\nu}$, and $S_{M}$ denotes the matter and energy contribution. Varying this action with respect to the metric, yields the modified field equations [2, 3]

$$FR_{\mu\nu}-\frac{1}{2}fg_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}F+g_{\mu\nu}\Box F=\kappa T_{\mu\nu},$$

(2)

where $f=f(R)$, $F=df/dR$, $\nabla_{\mu}$ denotes the covariant derivative, $\Box=\nabla_{\sigma}\nabla^{\sigma}=g^{\sigma\rho}\nabla_{\sigma}\nabla_{\rho}$ is the d’Alembert operator, and $T_{\mu\nu}$ is the stress-energy tensor, defined by

$$T_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta S_{M}}{\delta g^{\mu\nu}}.$$

(3)

The terms $\nabla_{\mu}\nabla_{\nu}F$ and $\Box F$ in Eq. (2), introduce fourth-order derivatives of the metric, distinguishing $f(R)$ gravity from second-order GR. Tracing this equation reveals a dynamical scalar degree of freedom, $F$, encoding the modified gravitational interaction. In this work, we consider the specific functional form

$$f(R)=-\frac{4\pi}{\beta}\left(4\gamma+2\alpha^{2}R+R^{2}\right),$$

(4)

where $\alpha$, $\beta$ and $\gamma$ are constant parameters, which encompasses several important limits of gravitational theories. In the limit $\alpha^{2}\to\infty$, $\beta/\alpha^{2}=-8\pi$ and $\gamma/\beta=\Lambda/(8\pi)$, the action reduces to GR with a cosmological constant $\Lambda$. Moreover, setting $\gamma=0$ recovers pure GR. Interestingly, for $\alpha^{2}=3m^{2}$, $\beta/\alpha^{2}=-8\pi$, $\gamma=0$, the model corresponds to the well-known Starobinsky quadratic gravity model [21],

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

(5)

where the parameter $m$ sets the mass scale of the scalaron, the additional scalar degree of freedom arising from the $R^{2}$ term, which is associated with the energy scale of inflation at the early universe, and acts as a coupling constant for higher-order curvature effects. In particular, the measured amplitude of primordial curvature perturbations in the Cosmic Microwave Background (CMB) fixes $m\approx 10^{-5}M_{pl}$, (where $M_{pl}$ is the reduced Planck mass), consistent with Planck data [39, 40]. Now, using Eq. (4), field equations can be written as

$$\alpha^{2}G_{\mu\nu}-\left(\gamma+\frac{1}{4}R^{2}\right)g_{\mu\nu}+\left(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}+R_{\mu\nu}\right)R=-\beta T_{\mu\nu}.$$

(6)

where we have used the Einstein tensor

$$G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}.$$

(7)

To study weak-field gravity in this framework, we linearize Eq. (6) around a Minkowski background metric $\eta_{\mu\nu}=\text{diag}(-1,1,1,1)$, by writing

$$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},\qquad h_{\mu\nu}\ll 1,$$

(8)

where $h_{\mu\nu}$ represents a small perturbation (dimensionless in geometrized units), and work in the Newtonian limit where gravitational fields are static, therefore time derivatives can be neglected ($\partial_{t}h_{00}=0$) and the d’Alembert operator $\Box$ reduces to the spatial Laplacian $\nabla^{2}$, moreover, in this regime we consider matter sources such that $T_{00}=\rho\gg|T_{ij}|$. Focusing on the time-time component of the perturbation and by using the trace reverse tensor

$$\bar{h}_{\mu\nu}=h_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}h,$$

(9)

the linearized equation takes the form

$$(\alpha^{2}\nabla^{2}+2\gamma)\bar{h}_{00}+(\nabla^{4}+\gamma)\bar{h}=2\left(\beta\rho+\gamma\right),$$

(10)

by neglecting the linear term due to Eq. (8) and reverting to the original perturbation $h_{\mu\nu}$, we arrive at the modified Poisson equation

$$\alpha^{2}\nabla^{2}h_{00}-\nabla^{4}h_{00}=\beta\rho+\gamma,$$

(11)

this is a fourth-order elliptic differential equation for the gravitational potential $h_{00}$ in the Newtonian limit. The term $\alpha^{2}\nabla^{2}h_{00}$ is associated to the standard Newtonian behavior, while the biharmonic operator $\nabla^{4}$ reflects the higher-derivative nature of the underlying $f(R)$ theory and introduces modifications relevant at short distances or high curvature. The source includes both, energy density, $\rho$, and a constant term, $\gamma$, acting as an effective background source related to the cosmological constant in the weak-field regime.

We now proceed to solve the field equation (11) in the presence of a static, spherically symmetric source. In this case, the metric perturbation depends only on the radial coordinate $r$, and we neglect the cosmological constant by setting $\gamma=0$ to isolate local gravitational effects. Therefore, Eq. (11) reduces to

$$\alpha^{2}\left(2h_{00}^{\prime}(r)+rh_{00}^{\prime\prime}(r)\right)-4h_{00}^{(3)}(r)-rh_{00}^{(4)}(r)=\beta r\rho(r),$$

(12)

this fourth-order ODE admits an analytical solution through successive integration, for this we can use the following identity

$$\frac{d}{dr}\left(r^{2}\frac{dh_{00}(r)}{dr}\right)=r\frac{d^{2}}{dr^{2}}\left[rh_{00}(r)\right],$$

(13)

which, when applied recursively, transforms the equation into

$$\frac{d^{2}}{dr^{2}}\left(\alpha^{2}rh_{00}(r)-\frac{d^{2}}{dr^{2}}\left[rh_{00}(r)\right]\right)=\beta r\rho(r),$$

(14)

this form is particularly convenient for integration, allowing the use of definite integrals with carefully chosen limits. The resulting solution is

$$h_{00}(r)=-\frac{\beta e^{\alpha r}}{r}\int_{r}^{\infty}e^{-2\alpha u}\int_{0}^{u}e^{\alpha q}\int_{0}^{q}\int_{p}^{\infty}s\rho(s)\,ds\,dp\,dq\,du+\frac{c_{4}e^{\alpha r}}{r}+\frac{c_{3}e^{-\alpha r}}{r}+\frac{c_{2}}{r}+c_{1},$$

(15)

where constants are subjected to the boundary conditions

$$c_{1}=\left[h_{00}(r)+rh_{00}^{\prime}(r)-\frac{1}{\alpha^{2}}\left(3h_{00}^{\prime\prime}(r)+rh_{00}^{(3)}(r)\right)\right]_{r\to\infty},$$

(16)

$$c_{2}=\left[rh_{00}(r)-\frac{1}{\alpha^{2}}\left(2h_{00}^{\prime}(r)+rh_{00}^{\prime\prime}(r)\right)\right]_{r=0},$$

(17)

$$c_{3}=\frac{c_{1}}{2\alpha^{3}}-\frac{1}{2\alpha}\left[h_{00}(r)+\left(r-\frac{2}{\alpha}\right)h_{00}^{\prime}(r)-\frac{1}{\alpha}rh_{00}^{\prime\prime}(r)\right]_{r=0},$$

(18)

and

$$c_{4}=\left[re^{-\alpha r}h_{00}(r)\right]_{r\to\infty}.$$

(19)

For physically admissible solutions, constant $c_{4}$ vanishes, as it quantifies exponentially growing modes incompatible with asymptotic flatness, this implies that $h_{00}(r)$ must decay as $O(1/r)$ or faster, and, in turn $rh_{00}^{\prime}\sim O(1/r)$, $h_{00}^{\prime\prime}\sim O(1/r^{3})$ and $rh_{00}^{(3)}\sim O(1/r^{3})$, therefore, for localized gravitational fields $c_{1}=0$. At the other hand, $c_{2}$ and $c_{3}$ are determined by the behavior of $h_{00}(r)$ at $r=0$. We are going to consider that the perturbation is analytical near $r=0$, so we can expand it in Taylor series

$$h_{00}(r)\simeq h_{00}(0)+rh_{00}^{\prime}(0)+\frac{1}{2}r^{2}h_{00}^{\prime\prime}(0)+\cdots,$$

(20)

so when replacing in Eq. (17), it is obtained

$$c_{2}=-\frac{2}{\alpha^{2}}h_{00}^{\prime}(0),$$

(21)

however if the potential is symmetric at the origin, $c_{2}$ vanishes. In the same way,

$$c_{3}=\frac{1}{\alpha}h_{00}^{\prime}(0)-\frac{1}{2}h_{00}(0),$$

(22)

thus, if we also assume that the potential vanishes at the origin, i.e., $h_{00}(0)=0$, then $c_{3}=0$. Nevertheless, an interesting observation arises when analyzing the behavior of the modified potential near the origin. By combining expressions for $c_{2}$ and $c_{3}$, we obtain

$$\left[h(r)(\alpha r-1)-rh^{\prime}(r)\right]_{r=0}=0,$$

(23)

thus, solutions with the form $c_{1}e^{\alpha r}/r$ are allowed, even though they are divergent at the origin, while still yielding vanishing integration constants. A series expansion reveals that potentials with the asymptotic form

$$h(r)=\frac{A}{r}+\alpha A+\frac{1}{2}\alpha^{2}Ar+O(r)^{2},$$

(24)

where $A$ is some constant, can satisfy $c_{2}=c_{3}=0$ and at the same time retain the $1/r$ singularity. In summary, when the potential is finite and analytic at the origin and vanishes at infinity, all integration constants vanish ($c_{1}=c_{2}=c_{3}=c_{4}=0$), and the solution is fully determined by the source $\rho(r)$. However, the conditions $c_{2}=c_{3}=0$ do not strictly require regularity at the origin, as they are also satisfied by divergent solutions of the form (24). Under these assumptions, and through successive integration by parts, we arrive at the explicit form of the metric perturbation

$$h_{00}(r)=-\frac{\beta}{\alpha^{2}}\left[\int_{r}^{\infty}\left(1-\frac{e^{\alpha(r-s)}}{2\alpha r}\right)s\rho(s)ds+\frac{1}{r}\int_{0}^{r}\left(s-\frac{e^{\alpha(s-r)}}{2\alpha}\right)s\rho(s)ds\right]+\frac{ce^{-\alpha r}}{\alpha r},$$

(25)

where the constant $c$, given by

$$c=-\frac{\beta}{2\alpha^{2}}\int_{0}^{\infty}s\rho(s)\,ds,$$

(26)

encodes a radial moment of the source distribution, weighted linearly in the radius. The expression for the potential, Eq. (25), consists of two integrals that encode the contributions from the matter distribution both inside and outside a given radius $r$, each modulated by exponential factors arising from the modified gravitational dynamics. The first integral is explicitly nonlocal, as it includes the contribution from the mass located beyond the radius $r$, modulated by an exponential factor that suppresses distant sources. In contrast, the second integral depends solely on the matter distribution inside the radius $r$, and is therefore a local contribution. While the exponential factor $e^{\alpha(s-r)}$ slightly enhances the influence of matter closer to the boundary (i.e., near $s\approx r$), it does not introduce any dependence on regions beyond the point of evaluation. As such, this integral represents a modified yet still local interior contribution to the gravitational potential. In the limit $\alpha\to\infty$, the exponential corrections vanish, and the integrals reduce to the standard expressions for the interior and exterior mass contributions in Newtonian gravity. Thus, potential (25) reflects the interplay between the radial structure of the source and the finite-range nature of gravity induced by the modified theory.

In this section, we evaluate the modified gravitational potential for a variety of spherically symmetric mass distributions. These include idealized models such as the spherical shell and homogeneous sphere, astrophysically motivated profiles like Plummer, Hernquist, and NFW, as well as new density models introduced for comparison. Each profile is examined for its qualitative behavior and its implications under the modified gravity framework. For all cases we set $\beta=-8\pi\alpha^{2}$, and it is found that integration constants vanish ($c_{1}=c_{2}=c_{3}=c_{4}=0$).

In this work, the gravitational potential was analyzed within the framework of $f(R)$ gravity, considering a quadratic form of the function $f(R)$. The resulting potential satisfies a fourth-order differential equation (11), which was solved analytically under the assumption of spherical symmetry. This equation exhibits two key features: (i) For $\alpha^{2}\to\infty$ and finite $\gamma$, the $\nabla^{4}$ term is suppressed, reducing to the standard Poisson equation with cosmological constant; and (ii) Yukawa-type corrections from the $\nabla^{4}$ operator introduces massive modes generating exponentially screened potentials.
By imposing physical boundary conditions - regularity at the origin and asymptotic flatness - all integration constants were shown to vanish. This guarantees that the potential, Eq. (25), is uniquely determined by the mass distribution that sources the field. Moreover, interestingly, it was demonstrated that even potentials exhibiting a mild divergence of the form $h(r)\sim A/r+\alpha A+\cdots$, near the origin, can still satisfy these boundary conditions, and thus lead to vanishing integration constants.
The solution for the modified gravitational potential, Eq. (25), was found to consist of the standard GR term supplemented by Yukawa-type corrections of the form $e^{-\alpha r}/r$, this term encodes the additional scalar interaction mediated by the scalaron field, with $\alpha^{-1}$ setting the screening scale. Accordingly, in the limit $\alpha\to\infty$, the potential smoothly reduces to the Newtonian form, fully recovering GR. Similarly, at large distances $r\gg\alpha^{-1}$, the potential exhibits classical behavior.
The potential was examined for a range of spherically symmetric mass distributions, including idealized configurations, astrophysically motivated profiles, and proposed novel density functions. In all cases, the solutions were found to be continuous and smooth throughout the entire domain. The boundary conditions consistently led to vanishing integration constants, confirming that the potential is entirely determined by the matter distribution. While the Newtonian limit is recovered as $\alpha\to\infty$, finite values of $\alpha$ introduce nontrivial corrections at intermediate scales $r\sim\alpha^{-1}$, which could have important phenomenological implications in gravitational systems.
Our analysis reveals a fundamental characteristic of the modified potentials in $f(R)$ gravity: while all solutions exhibit a $1/r$-type divergence at the origin ($r\to 0$) for finite $\alpha$, this singularity is naturally regularized in the GR limit ($\alpha\to\infty$), where the potentials approach finite, constant values.
In the case of diffuse and regular profiles like the Plummer model, the corrections are more pronounced at intermediate radii. In contrast, for cuspy profiles such as Hernquist and NFW, which exhibit a $1/r$ divergence in density near the origin, the modified potential reflects the central steepness, leading to sharper deviations at small radii. Meanwhile, for compact configurations like the spherical shell or exponentially decaying distributions, the corrections are localized near the origin or the outer boundary.
Furthermore, we introduced new density profiles motivated by analytical simplicity and mathematical tractability. These models allowed us to observe how different radial decay and central slopes influence the gravitational potential. For the exponential cutoff model, Eq. (44), the overall shape of the potential, Fig. 5, closely resembles that of the homogeneous sphere, Fig. 1. In contrast, the exponential singular model, Eq. (55), like Hernquist and NFW profiles, diverges at the origin and decays rapidly yet smoothly to zero at large radii. Consequently, its gravitational potential exhibits similar features: a divergence at $r=0$ and a fast, smooth convergence toward the classical Newtonian form.

The novel density profiles introduced in this work, characterized by simple analytical forms and tunable parameters such as decay rates or core slopes, offer alternatives for modeling astrophysical systems.

Although this work has established the mathematical framework for the gravitational potential solutions, we emphasize that our focus has been on their theoretical properties rather than observational detection. Future research could extend this study by exploring numerical simulations and strong gravity regimes (considered the most promising laboratories) where shallow potentials enhance the effects of modified gravity.