I've been reading about the equation, and all of the sources I found state that the equation preserves the trace, self-adjointness, and positivity of the density matrix.
The first two properties are easily verified, but I can't seem to figure out how to prove the positivity ($\langle\alpha|\rho|\alpha\rangle\geq0$) is maintained. So far I wrote down:
$$ \langle\alpha|\frac{d\rho}{dt}|\alpha\rangle= \frac{d}{dt}\left(\langle\alpha|\rho|\alpha\rangle\right)= \frac{d\rho_{\alpha\alpha}}{dt}= -\frac{i}{\hbar}\langle\alpha|H\rho-\rho H|\alpha\rangle=-\frac{i}{\hbar} \sum_\mu \left(H_{\alpha\mu}\rho_{\mu\alpha} - \rho_{\alpha\mu}H_{\mu\alpha} \right) $$
Since both $H$ and $\rho$ are self-adjoint the last expression may be reduced to:
$$ \frac{d\rho_{\alpha\alpha}}{dt}= -\frac{i}{\hbar} \sum_\mu \left(H_{\alpha\mu}\rho_{\mu\alpha} - \left(\rho_{\mu\alpha}H_{\alpha\mu}\right)^* \right) =\frac{2}{\hbar}\text{Im}\left[ \sum_\mu H_{\alpha\mu}\rho_{\mu\alpha} \right] $$
I'm kind of stuck from here, any Help?
Edit:
This part of the question was edited for clarification.
I noticed another peculiar thing, when writing this down. My expressions don't relate to any particular basis. However, if I assume that the basis of choice is the eigenbasis of $\rho$ (or $H$), I trivially get $\frac{d\rho_{\alpha\alpha}}{dt}=0$ as:
$$ \frac{d\rho_{\alpha\alpha}}{dt}= -\frac{i}{\hbar}\langle\alpha|H\rho-\rho H|\alpha\rangle= -\frac{i}{\hbar}\left( \rho_{\alpha\alpha}\langle\alpha|H|\alpha\rangle - \rho_{\alpha\alpha}^*\langle\alpha|H|\alpha\rangle \right)= -\frac{i}{\hbar}\rho_{\alpha\alpha} \left( H_{\alpha\alpha} - H_{\alpha\alpha} \right)=0 $$
What is the meaning of this? As I didn't make any assumptions regarding the eigenbasis of $H$, it means the the diagonal terms of the commutator are zero, and therefore to populations are constant of motion.
This would make sense as for a statistical ensemble in equilibrium we claim $\rho\propto e^{-\beta H}$, which means $\left[ H,\rho\right]=0$, but if so what is the whole point of writing the equation to begin with?
This reasoning is of course flawed as it implies that for every system the diagonal of the density matrix is a constant of motion, which means we cannot alter populations of various states. Where does my reasoning fail? Is the source of the problem the fact that I wasn't careful enough ragarding the commutation of the Hamiltonian with itself at different times?
Edit 2:
I leave the question as it is for the sake others who might repeat my mistake. I found my error in the first step of the derivation. I wrote that:
$$ \langle\alpha|\frac{d\rho}{dt}|\alpha\rangle= \frac{d}{dt}\left(\langle\alpha|\rho|\alpha\rangle\right)= \frac{d\rho_{\alpha\alpha}}{dt} $$
This is of course wrong! Since the Liouville equation is written in the Schrödinger picture the states are time dependent and cannot be incorporated into the time derivative.
What I implicitly did by assuming the it could be incorporated into the time derivative, is I assumed that the state $|\alpha\rangle$, is a stationary state, i.e. that diagonalizes simultaneously $H$ and $\rho$, thus implicitly assuming that they commute, in which case indeed the populations are constants of motion.
In short the calculation of the RHS in the first equation I wrote is correct. However, it has nothing to do the the rate of change of the diagonal elements of $\rho$, and isn't related to population in any way.
However, instead I could write:
$$ \frac{d\rho_{\alpha\alpha}}{dt}= \frac{d}{dt}\left(\langle\alpha|\rho|\alpha\rangle\right)= \left(\frac{d\langle\alpha|}{dt}\right)\rho|\alpha\rangle + \langle\alpha| \left(\frac{d\rho}{dt}\right)|\alpha\rangle + \langle\alpha|\rho \left(\frac{d|\alpha\rangle}{dt}\right)= % \\ % =\frac{i}{\hbar}\left( \langle\alpha|H\rho|\alpha\rangle - \langle\alpha|\rho H|\alpha\rangle \right)+ \langle\alpha| \left(\frac{d\rho}{dt}\right)|\alpha\rangle= % \\ % =\frac{i}{\hbar} \langle\alpha|\left[H,\rho\right]|\alpha\rangle - \frac{i}{\hbar} \langle\alpha|\left[H,\rho\right]|\alpha\rangle=0 $$
Which essentially gives the same information as in the answer below. The last expression however shouldn't be interpreted the time evolution of any specific population. The time propogation of occupation probability is given by:
$$ P_{\alpha\alpha}(t) = |c_\alpha(t)|^2 = \langle\alpha(0)|\rho(t)|\alpha(0)\rangle= \left(G(t)\rho(0)G^\dagger(t)\right)_{\alpha\alpha} $$
where $c_\alpha(t)$ are the time dependent coefficients of some spectral decomposition: $|\psi(t)\rangle = \sum_\alpha c_\alpha(t)|\alpha\rangle$, and $G(t)$ is the time propagator between the states with time difference $t$.
