  Differential rotation between the neutron star crust and a more rapidly
rotating interior superfluid leads to frictional heating that affects the
star's long-term thermal evolution and resulting surface emission. Here we
present the results of thermal evolution simulations based on two models of the
vortex pinning forces that bracket a range of plausible pinning strengths. We
include the effects of superfluidity, magnetic fields, and temperature
gradients. As representative standard and accelerated neutrino emission
processes taking place in the core, we consider the modified Urca process in
normal baryonic matter, and the much faster quark Urca process. Comparison of
our results with neutron star surface temperature data, including the recent
temperature measurement of the Geminga pulsar, shows that stars with soft
equations of state and modest frictional heating are in closest agreement with
the data; stars with stronger frictional heating have temperatures inconsistent
with the upper limit of PSR 1929+10. Stiffer stars undergoing standard cooling
generally have temperatures lying above the Vela detection, a situation
worsened by the inclusion of frictional heating. Stars undergoing accelerated
cooling without frictional heating have temperatures that fall far below most
temperature measurements; the Vela and Geminga detections being the most
compelling examples. Only in stiff stars, which have thick crusts, can the
inclusion of strong frictional heating raise the temperature at late stages in
the evolution to a level consistent with the data. However, such a large amount
of heating leads to a temperature at $\sim 1000$ yr in excess of the Crab upper
limit. Suppression of accelerated neutrino emission processes, perhaps by
superfluid pairing in the core, may yield acceptable cooling models.
