  Theoretical modeling of surface emission from magnetized neutron stars (NSs)
requires proper treatment of QED effects, in particular the effect of photon
mode conversion due to the ``vacuum resonance'' between plasma and vacuum
polarization. Previous NS atmosphere models incorporated this effect
approximately, using transfer equations for the photon modes which are
inadequate near the vacuum resonance, particularly for field strengths around
$B_l\simeq 7\times 10^{13}$ G, where the vacuum resonance occurs near the
photosphere. In this paper, we provide an accurate treatment of the QED-induced
mode conversion effect in NS atmosphere models, employing both the modal
radiative transfer equations, coupled with an accurate mode conversion
probability at the resonance, and the full evolution equations for the photon
Stokes parameters. In doing so, we are able to quantitatively calculate the
effects of vacuum polarization on the atmosphere emission spectra, beam
patterns, and polarizations for the entire range of field strengths,
$B=10^{12}-10^{15}$ G. We find that for NSs with $B\ga 2 B_l$, vacuum
polarization reduces the widths of spectral features, and softens the hard
spectral tail typical of magnetized atmosphere models. For $B\la B_l/2$, vacuum
polarization does not change the emission spectra, but can significantly affect
the polarization signals. We show that vacuum polarization induces a unique
energy-dependent linear polarization signature, and that circular polarization
can be generated in the magnetospheres of rapidly rotating NSs. We discuss the
implications of our results for recent observations of thermally emitting
isolated NSs and magnetars, as well as the prospects for future spectral and
polarization observations.
