Tidal Love numbers of neutron stars in f(R) gravity

The recent detection of gravitational waves from a neutron star merger was a significant step towards constraining the nuclear matter equation of state by using the tidal Love numbers (TLNs) of the merging neutron stars. Measuring or constraining the neutron star TLNs allows us in principle to exclude or constraint many equations of state. This approach, however, has the drawback that many modified theories of gravity could produce deviations from General Relativity similar to the deviations coming from the uncertainties in the equation of state. The first and the most natural step in resolving the mentioned problem is to quantify the effects on the TLNs from the modifications of General Relativity. With this motivation in mind, in the present paper we calculate the TLNs of (non-rotating) neutron stars in $f(R)$ gravity. For this purpose, we first derived the equations describing both the polar and the axial stationary perturbations of neutron stars in a particular class of $f(R)$ gravity, the so-called $R^2$-gravity. Then, by solving numerically the perturbation equations, we calculate explicitly the polar and the axial $l=2$ TLNs of the neutron stars in $R^2$-gravity for three characteristic realistic equations of state. Our results show that while the polar TLNs are slightly influenced by the $R^2$ modification of General Relativity, the axial TLNs can be several times larger (in terms of the absolute value) compared to the general relativistic case.