  We computed accurate atmosphere models of rotation-powered millisecond
pulsars in which the polar caps of a neutron star (NS) are externally heated by
magnetospheric return currents. The external ram pressure, energy losses, and
stopping depth of the penetrating charged particles were computed
self-consistently with the atmosphere model, instead of assuming a simplified
deep-heated atmosphere in radiative equilibrium. We used exact Compton
scattering formalism to model the properties of the emergent X-ray radiation.
The deep-heating approximation was found to be valid only if most of the heat
originates from ultra-relativistic bombarding particles with Lorentz factors of
$\gamma \gtrsim 100$. In the opposite regime, the atmosphere attains a distinct
two-layer structure with an overheated optically thin skin on top of an
optically thick cool plasma. The overheated skin strongly modifies the emergent
radiation: it produces a Compton-upscattered high-energy tail in the spectrum
and alters the radiation beaming pattern from limb darkening to limb
brightening for emitted hard X-rays. This kind of drastic change in the
emission properties can have a significant impact on the inferred NS pulse
profile parameters as performed, for example, by Neutron star Interior
Composition ExploreR. Finally, the connection between the energy distribution
of the return current particles and the atmosphere emission properties offers a
new tool to probe the exact physics of pulsar magnetospheres.
