In many cases, measurements are performed on “mixed” ensembles of realizations of a system, which cannot be associated with a single vector in its Hilbert space. In these cases the state of the system is represented by a proper “density operator.” It is instructive to associate density operators with state vectors in a vector space, termed Liouville’s space, where the Schrödinger equation is reformulated as the Liouville–Von Neumann equation. An important consequence of this equation is that the density operator of a system at equilibrium must commute with its Hamiltonian. Namely, the matrix representation of the density operator in the basis of Hamiltonian eigenstates is diagonal, where the diagonal elements are the relative populations of the system’s Hamiltonian eigenstates. The equilibrium populations are revealed by maximizing the (Von Neumann) entropy, subject to given constraints. We derive the equilibrium density operator explicitly for the cases of canonical and grand canonical ensembles.