  X-ray emission from the surface of isolated neutron stars (NSs) has been now
observed in a variety of sources. The ubiquitous presence of pulsations clearly
indicates that thermal photons either come from a limited area, possibly heated
by some external mechanism, or from the entire (cooling) surface but with an
inhomogeneous temperature distribution. In a NS the thermal map is shaped by
the magnetic field topology, since heat flows in the crust mostly along the
magnetic field lines. Self-consistent surface thermal maps can hence be
produced by simulating the coupled magnetic and thermal evolution of the star.
We compute the evolution of the neutron star crust in three dimensions for
different initial configurations of the magnetic field and use the ensuing
thermal surface maps to derive the spectrum and the pulse profile as seen by an
observer at infinity, accounting for general-relativistic effects. In
particular, we compare cases with a high degree of symmetry with inherently 3D
ones, obtained by adding a quadrupole to the initial dipolar field. Axially
symmetric fields result in rather small pulsed fractions ($\lesssim 5\%$),
while more complex configurations produce higher pulsed fractions, up to
$\sim25\%$. We find that the spectral properties of our axisymmetric model are
close to those of the bright isolated NS RX~J1856.5-3754 at an evolutionary
time comparable with the inferred dynamical age of the source.
