Quiescent hard X-ray and soft gamma-ray emission from
neutron stars constitute a promising frontier to explore
axion-like-particles (ALPs). ALP production in the core peaks at
energies of a few keV to a few hundreds of keV; subsequently, the
ALPs escape and convert to photons in the magnetosphere. The
emissivity goes as ∼ T
6 while the conversion probability is
enhanced for large magnetic fields, making magnetars, with their
high core temperatures and strong magnetic fields, ideal targets for
probing ALPs. We compute the energy spectrum of photons resulting
from conversion of ALPs in the magnetosphere and then compare it
against hard X-ray data from NuSTAR, INTEGRAL, and XMM-Newton for a set
of eight magnetars for which such data exists. Upper limits are
placed on the product of the ALP-nucleon and ALP-photon
couplings. For the production in the core, we perform a calculation
of the ALP emissivity in degenerate nuclear matter modeled by a
relativistic mean field theory. The reduction of the emissivity due
to improvements to the one-pion exchange approximation is
incorporated, as is the suppression of the emissivity due to proton
superfluidity in the neutron star core. A range of core temperatures
is considered, corresponding to different models of the steady heat
transfer from the core to the stellar surface. For the subsequent
conversion, we solve the coupled differential equations mixing ALPs
and photons in the magnetosphere. The conversion occurs due to a
competition between the dipolar magnetic field and the photon
refractive index induced by the external magnetic
field. Semi-analytic expressions are provided alongside the full
numerical results. We also present an analysis of
the uncertainty on the axion limits we derive due to the
uncertainties in the magnetar masses, nuclear matter equation of
state, and the proton superfluid critical temperature.