  Cooling simulations of neutron stars and their comparison with the data from
thermally emitting x-ray sources put constraints on the properties of axions,
and by extension of any light pseudoscalar dark matter particles, whose
existence has been postulated to solve the strong-CP problem of QCD. We
incorporate the axion emission by pair-breaking and formation processes by $S$-
and $P$-wave nucleonic condensates in a benchmark code for cooling simulations
as well as provide fit formulas for the rates of these processes. Axion cooling
of neutron stars has been simulated for 24 models covering the mass range 1 to
1.8 solar masses, featuring nonaccreted iron and accreted light-element
envelopes, and a range of nucleon-axion couplings. The models are based on an
equation state predicting conservative physics of superdense nuclear matter
that does not allow for the onset of fast cooling processes induced by phase
transitions to non-nucleonic forms of matter or high proton concentration. The
cooling tracks in the temperature vs age plane were confronted with the
(time-averaged) measured surface temperature of the central compact object in
the Cas A supernova remnant as well as surface temperatures of three nearby
middle-aged thermally emitting pulsars.We find that the axion coupling is
limited to $f_a/10^{7}\textrm{GeV} \ge (5$--$10)$, which translates into an
upper bound on axion mass $m_a \le (0.06$--$0.12)~\textrm{eV}$ for Peccei-Quinn
charges of the neutron $\vert C_n \vert \sim 0.04$ and proton $\vert C_p \vert
\sim 0.4$ characteristic for hadronic models of axions.
