  We investigated realistic neutron stars in axion $R^{2}$ gravity. The
coupling between curvature and axion field $\phi$ is assumed in the simple form
$\sim R^2\phi$. For the axion mass in the range $m_{a}\sim 10^{-11}-10^{-10}$
eV the solitonic core within neutron star and corresponding halo with size
$\sim 100$ km can exist. Therefore the effective contribution of $R^2$ term
grows inside the star and it leads to change of star parameters (namely, mass
and radius). We obtained the increase of star mass independent from central
density for wide range of masses. Therefore, maximal possible mass for given
equation of state grows. At the same time, the star radius increases not so
considerably in comparison with GR. {Hence, our model may predict possible
existence of supermassive compact stars with masses $M\sim 2.2-2.3M_\odot$ and
radii $R_{s}\sim 11$ km for realistic equation of state (we considered APR
equation of state). In General Relativity one can obtain neutron stars with
such characteristics only for unrealistic, extremely stiff equations of state.}
Note that this increase of mass occurs due to change of solution for scalar
curvature outside the star. In GR curvature drops to zero on star surface where
$\rho=p=0$. In the model under consideration the scalar curvature dumps more
slowly in comparison with vacuum $R^2$ gravity due to axion "galo" around the
star.
