  In this work we consider an axionic scalar-tensor theory of gravity and its
effects on static neutron stars. The axionic theory is considered in the regime
in which the axion oscillates around its potential minimum, which
cosmologically occurs post-inflationary, when the Hubble rate is of the same
order as the axion mass. We construct the TOV equations for this axionic theory
and for a spherically symmetric static spacetime and we solve these numerically
using a quite robust double shooting LSODA based python integration method.
Regarding the equations of state, we used nine mainstream and quite popular
ones, namely, the WFF1, the SLy, the APR, the MS1, the AP3, the AP4, the ENG,
the MPA1 and the MS1b, using the piecewise polytropic description for each.
From the extracted data we calculate the Jordan frame masses and radii, and we
confront the resulting phenomenology with five well-known neutron star
constraints. As we demonstrate, the AP3, the ENG and the MPA1 equations of
state yield phenomenologically viable results which are compatible with the
constraints, with the MPA1 equation of state enjoying an elevated role among
the three. The reason is that the MPA1 fits well the phenomenological
constraints. A mentionable feature is the fact that all the viable
phenomenologically equations of state produce maximum masses which are in the
mass-gap region with $M_{max}>2.5M_{\odot}$, but lower that the causal 3 solar
masses limit. We also compare the neutron star phenomenology produced by the
axionic scalar-tensor theory with the phenomenology produced by inflationary
attractors scalar-tensor theories.
