  Axions can be copiously produced in localized regions of neutron star
magnetospheres where the ambient plasma is unable to efficiently screen the
induced electric field. As these axions stream away from the neutron star they
can resonantly transition into photons, generating a large broadband
contribution to the neutron star's intrinsic radio flux. In this work, we
develop a comprehensive end-to-end framework to model this process from the
initial production of axions to the final detection of radio photons, and
derive constraints on the axion-photon coupling, $g_{a\gamma\gamma}$, using
observations of 27 nearby pulsars. We study the modeling uncertainty in the
sourced axion spectrum by comparing predictions from 2.5 dimensional
particle-in-cell simulations with those derived using a semi-analytic model;
these results show remarkable agreement, leading to constraints on the
axion-photon coupling that typically differ by a factor of no more than $\sim
2$. The limits presented here are the strongest to date for axion masses
$10^{-8} \, {\rm eV} \lesssim m_a \lesssim 10^{-5} \, {\rm eV}$, and crucially
do not rely on the assumption that axions are dark matter.
