Abstract We present a complete analytic and semi-analytic study of gravitational collapse and primordial black hole (PBH) formation in the quadratic f (R) model f (R) =R+ R² f (R) = R + α R 2. We first derive the perturbative expansion around General Relativity (GR), working to first order in the small parameter α. For a collapsing flat FLRW dust interior we compute the explicit first-order corrections to the scale factor, the stellar radius, and the horizon formation time. We then use these results to infer the expected trend in the PBH formation threshold c δ c, rather than a direct quantitative determination. Within this perturbative framework, the quadratic correction modifies the dust collapse dynamics at first order, while the flat radiation-dominated background does not receive a nontrivial correction at the same order. As result, any modification of PBH formation in the radiation era must arise from nonlinear or non-perturbative effects. The perturbative analysis therefore provides qualitative insight into how curvature corrections influence collapse, particularly in high-curvature regimes. To access this regime we reformulate the theory in the Einstein frame, where the model becomes GR plus the scalaron field ϕ with the Starobinsky potential. We provide the complete ODE system governing both the cosmological background and the evolution of an overdense closed FLRW patch. This system can be numerically integrated to obtain the critical overdensity c (k) δ c (k) for PBH formation near the end of inflation.
Nashed et al. (Sun,) studied this question.