  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 $\sim 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 careful 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 strong suppression of the emissivity
due to nucleon 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. Our treatment also includes
the first calculation of the emissivity due to $n+p \rightarrow n+p+a$
processes in the limit of strongly degenerate nuclear matter. 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.
