  The thermal evolution of neutron stars (NSs) is investigated by coupling with
the evolution of $\textit{r}$-mode instability that is described by a second
order model.The heating effect due to shear viscous damping of the
$\textit{r}$-modes enables us to understand the high temperature of two young
pulsars (i.e., PSR B0531+21 and RX J0822-4300) in the framework of the simple
$npe$ NS model, without superfluidity or exotic particles.Moreover, the light
curves predicted by the model within an acceptable parameter regime may
probably cover all of the young and middle-aged pulsars in the $\lg
T_s^{\infty}-\lg t$ panel, and an artificially strong $p$ superfluidity invoked
in some early works is not needed here. Additionally, by considering the
radiative viscous damping of the $\textit{r}$-modes, a surprising extra cooling
effect is found, which can even exceed the heating effect sometimes although
plays an ignorable role in the thermal history.
