  Thermal evolution of neutron stars is known to depend on the properties of
superdense matter in neutron star cores. We suggest a statistical analysis of
isolated cooling middle-aged neutron stars and old transiently accreting
quasi-stationary neutron stars warmed up by deep crustal heating in low-mass
X-ray binaries. The method is based on simulations of the evolution of stars of
different masses and on averaging the results over respective mass
distributions. This gives theoretical distributions of isolated neutron stars
in the surface temperature--age plane and of accreting stars in the photon
thermal luminosity--mean mass accretion rate plane to be compared with
observations. This approach permits to explore not only superdense matter but
also the mass distributions of isolated and accreting neutron stars. We show
that the observations of these stars can be reasonably well explained by
assuming the presence of the powerful direct Urca process of neutrino emission
in the inner cores of massive stars, introducing a slight broadening of the
direct Urca threshold (for instance, by proton superfluidity), and by tuning
mass distributions of isolated and accreted neutron stars.
