  When a rotating neutron star loses angular momentum, the progressive
reduction of the centrifugal force makes it contract. This perturbs each fluid
element, raising the local pressure and originating deviations from beta
equilibrium, inducing reactions that release heat (rotochemical heating). This
effect has previously been studied by Fern\'andez & Reisenegger (2005) for
non-superfluid neutron stars and by Petrovich & Reisenegger (2010) for
superfluid millisecond pulsars. Both studies found that pulsars reach a
quasi-steady state in which the compression driving the matter out of beta
equilibrium is balanced by the reactions trying to restore the equilibrium. We
extend previous studies by considering the effect of density-dependence and
anisotropy of the superfluid energy gaps, for the case in which the dominant
reactions are the modified Urca processes, the protons are non-superconducting,
and the neutron superfluidity is parametrized by models proposed in the
literature. By comparing our predictions with the surface temperature of the
millisecond pulsar PSR J0437-4715 and upper limits for twenty-one classical
pulsars, we find the millisecond pulsar can be only explained by the models
with the effectively largest energy gaps (type B models), the classical pulsars
require with the gap models that vanish for some angle (type C) and two
different envelope compositions. Thus, no single model for neutron
superfluidity can simultaneously account for the thermal emission of all
available observations of non-accreting neutron stars, possibly due to our
neglect of proton superconductivity.
