We apply quadratic f(R) = R + ϵR
2 field equations, where ϵ has a dimension [L2], 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 ϵ ≈ ± 3 km2. 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, pr
, and tangential, pt
, pressures and density ρ which can be expressed as pr
≈ vr
2 (ρ-ρ
1) and pr
≈ vt
2 (ρ-ρ
2), where vr
(vt
) is the sound speed in radial (tangential) direction, ρ
1 = ρs
(surface density) and ρ
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, C = 2GMRs
-1
c
-2, which satisfies the f(R) modified Buchdahl limit. Remarkably, the quadratic f(R) gravity with negative ϵ naturally restricts the maximum compactness to values lower than Buchdahl limit, unlike the GR or f(R) gravity with positive ϵ where the compactness can arbitrarily approach the black hole limit C → 1. The model predicts a core density a few times the saturation nuclear density ρ
nuc = 2.7 × 1014 g/cm3, and a surface density ρs
> ρnuc
. We provide the mass-radius diagram corresponding to the obtained boundary density which has been shown to be in agreement with other observations.