In quantum statistical mechanics, the equilibrium state is characterized by a density matrix $\rho$. Let me focus on the grand canonical ensemble, although the question also holds for the canonical ensemble. In the grand canonical ensemble, this density matrix is: $$\rho(\beta,\mu) = \frac{1}{Z(\beta,\mu)}e^{-\beta (H-\mu N)}$$ where $H$ is the Hamiltonian operator of the system, $N$ is the number operator and $\mu$ is the chemical potential. Of course, we are assuming here that $H-\mu N$ is such that $\rho(\beta,\mu)$ is trace class.
Suppose $H$ has a pure point spectrum and we can find a (countable) basis of eigenstates of $H-\mu N$, each with eigenvalue $\epsilon_{n}$. In this case, at equilibrium the state can be in any of these states with probability $\frac{1}{Z(\beta,\mu)}e^{-\beta \epsilon_{n}}$. This is in agreement with classical statistical mechanics.
However, as far as I understand, the principle that $\rho(\beta,\mu)$ determines the equilibrium state is more general than that. What I mean is that, even if $H-\mu N$ does not have a basis of eigenstates, but $e^{-\beta (H-\mu N)}$ is trace class, then the equilibrium state is described by $\rho(\beta,\mu)$. Being $\rho(\beta,\mu)$ trace class (thus, compact) it has an expansion of the form: $$\rho(\beta,\mu) = \sum_{n\in \mathbb{N}}\lambda_{n}|\psi_{n}\rangle\langle \psi_{n}|$$ where $\rho(\beta,\mu)\psi_{n} = \lambda_{n}\psi_{n}$, $\lambda_{n} \to 0$ and: $$\sum_{n\in \mathbb{N}}\lambda_{n} = 1.$$ In this case, the interpretation of $\lambda_{n}$ as the probability of the system to be at the state $\psi_{n}$ still holds. However, what physical meaning does each $\lambda_{n}$ have? In other words, what it represents being an eigenvalue of $e^{-\beta(H-\mu N)}$? Is it even happen in quantum mechanics (I mean, at least this is possible mathematically).
