Cooling rates are first calculated for neutron stars of about 1 M
⊙ and 10 km radius, with magnetic fields from zero to about 1014 G, for two extreme cases of maximum and no superfluidity. The results show that most pulsars are so cold that thermal ionization of surface atoms would be negligible. Next, nucleon superfluidity and crystallization of heavy nuclei are treated more quantitatively, and more realistic hadron star models are chosen. Cooling rates are thus calculated for a stable hyperon star near the maximum mass limit, a medium weight neutron star, and a light neutron star with neutron-rich heavy nuclei near the minimum mass limit. Results show that cooling rates are a sensitive function of density. The lightest star is cooler than others in earlier stages but the trend is reversed later. The medium weight star is generally the coldest of all in lower temperature regions where the effect of superfluidity becomes significant. However, if a heavy star contains pions, its cooling will be even faster. The Crab pulsar and Vela pulsar, expected to be the two youngest, can be as hot as (2 ∼ 4) x 106 K (on the surface), comparable with the results obtained from internal frictional heating by Greenstein, if they are medium weight to heavy hadron stars. However, older pulsars are cold. In fact, at about a few million years, the age of average radio pulsars, the surface temperature becomes only several hundred to several thousand degrees. Thus, the earlier conclusions about cold pulsars are still valid. Cooling of a massive white dwarf star is also shown.