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.3 M _ { \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 .