density operator - Does a statistical system go into a pure state as the temperature $Tto 0$ (or $betatoinfty$)?

The density matrix $\hat{\rho}$ for a canonical ensemble is given by $$\hat{\rho}=\frac{\sum\limits_{n}e^{-\beta E_n}|n\rangle\langle n|}{Z}\tag{1}$$ $$=\frac{e^{-\beta E_0}}{Z}\Big[\sum\limits_{n=0}|0\rangle\langle 0|+\sum\limits_{n=1}e^{-\beta(E_1-E_0)}|1\rangle\langle 1|+\sum\limits_{n=2}e^{-\beta(E_2-E_0)}|2\rangle\langle 2|+...\Big]\tag{2}$$ where the sum over $n$, represents the sum over quantum states, not energy levels. Also note that since $E_0$ is the ground state enrgy, the difference $(E_n-E_0)>0$, for all $n>0$.

As $T=0$, or equivalently, $\beta\to\infty$, only the first term in (2)contributes dominantly. The density matrix $\hat{\rho}$ goes to $$\hat{\rho}\to\frac{1}{Z}g(E_0)e^{-\beta E_0}|E_0\rangle\langle E_0|\tag{3}$$ where partition function $Z$ goes to $$Z\to g(E_0)e^{-\beta E_0}\tag{4}$$ Here, $g(E_n)$ represents the degeneracy of the $n^{th}$ energy level with energy eigenvalue $E_n$. With (1) and (2), the density matrix simplifies to $\hat{\rho}=|E_0\rangle\langle E_0|$ so that $\hat{\rho}^2=\hat{\rho}$. Therefore, the calculation implies that at $T=0$, the system goes to a pure state.

I don't understand this physically. There will be degeneracies in the ground state. Then why should the system settle down in a pure state? And which degenerate state will it settle down? Intuitively, at $T=0$, the system should be described by an ensemble every member of which are in the ground energy level but in different degenerate states.