$$\frac{\partial \rho }{\partial t}= -\sum_{i}\left(\frac{\partial \rho}{\partial q_i}\,\dot{q_i}+\frac{\partial\rho}{\partial p_i}\,\dot p_i\right)$$
This means that if we have a function of $t, p, q$ namely $\rho(t,\vec p,\vec q)$ and we have a trajectory that is a curve in $(p,q)$ space, namely $q_i(t), p_i(t), i=1\ldots N,$ then:
$$
\frac{\mathrm d}{\mathrm dt} \rho(t, \vec p(t), \vec q(t)) =\frac{\partial \rho }{\partial t}+ \sum_{i}\left(\frac{\partial \rho}{\partial q_i}\,\dot{q_i}+\frac{\partial\rho}{\partial p_i}\,\dot p_i\right)
$$
How if $\rho$ is constant along trajectories, then LHS is 0 and the equation you have written follows directory.
So:
- a trajectory is any curve in 2N dimensional space described in $q_i$ and $p_i$ coordinates
- $\rho$ is a function of both time and $\vec q$ and $\vec p$
- whole concept is just an application of a chain rule.
Why There is a Need for a Further Axiom
To derive Liouville's equation, you indeed need another axiom further to your assumptions. Something like: "there is no nett creation or destruction of any particle of any species throughout the particle system state evolution". The easiest way to understand the need for this axiom is to cite a system wherefor Liouville's equation cannot hold, even though particles undergo dynamical evolution described by Hamilton's equations throughout their lifetimes: a system of particles undergoing a far-from-equilibrium chemical reaction. In such a system, reactant particle species are consumed by the reaction, and disappear from phase space. Reaction product particles appear in phase space in their place. Moreover, chemical energy is converted to kinetic energy (or contrariwise), so that a product species will "suddenly" appear at a different point in phase space from the one where the correspondingly consumed reactant species particles vanished. Liouville's equations would conceptually be replaced by a coupled system of equations, one for each species $j$, of the form:
$$\frac{\partial\,\rho_j(X,\,t)}{\partial\,t} = \{H,\,\rho_j\} + \sum\limits_k \int_\mathcal{P} M_{j\,k}(X,\,X^\prime)\,\rho_k(X^\prime,\,t)\,\mathrm{d}\Gamma^\prime$$
where the integral is over all phase space $\mathcal{P}$, $\Gamma$ is the measure defined by the volume form and the kernel $M_{j\,k}$ expresses detailed stochimetric balance between the chemically reacting species as well as other physical principles such as conservation of energy, momentum and strict increase with time of entropy. Note that I said "nett" creation or destruction: Liouville's equation would work if the reaction were at equilibrium.
Complete Axioms
The following axioms (1. and 2. are equivalent to yours) will get you Liouville's equation:
- Axiom 1: Phase space is a $2\,N$ dimensional $C^2$ manifold $\mathcal{P}$;
- Axiom 2: Points in phase space always and only evolve with a flow parameter $t$ through Hamilton's equations defined by a $C^2$ Hamiltonian $H:\mathcal{P}\times \mathbb{R}\to\mathbb{R}$, the latter possibly time varying (hence the domain $\mathcal{P}\times\mathbb{R}$);
- Axiom 3: The full states of particles are points in $\mathcal{P}$ evolving according to axiom 2 and there is no nett creation or destruction of any particle of any species throughout the particle system state evolution.
From Complete Axioms to Liouville's Equation
From these axioms, the chain of inference you need is as follows:
- Inference 1: From axioms 1. and 2., deduce that any $X\in T_p\,\mathcal{P};\;\forall p\in \mathcal{P}$ expressed in canonical co-ordinates (i.e. ones for which the Hamilton equations hold) that is Lie-dragged by the Hamiltonian flow evolves according to $\dot{X} = A(t)\,X$ where $A(t)\in\mathfrak{sp}(N,\,\mathbb{R})$, thus the symplectic 2-form $\omega(X,\,Y)\stackrel{def}{=} X^T\,\Omega\,Y$ where, for the special case of canonical co-ordinates, $\Omega =\left(\begin{array}{cc}0&-1_N\\1_N&0\end{array}\right)\;\forall p\in\mathcal{P}$ is conserved under the mapping $\mathcal{P}\mapsto \Phi(H,\,t)\,\mathcal{P},\,\forall t\in\mathbb{R}$ induced by the Hamiltonian flow. (Indeed, at any given point $p\in\mathcal{P}$ find $N$ different $C^2$ Hamiltonians such that the tangents to their flows span $T_p\,\mathcal{P}$ to deduce that the Lie derivative of $\omega$ in any direction is nought, thus $\mathrm{d}\omega=0$ from Cartan's formula relating Lie and Exterior derivatives, but this information is further to our immediate needs). Take heed that inference 1 holds whether or not the Hamiltonian be time-varying. In the latter case, the Hamiltonian is not constant along the flow, but the flow still conserves the symplectic form.
- Inference 2: From inference 1, we have immediately that the volume form $\Gamma = \omega^N$ ($N^{th}$ exterior power) is conserved under Hamiltonian flows. Thus deduce Liouville's theorem (as opposed to equation). Alternatively, the conservation of the symplectic form shown in Inference 1 implies that the Jacobi matrix of the transformation $\mathcal{P}\mapsto \Phi(H,\,t)\mathcal{P}$ is a symplectic matrix (member of $\mathrm{Sp}(N,\,\mathbb{R})$), which always has a unit determinant. Thus the volume form is conserved.
- Inference 3: But the volume form is also the Jacobian of the transformation $\mathcal{P}\mapsto \Phi(H,\,t)\,\mathcal{P}$ and $J(p,\,\Phi(H,\,0))=J(p,\,\mathrm{id})=1$. Since the volume form is conserved, the Jacobian $J(p,\,\Phi(H,\,t))=1,\forall p\in\mathcal{P},\forall t\in \mathbb{R}$. Thus $\Phi$ is everywhere a local bijection (inverse function theorem). Alternatively, we can make the same deduction this straight from Inference 1, which implies that the Jacobi matrix of the transformation $\mathcal{P}\mapsto \Phi(H,\,t)\mathcal{P}$ is a symplectic matrix (member of $SP(N,\,\mathbb{R})$), which is never singular and indeed always has a unit determinant.
- Inference 4: From Axioms 2 and 1, deduce that the distance function defined in canonical co-ordinates by $d(p_1,\,p_2) = (p_1-p_2)^T\,(p_1-p_2)$, zero iff $p_1=p_2$, between any pair of points $p_1,\,p_2\in\mathcal{P}$ must be a continuous function of the flow parameter $t$ (continuous with respect to the topology with basis of open balls defined by this distance function);
- Inference 5: From Inference 3, at any $p\in\mathcal{P}$ and $t\in \mathbb{R}$, there is an open set $\mathcal{U}_p$ small enough such that $\Phi(H,\,t):\mathcal{U}_p\to\Phi(H,\,t)(\mathcal{U}_p)$ is a bijection. The question now arises as to whether $\Phi(H,\,t)$ can map any point outside $\mathcal{U}_p$ into $\Phi(H,\,t)(\mathcal{U}_p)$ (which situation would make $\Phi(H,\,t)$ a local bijection but globally many to one for some $t\in\mathbb{R}$). However, if two or more points are mapped to one point $\tilde{p}\in\Phi(H,\,t)(\mathcal{U}_p)$, from inference 4. deduce that $\exists t$ small enough that the two chosen preimages of $\tilde{p}$ both lie in $\mathcal{U}_p$, thus contradicting local bijectivity. (Informally, from inference 4, multiple points of a function can only arise from mappings along connected "forked" flow lines, so zoom in near enough in on the fork point and thus contradict local bijectivity, showing that forks are impossible). Repeating the reasoning for $-t$ lets us deduce that multiple points are impossible and $\Phi(H,\,t):\mathcal{P}\to\mathcal{P}$ is a global bijection (indeed a symplectomophism in the light of inference 1, but, again, this information is further to our needs);
- Inference 6: From inference 5 and axiom 3, deduce that if there is some number $M$ of particles in any subset $\mathcal{V}\subseteq\mathcal{P}$, then there are precisely $M$ particles in $\Phi(H,\,t)(\mathcal{V})$. From inference 2. deduce that $\mathcal{V}$ and $\Phi(H,\,t)(\mathcal{V})$ have the same volumes. Therefore infer that the average particle density in any subset $\mathcal{V}\subseteq\mathcal{P}$ is constant if the particle states and subsets evolve by Hamiltonian flows;
- Inference 7: Apply inference 6 to a small open set that is shrunken according to an appropriate limiting process to deduce that the density function $\rho(p,\,t)$ at point $p$ and at time $t$ must be the same as the density at point $\Phi(H,\,-\mathrm{d}t)\,p$ at time $t-\mathrm{d}t$. Putting these words into symbols: $\mathcal{L}_{-X}\rho=\frac{\partial\,\rho}{\partial\,t}$, where $X$ is the vector field tangent to the Hamiltonian flow $\Phi$. This is, of course, $\{H,\,\rho\}=\frac{\partial\,\rho}{\partial\,t}$, or Liouville's equation.
Circularity of Other Proofs
Ultimately, I don't believe that proofs of Liouville's equation grounded on the divergence theorem are different from the above: I think that they are tacitly introducing Axiom 3 as "obvious" (even though I hope I have shown at the beginning of my answer that it doesn't always hold) and then the continuity equation and incompressible flows are simply an expression of this tacitly assumed axiom. So I don't think that these "proofs" are circular, just somewhat badly written in making use of tacit assumptions.
Summary
User Image sums all this up nicely (I was perhaps too brainfried to make the last step):
For Axiom 3 however, you showed that Axiom 3 $\Rightarrow \frac{d \varrho}{d t} = 0$. The other direction $\frac{d \varrho}{d t} = 0 \Rightarrow$ Axiom 3 is readily discussed in any textbook (trajectories do not start, end or cross etc.). So in fact we have Axiom 3 $\Leftrightarrow \frac{d \varrho}{d t} = 0$ when we are in the context of Axiom 1+2, e.g. classical mechanics. Hence, Liouville's equation is an axiom.
and indeed, in the presence of the other two, my axiom 3 is logically equivalent to Liouville's equation. My version is perhaps more physically transparent, but open to interpretation, and so the assertion of Liouville's equation as an axiom is perhaps more succinct and precise. So the answer to the title question is that Liouville's Equation must indeed be added as an axiom, and, in the presence of Axioms 1 and 2, it has the meaning that particle number of all species is conserved.
Best Answer
Introduction
Let us define the density of particles of species $s$ in a volume element, $d\mathbf{x} \ d\mathbf{v}$, at a fixed time, $t$, centered at $(\mathbf{x}, \mathbf{v})$ as the quantity $f_{s}(\mathbf{x},\mathbf{v},t)$. I assume this function is non-negative, contains a finite amount of matter, and it exists in the space of positive times and $\mathbb{R}^{3}$ and $\mathbb{R}_{\mathbf{v}}^{3}$, where $\mathbb{R}_{\mathbf{v}}^{3}$ is the space of all possible 3-vector velocities. Then one can see that there are two ways to interpret $f$: (1) it can be an approximation of the true phase space density of a gas (large scale compared to inter-particle separations); or (2) it can reflect our ignorance of the true positions and velocities of the particles in the system. The first interpretation is deterministic while the second is probabilistic. The latter was used implicitly by Boltzmann. Let us assume that $f_{s}(\mathbf{x},\mathbf{v},t)$ $\rightarrow$ $\langle f \rangle + \delta f$, where $\langle f \rangle$ is an ensemble average of $f_{s}$ and I have dropped the subcript out of laziness.
Liouville's Equation
I know that $\langle f \rangle$ satisfies Liouville's equation, or more appropriately, $\partial \langle f \rangle$/$\partial t = 0$. In general, the equation of motion states: $$ \begin{equation} \frac{ \partial f }{ \partial t } = f \left[ \left( \frac{ \partial }{ \partial \textbf{q} } \frac{ d\textbf{q} }{ dt } \right) + \left( \frac{ \partial }{ \partial \textbf{p} } \frac{ d\textbf{p} }{ dt } \right) \right] + \left[ \frac{ d\textbf{q} }{ dt } \cdot \frac{ \partial f }{ \partial \textbf{q} } + \frac{ d\textbf{p} }{ dt } \cdot \frac{ \partial f }{ \partial \textbf{p} } \right] \tag{1} \end{equation} $$ where I have defined the canonical phase space of $(\mathbf{q}, \mathbf{p})$. If I simplify the terms dA/dt to $\dot{A}$ and let $\boldsymbol{\Gamma} = (\mathbf{q}, \mathbf{p})$, then I find: $$ \begin{align} \frac{ \partial f }{ \partial t } & = - f \frac{ \partial }{ \partial \boldsymbol{\Gamma} } \cdot \dot{\boldsymbol{\Gamma}} - \dot{\boldsymbol{\Gamma}} \cdot \frac{ \partial f }{ \partial \boldsymbol{\Gamma} } \tag{2a} \\ & = - \frac{ \partial }{ \partial \boldsymbol{\Gamma} } \cdot \left( \dot{\boldsymbol{\Gamma}} f \right) \tag{2b} \end{align} $$ where one can see that the last form looks like the continuity equation. If I define the total time derivative as: $$ \begin{equation} \frac{ d }{ dt } = \frac{ \partial }{ \partial t } + \dot{\boldsymbol{\Gamma}} \cdot \frac{ \partial }{ \partial \boldsymbol{\Gamma} } \tag{3} \end{equation} $$ then I can show that the time rate of change of the distribution function is given by: $$ \begin{align} \frac{ d f }{ dt } & = \frac{ \partial f }{ \partial t } + \dot{\boldsymbol{\Gamma}} \cdot \frac{ \partial f }{ \partial \boldsymbol{\Gamma} } \tag{4a} \\ & = - \left[ f \frac{ \partial }{ \partial \boldsymbol{\Gamma} } \cdot \dot{\boldsymbol{\Gamma}} + \dot{\boldsymbol{\Gamma}} \cdot \frac{ \partial f }{ \partial \boldsymbol{\Gamma} } \right] + \dot{\boldsymbol{\Gamma}} \cdot \frac{ \partial f }{ \partial \boldsymbol{\Gamma} } \tag{4b} \\ & = - f \frac{ \partial }{ \partial \boldsymbol{\Gamma} } \cdot \dot{\boldsymbol{\Gamma}} \tag{4c} \\ & \equiv - f \Lambda\left( \boldsymbol{\Gamma} \right) \tag{4d} \end{align} $$ where $\Lambda \left( \boldsymbol{\Gamma} \right)$ is called the phase space compression factor. Note that Equations 4a through 4d are different forms of Liouville's equation, which have been obtained without reference to the equations of motion and they do not require the existence of a Hamiltonian. I can rewrite Equation 4d in the following form: $$ \begin{equation} \frac{ d }{ dt } \ln \lvert f \rvert = - \Lambda\left( \boldsymbol{\Gamma} \right) \tag{5} \end{equation} $$
Relation to Hamiltonian
Most readers might not recognize Equations 4d and 5 as Liouville's equation because one usually derives it from a Hamiltonian. If the equations of motion can be generated from a Hamiltonian, then $\Lambda \left( \boldsymbol{\Gamma} \right) = 0$, even in the presence of external fields that act to drive the system away from equilibrium. Note that the existence of a Hamiltonian is a sufficient, but not necessary condition for $\Lambda \left( \boldsymbol{\Gamma} \right) = 0$. For incompressible phase space, I recover the simple form of Liouville's equation: $$ \begin{equation} \frac{ d f }{ dt } = 0 \end{equation} $$ However, Liouville's theorem can be violated by any of the following:
Source of Irreversibility
Irreversibility is somewhat of a conundrum because it arises largely due to our choice of boundary conditions, smoothing assumptions (e.g., coarse graining or mean field theory), and limits. For instance, if I assume a velocity distribution of particles can be represented by a continuous model function, the use of a continuous distribution function inserts irreversibility into the equation. One can argue that this is splitting hairs because it is obvious that irreversibility exists in nature. However, I think it is important because your question points at a deeper issue.
If I assumed perfectly elastic binary particle collisions and ignore quantum uncertainties, one could, in principle, follow the trajectories of all particles in a system forward and backward in time. There would be no irreversibility in this model, if I had strong enough computers. However, binary particle collisions are not truly elastic, so our assumption of elasticity has created a loss of information.
Another subtle point is that Boltzmann a priori defined his, now famous, H-theorem such that time would increase in the correct direction (i.e., positive time). He did not originally relate the H-theorem to entropy, that interpretation came later (I believe with Gibbs, but someone correct me if I am wrong here).
The point is that the concepts of irreversibility and entropy are coupled, but not necessarily through direct means. I am inclined to think that the irreversibility to which you refer arises from our methods of solving the math necessary for modeling dynamical statistical systems.
References