  There are many factors that contribute to the breaking of the spherical
symmetry of a neutron star. Most notably is rotation, magnetic fields, and/or
accretion of matter from companion stars. All these phenomena influence the
macroscopic structures of neutron stars, but also impact their microscopic
compositions. The purpose of this paper is to investigate the cooling of
rotationally deformed, two-dimensional (2D) neutron stars in the framework of
general relativity theory, with the ultimate goal of better understand the
impact of 2D effects on the thermal evolution of such objects. The equations
that govern the thermal evolution of rotating neutron stars are presented in
this paper. The cooling of neutron stars with different frequencies is computed
self-consistently by combining a fully general relativistic 2D rotation code
with a general relativistic 2D cooling code. We show that rotation can
significantly influence the thermal evolution of rotating neutron stars. Among
the major new aspects are the appearances of hot spots on the poles, and an
increase of the thermal coupling times between the core and the crust of
rotating neutron stars. We show that this increase is independent of the
microscopic properties of the stellar core, but depends only on the frequency
of the star.
