Constraining Quadratic $f(R)$ Gravity from Astrophysical Observations of the Pulsar J0704+6620

We apply quadratic $f(R)=R+\epsilon R^2$ field equations, where $\epsilon$ has a dimension [L$^2$], to static spherical stellar model. We assume the interior configuration is determined by Krori-Barua ansatz and additionally the fluid is anisotropic. Using the astrophysical measurements of the pulsar PSR J0740+6620 as inferred by NICER and XMM observations, we determine $\epsilon\approx \pm 3$ km$^2$. We show that the model can provide a stable configuration of the pulsar PSR J0740+6620 in both geometrical and physical sectors. We show that the Krori-Barua ansatz within $f(R)$ quadratic gravity provides semi-analytical relations between radial, $p_r$, and tangential, $p_t$, pressures and density $\rho$ which can be expressed as $p_r\approx v_r^2 (\rho-\rho_1)$ and $p_r\approx v_t^2 (\rho-\rho_2)$, where $v_r$ ($v_t$) is the sound speed in radial (tangential) direction, $\rho_1=\rho_s$ (surface density) and $\rho_2$ are completely determined in terms of the model parameters. These relations are in agreement with the best-fit equations of state as obtained in the present study. We further put the upper limit on the compactness, which satisfies the $f(R)$ modified Buchdahl limit. Interestingly, the quadratic $f(R)$ gravity with negative $\epsilon$ naturally restricts the maximum compactness to values lower than Buchdahl limit, unlike the GR or $f(R)$ gravity with positive $\epsilon$ where the compactness can arbitrarily approach the black hole limit $C\to 1$. The model predicts a core density a few times the saturation nuclear density $\rho_{\text{nuc}} = 2.7\times 10^{14}$ g/cm$^3$, and a surface density $\rho_s > \rho_{\text{nuc}}$. We provide the mass-radius diagram corresponding to the obtained boundary density which has been shown to be in agreement with other observations.