  Photon-axion mixing can create observable signatures in thermal spectra of
isolated, cooling neutron stars. Their shape depends on the polarization
properties of the radiation, which, in turn, are determined by the structure of
stellar outermost layers. Here we investigate the effect of mixing on the
spectrum and polarimetric observables, polarization fraction and polarization
angle, using realistic models of surface emission. We focus on RX J1856.5-3754,
the only source among the X-ray dim isolated neutron stars for which
polarimetric measurements in the optical band were performed. Our results show
that in the case of a condensed surface in both fixed and free-ion limits, the
mixing can significantly limit the geometric configurations which reproduce the
observed linear polarization fraction of 16.43%. In the case of an atmosphere,
the mixing does not create any noticeable signatures. Complementing our
approach with the data from upcoming soft X-ray polarimetry missions will allow
to obtain constraints on $g_{\gamma a} \sim 10^{-11}$ GeV$^{-1}$ and $m_a
\lesssim 10^{-6}$ eV, improving the present experimental and astrophysical
limits.
