  We study the thermal structure of neutron stars with magnetized envelopes
composed of accreted material, using updated thermal conductivities of plasmas
in quantizing magnetic fields, as well as equation of state and radiative
opacities for partially ionized hydrogen in strong magnetic fields. The
relation between the internal and local surface temperatures is calculated and
fitted by an analytic function of the internal temperature, magnetic field
strength, angle between the field lines and the normal to the surface, surface
gravity, and the mass of the accreted material. The luminosity of a neutron
star with a dipole magnetic field is calculated for various values of the
accreted mass, internal temperature, and magnetic field strength. Using these
results, we simulate cooling of superfluid neutron stars with magnetized
accreted envelopes. We consider slow and fast cooling regimes, paying special
attention to very slow cooling of low-mass superfluid neutron stars. In the
latter case, the cooling is strongly affected by the combined effect of
magnetized accreted envelopes and neutron superfluidity in the stellar crust.
Our results are important for interpretation of observations of isolated
neutron stars hottest for their age, such as RX J0822-43 and PSR B1055-52.
