On page 174 of Townsend's "A Modern Approach to Quantum Mechanics", 2nd edition, it says the following:
"For a mixed state, one for which $p_k$ is the probability that a particle is in the state $|\psi^{(k)}\rangle$, then $\hat{\rho}=\sum_kp_k|\psi^{(k)}\rangle\langle\psi^{(k)}|$ where $\sum_kp_k=1$.
[...]
Since the density matrix is Hermitian ($\rho_{ij} = \rho_{ji}^*)$ [which I agree], the density matrix can always be diagonalized with diagonal matrix elements given by the probability $p_k$ [?!]. Thus $tr \hat{\rho}^2=\sum_k p_k^2 \leq 1$."
I do not believe the clause marked [?!] to be true. That would require choosing the states $|\psi^{(k)}\rangle$ as the "basis", but we know that the states $|\psi^{(k)}\rangle$ are neither orthonormal nor is it even a basis (for example, we could have more states than the dimension of the Hilbert space).
Am I correct in my critique above? Moreover, given that the density matrix is indeed Hermitian, what then is the "correct orthonormal basis of eigenvectors", and what do these eigenvectors and their corresponding eigenvalues physically mean?