  We present a theoretical model for detecting axions from neutron stars in a
QCD phase of quark matter, in which the gamma-rays are produced from decaying
axions. The axions would be produced from a quark-antiquark pair $u\bar{u}$ or
$d\bar{d}$. The chiral anomaly of QCD and the spontaneously broken symmetry are
invoked to explain the non-conservation of the axion current. The axion
quark-antiquark coupling can be computed via the Goldberger-Treiman relation,
relying on the two-loop coupling form factors from the QCD axion current. From
the coupling form factors, the axion emissivities $\epsilon_a$ can be derived,
from which fluxes can be determined. A Monte Carlo simulation is used to
compute the emissivity for axions from said neutron stars, from which the
predicted gamma-ray signal from the Fermi Large Area Telescope (Fermi LAT) can
be computed. We predict a photon flux which may be detectable by Fermi LAT, and
limits on the QCD mass $m_a$. A signal from LIGO GWS 170817 could be produced
from the neutron star merger as well, which is compared to the expected signal
from the Fermi LAT. An upper bound $m_a\lesssim 10^{-10}$eV is thereby placed.
Very stringent limits on the axion mass of nearly $m_a \simeq 10^{-14}$eV could
be placed using this model with the closest neutron stars near 100 kpc. Axions
could thus be detectable in gamma-rays for neutron stars as distant as 1 Mpc.
