  In previous papers, we have shown that, as the rotation of a neutron star
slows down, it will be internally heated as a consequence of the progressively
changing mix of particles (rotochemical heating). In previously studied cases
(non-superfluid neutron stars or superfluid stars with only modified Urca
reactions), this leads to a quasi-steady state in which the star radiates
thermal photons for a long time, possibly accounting for the ultraviolet
radiation observed from the millisecond pulsar J0437-4715. For the first time,
we explore the phenomenology of rotochemical heating with direct Urca reactions
and uniform and isotropic superfluid energy gaps of different sizes. We first
do exploratory work by integrating the thermal and chemical evolution equations
numerically for different energy gaps, which uncovers a rich phenomenology of
stable and unstable solutions. To understand these, we perform a stability
analysis around the quasi-steady state, identifying the characteristic times of
growing, decaying, and oscillating solutions. For small gaps, the phenomenology
is similar to the previously studied cases, in the sense that the solutions
quickly converge to a quasi-steady state. For large gaps ($\gtrsim 0.05$ MeV),
these solutions become unstable, leading to a limit-cycle behavior of
periodicity $\sim 10^{6-7} yr, in which the star is hot ($T_s\gtrsim 10^5$ K)
for a small fraction of the cycle ($\sim 5- 20 %$), and cold for a longer time.
