Statistical mechanics is used to describe systems with large number of particles ~$10^{23}$.
The observable universe contains between $10^{22}$ to $10^{24}$ stars. Can we treat those many stars as a statistical mechanical system (for which one can define an entropy, temperature..etc)?
.
Yes. In general relativity and cosmology, the collection of galaxies is often even treated as an ideal fluid, in the thermodynamic sense, with temperature and pressure.
The standard textbook on general relativity, the book Gravitation by Misner, Thorne, and Wheeler, discusses this model in Section 27.2, though only on the thermodynamic level.
For the statistical mechanics of gravitation, see, e.g.,
W. Thirring,
Systems with negative specific heat,
Z. Physik 235 (1970), 339-352.
This is for inside a star, but a similar analysis works on the level of galaxies. See, e.g.,
Ahmad et al.,
Statistical Mechanics of the Cosmological Many-Body Problem,
The Astrophysical Journal 571 (2002), 576-584.
A free online copy is at
http://iopscience.iop.org/0004-637X/571/2/576
If you're prepared to regard the application of the virial theorem to large numbers of particles as statistical mechanics, then the answer to your question is definitely yes, because that's how Zwicky worked out that dark matter must be present in the Coma cluster of galaxies. See this article for details, or a quick Google on "Zwicky virial theorem" will find many such articles.
