  Neutron stars harbour extremely strong magnetic fields within their solid
outer crust. The topology of this field strongly influences the surface
temperature distribution, and hence the star's observational properties. In
this work, we present the first realistic simulations of the coupled crustal
magneto-thermal evolution of isolated neutron stars in three dimensions with
account for neutrino emission, obtained with the pseudo-spectral code Parody.
We investigate both the secular evolution, especially in connection with the
onset of instabilities during the Hall phase, and the short-term evolution
following episodes of localised energy injection. Simulations show that a
resistive tearing instability develops in about a Hall time if the initial
toroidal field exceeds ~$10^{15}$ G. This leads to crustal failures because of
the huge magnetic stresses coupled with the local temperature enhancement
produced by dissipation. Localised heat deposition in the crust results in the
appearance of hot spots on the star surface which can exhibit a variety of
patterns. Since the transport properties are strongly influenced by the
magnetic field, the hot regions tend to drift away and get deformed following
the magnetic field lines while cooling. The shapes obtained with our
simulations are reminiscent of those recently derived from NICER X-ray
observations of the millisecond pulsar PSR J0030+0451.
