Conversion of photons into axions under the presence of a strong magnetic field can dim the radiation from magnetized astrophysical objects. Here we perform a detailed calculation aimed at quantifying the signatures of photon-axion conversion in the spectra, light curves, and polarization of neutron stars (NSs). We take into account the energy and angle dependence of the conversion probability and the surface thermal emission from NSs. The latter is computed from magnetized atmosphere models that include the effect of photon polarization mode conversion due to vacuum polarization. The resulting spectral models, inclusive of the general-relativistic effects of gravitational redshift and light deflection, allow us to make realistic predictions for the effects of photon to axion conversion on observed NS spectra, light curves, and polarization signals. We identify unique signatures of the conversion, such as an increase of the effective area of a hot spot as it rotates away from the observer line of sight. For a star emitting from the entire surface, the conversion produces apparent radii that are either larger or smaller (depending on axion mass and coupling strength) than the limits set by NS equations of state. For an emission region that is observed phase-on, photon-axion conversion results in an inversion of the plane of polarization with respect to the no-conversion case. While the quantitative details of the features that we identify depend on NS properties (magnetic field strength and temperature) and axion parameters, the spectral and polarization signatures induced by photon-axion conversion are distinctive enough to make NSs very interesting and promising probes of axion physics.