  Recent observations of X-ray pulsars at low luminosities allow, for the first
time, to compare theoretical models for the emission from highly magnetized
neutron star atmospheres at low mass accretion rates ($\dot{M} \lesssim
10^{15}$ g s$^{-1}$) with the broadband X-ray data. The purpose of this paper
is to investigate the spectral formation in the neutron star atmosphere at low
$\dot{M}$ and to conduct a parameter study of physical properties of the
emitting region. We obtain the structure of the static atmosphere, assuming
that Coulomb collisions are the dominant deceleration process. The upper part
of the atmosphere is strongly heated by the braking plasma, reaching
temperatures of 30-40 keV, while its denser isothermal interior is much cooler
(~2 keV). We numerically solve the polarized radiative transfer in the
atmosphere with magnetic Compton scattering, free-free processes, and
non-thermal cyclotron emission due to possible collisional excitations of
electrons. The strongly polarized emitted spectrum has a double-hump shape that
is observed in low-luminosity X-ray pulsars. A low-energy "thermal" component
is dominated by extraordinary photons that can leave the atmosphere from deeper
layers due to their long mean free path at soft energies. We find that a
high-energy component is formed due to resonant Comptonization in the heated
non-isothermal part of the atmosphere even in the absence of collisional
excitations. The latter, however, affect the ratio of the two components. A
strong cyclotron line originates from the optically thin, uppermost zone. A fit
of the model to NuSTAR and Swift/XRT observations of GX 304-1 provides an
accurate description of the data with reasonable parameters. The model can thus
reproduce the characteristic double-hump spectrum observed in low-luminosity
X-ray pulsars and provides insights into spectral formation.
