In various texts, I see Liouville's theorem stated both verbally and as an equation. But it seems to me that these two formulations don't agree.
For example, in Wikipedia:
The distribution function is constant along any trajectory in phase space.
And:
Consider a Hamiltonian dynamical system with canonical coordinates $q_i$ and conjugate momenta $p_i$, where $i=1,\dots,n$. Then the phase space distribution $\rho(p,q)$ determines the probability $\rho(p,q)\; > \mathrm{d}^nq\,\mathrm{d}^n p$ that the system will be found in the infinitesimal phase space volume $\mathrm{d} ^nq\,\mathrm{d}^n p$. The Liouville equation governs the evolution of $\rho(p,q;t)$ in time $t$:
$\frac{d\rho}{dt}=\frac{\partial\rho}{\partial t}+\sum_{i=1}^n\left(\frac{\partial\rho}{\partial q_i}\dot{q}_i+\frac{\partial\rho}{\partial p_i}\dot{p}_i\right)=0.$
I think that the equation is not right, whereas the verbal description is correct.
To see this, first please confirm that you agree with the following statement, because this might be the source of my confusion:
Liouville's theorem applies to any distribution function $\rho (p,q)$. No assumptions are made on $\rho (p,q)$, as long it is a well behaved function that can represent a distribution. For example, it can be any probability density function (even though I'm not sure it is necessary that it is normalized). In particular, $\rho (p,q)$ doesn't need to be a canonical ensemble, micro-canonical ensemble, etc.. It doesn't need to correspond to thermodynamic equilibrium, or to have well-defined macroscopic values such as volume. $\rho$ is also allowed to have (explicit) time dependence but we ignore this here for simplicity, so $\frac{\partial\rho}{\partial t}=0$.
If you agree with the last paragraph, then it's easy to see that Liouville's equation is wrong:
Let's look at a system with a single point particle with mass $m$ in one dimension, without any forces. So, a phase space point is given by two numbers $(q,p)$. And, given that at time $t_0$ the system is in state $(q_0,p_0)$, the state at any other time $t$ is given by: $\left(q_0+\frac{p_0}{m}(t-t_0),p_0\right)$. Now, let's take $\rho (p,q)$ to be a 2d Gaussian distribution: $\rho (p,q)=Ce^{-(q^2+p^2)}$, ($C$ is a normalization constant).
Let's look at a trajectory with $q_0=0,p_0=1,t_0=0$, and let's take $m=1$. For this trajectory, we have:
$\rho(p(t), q(t))=Ce^{-((q_0+p_0 t)^2+p_0^2)}$
Obviously, the derivative of $Ce^{-((q_0+p_0 t)^2+p_0^2)}$ with respect to $t$ does not vanish, in contradiction with Liouville's equation.
Perhaps, my mistake is that I simply plug-in the trajectories $q(t),p(t)$ into $\rho(q,p)$? I understand that this plugging-in does not represent what the theorem is trying to say, because this misses the flow of the distribution itself, i.e., it doesn't take into account the dynamics of the swarm of representative points that make up the distribution. But, taking this into account would require a different notation for the equation, right? What would that notation be? Or is the notation used in the Wikipedia page correct and I'm just not interpreting it properly?