and the “apparent” effective temperature $T_{\rm s}^\infty = T_{\rm s} \, \sqrt{1 - r_{\rm g}/R}$, where $r_{\rm g}=2GM/c^2$ is the Schwarzschild radius and $M$ is the gravitational neutron star mass. If the surface temperature distribution is anisotropic (for instance, owing to a strong magnetic field), $L_\gamma^\infty$ is determined by a properly averaged surface temperature (e.g., Ref. [@py01]). The theory gives cooling curves, $T_s^\infty(t)$, and predicts three main cooling stages. The initial stage of internal thermal relaxation lasts $t \la 10-50$ yr. At this stage, the neutron star crust stays hotter than the core and thermally decoupled from it (because of much stronger neutrino cooling in the core). The surface temperature is then insensitive to the physics of the core but strongly depends on physical properties of the crust [@lattimeretal94; @gyp01]. Since no neutron star has been observed at this stage, we will not discuss it in detail. The next stage of neutrino cooling with isothermal interior lasts $t \la 10^5-10^6$ yr. At this stage the neutrino luminosity $L_\nu \gg L_\gamma$; thermal conduction is high and makes the stellar interior isothermal, with the main temperature gradient located in the heat blanketing envelope. The cooling is mostly regulated by a strong neutrino 