Supermassive Neutron Stars in Starobinsky Gravity with Causal Hybrid Stellar Matter

We investigate the stellar structure of neutron stars in the framework of Starobinsky gravity, characterized by a quadratic correction to the Einstein-Hilbert action, $f(R) = R + αR^2$. In order to {\em preserve causality throughout\/} the star, we adopt a two-phase hybrid construction for the stellar matter, in which the core region consists of deconfined quark matter described by the MIT bag model, while the outer layers are composed of hadronic matter represented by unified equations of state such as SLy4, BSk20, and BSk21. Within this framework, we derive the modified field equations with static spherical symmetry, and numerically integrate the corresponding Tolman-Oppenheimer-Volkoff (TOV) equations with the chosen hybrid equations of state. Our analyses show that, unlike in general relativity, the Ricci scalar remains nonzero outside the stellar surface, and gradually falls to zero beyond 50 km, while the stellar surface remains within 10--12 km for the hybrid equations of state considered. This extended Ricci scalar profile arises from the extra degree of freedom (scalaron) inherent in Starobinsky gravity, which also contributes to the gravitational mass outside the star, causing the ADM mass measured at infinity to exceed the stellar mass at the surface. Nevertheless, the key physical relationships, such as the mass-central density and mass-radius curves, remain consistent with what is expected physically. Notably, we find that the maximum stable mass of neutron stars increases with the Starobinsky parameter $α$, with the combined MIT-BSk21 model supporting an ADM mass of up to $2.07 \, M_\odot$ for $α= 10\,r_g^2$. This theoreical limit for a nonrotating neutron star suggests that a rotating configuration could reach mass thresholds in the range of $2.48$ to $2.59 \, M_\odot$, considering that rapid rotation can enhance the maximum mass by approximately 20--25\%.