  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.
