Theorem: Liouville's Theorem (Hamiltonian Mechanics)

Consider a Hamiltonian system  characterised by the Hamiltonian $\map \HH {t, \mathbf q, \mathbf p}$, where the state of the system, $\mathbf r$, is described by generalized coordinates $\mathbf q$  and $\mathbf p$, which correspond to $N$ generalized positions and $N$ generalized momenta, respectively.
Let the trajectory, $\map {\mathbf r} t$, be the solution of the intial-value problem given by Hamilton's equations together with the initial condition.
Consider the Hamiltonian system's joint distribution function $\map f {t, \map {\mathbf r} t}$.
Allowing that:
:(i) the flux of the number of particles per unit phase-space volume is a continuously differentiable vector field defined on each and every neighborhood of the phase-space volume
:(ii) the trajectory $\map {\mathbf r} t$ is continuously differentiable with respect to time for all time after the initial time
:(iii) the quantity of particles per unit phase-space volume is a conserved quantity that cannot be created nor destroyed in any neighborhood of the phase-space volume
:(iv) the Hamiltonian has continuous second partial derivatives
then the total time derivative of the joint distribution function is identically zero.
In other words:
:$\dfrac {\map {\d f} {t, \map {\mathbf r} t} } {\d t} = 0$
That is:
:Volume integrals are time invariant in phase space.

Proof:
{{tidy}}
{{MissingLinks}}
{{refactor|level = advanced|There are a number of definitions and constructions in the below which really need to be put onto their own pages, where they can be evoked as need be by means of referencing}}
In this proof we first describe the measure space for the Hamiltonian system.
This formalism is necessary for the remainder of the proof.
Next, we invoke the continuity equation and repeatedly adjust the continuity equation based on the several premises of the theorem.  
We describe the measure space for the Hamiltonian system, which is a triple $\struct {\Omega, \FF, \mu}$.
1.  Let $\Omega_{p_k}$ and $\Omega_{q_k}$ denote the sets:
{{improve|Define the below as a subdivision}}
{{begin-eqn}}
{{eqn | l = \Omega_{q^k}
      | r = \set {\hointr {q^k_0} {q^k_1}, \hointr {q^k_1} {q^k_2}, \ldots, \hointr {q^k_{i_{k - 2} } } {q^k_{i_{k - 1} } }, \closedint {q^k_{i_{k - 1} } } {q^k_{i_k} } }
      | c = 
}}
{{eqn-intertext|and:
}}
{{eqn | l = \Omega_{p_k}
      | r = \set {\hointr {p_k^0} {p_k^1}, \hointr {p_k^1} {p_k^2}, \ldots, \hointr {p_k^{j_{k - 2} } } {p_k^{j_{k - 1} } }, \closedint {p_k^{j_{k - 1} } } {p_k^{j_k} } }
      | c = 
}}
{{end-eqn}}
where:
:for $k = 1, \ldots N$, $i_k$ and $j_k$ are real numbers greater than or equal to $1$
:$q_\text{min}, q^k_0, q^k_1, \ldots, q^k_{i_{k - 2} }, q^k_{i_{k - 1} }, q^k_{i_k}, q_\text{max}$ are fixed real numbers such that:
::$q_\text{min} = q^k_0 < q^k_1 < \cdots < q^k_{i_{k - 2} } < q^k_{i_{k - 1} } < q^k_{i_k} = q_\text{max}$
:$p_\text{min}, p_k^0, p_k^1, \ldots, p_k^{j_{k - 2} }, p_k^{j_{k - 1} }, p_k^{j_k}, p_\text{max}$ are fixed real numbers such that:
::$p_\text{min} = p_k^0 < p_k^1 < \cdots < p_k^{j_{k - 2} } < p_k^{j_{k - 1} } < p_k^{j_k} = p_\text{max}$
Let $\Omega$ denote the sample space, which is the set of all possible outcomes, that is given by the cartesian product in the equation:
:$\Omega = \Omega_{q^1} \times \Omega_{q^2} \times \cdots \times \Omega_{q^n} \times \Omega_{p_1} \times \Omega_{p_2} \cdots \times \Omega_{p_n}$
2. Let $\FF$ denote the $\sigma$-algebra of $\Omega$, which is a set of subsets  of $\Omega$.
In this proof, let: 
:$\FF = \powerset \Omega$
where $\powerset \Omega$ is the power set of $\Omega$.
3. Let $\mu$ denote a finite  measure on $\struct {\Omega, \FF}$ with total measure $N^*$, where $N^*$ is the quantity of particles in the Hamiltonian system.
In other words:
:$\map \mu \Omega = N^*$
For example, if our Hamiltonian is an ideal gas, each particle has 6 degrees of freedom (that is the spatial coordinates $q^x, q^y, q^z$ and the momenta $p_x, p_y, p_z$), then $N$ and $N^*$ are related by the equation $2 N = 6 N^*$.
We use the continuity equation.
Let $\omega_o$ denote an arbitrary but fixed sample point in $\Omega$.
Let $V_o$, where $V_o \subseteq \R^{2 N}$, denote the phase-space volume of $\omega_o$.
In other words:
:$\ds V_o = \int_{q^1_{\paren {\ell_1 - 1} } }^{q^1_{\ell_1} } \cdots \int_{q^N_{\paren {\ell_N - 1} } }^{q^N_{\ell_N} } \int_{p_1^{\paren {m_1 - 1} } }^{p_1^{m_1} } \cdots \int_{p_N^{\paren {m_N - 1} } }^{p_N^{m_N} } \rd^{2 N} \mathbf r$
Let $S_o$ denote the piecewise smooth boundary that encloses the volume $V_o$.
Let $\oiint_{S_o} \rd \mathbf S$ denote a surface integral over the closed surface $S_o$.
Let $\map {N_o} t$ denote the measurable quantity of particles within $V_o$ at time $t$.
This can be written as as $\map \mu {\omega_o}$.
Or, this can be written in terms of a joint distribution function, which can be thought of intuitively as the number of particles per unit phase-space volume that fall within each and every $2N$-dimensional infinitesimal interval:
:$(1): \quad \ds \map {N_o} t = \int_{q^1_{\paren {\ell_1 - 1} } }^{q^1_{\ell_1} } \cdots \int_{q^N_{\paren {\ell_N - 1} } }^{q^N_{\ell_N} } \int_{p_1^{\paren {m_1 - 1} } }^{p_1^{m_1} } \cdots \int_{p_N^{\paren {m_N - 1} } }^{p_N^{m_N} } \map f {t, \mathbf r} \rd^{2 N} \mathbf r$
By $\map {\mathbf j_o} t$, I mean the flux of $N_o $ at time $t$, which I assume to be a continuously differentiable vector field defined on a neighborhood of $V_o$.
By $\Sigma_o$ I mean the net rate that $N_o$ is being generated inside the volume $V$ per unit time.
With these definitions, the integral form of the continuity equation expressing the rate of increase of $N_o$ within the volume $V_o$ is: 
:$\ds \dfrac {\partial N_o} {\partial \d t} + \oiint_S \mathbf j \cdot \d \mathbf S = \Sigma_o$
Since, by construction, $N_o$ is a conserved quantity that cannot be created nor destroyed, therefore, the net rate per unit time that the measurable quantity of particles $N_o$ is being generated in the volume $V_o$ per unit time is zero.
Further, since the trajectory $\map {\mathbf r} t$ is continuously differentiable with respect to time for all time after the initial time, there exists a velocity field $\dot {\mathbf r}$ that describes the relevant flux.
In other words:
:$\mathbf j = \map f {t, \map {\mathbf r} t} \dot {\mathbf r}$
Then, the continuity equation becomes:
:$\ds \frac {\partial N_o} {\partial \d t} + \oiint_S \map f {t, \map {\mathbf r} t} \dot {\mathbf r} \cdot \d \mathbf S = 0$
By taking the partial derivative of Eq. 1 with respect to time, considering that:
:$\d \mathbf S = \hat {\mathbf n} \d S$
where $\mathbf {\hat n}$ is the outward pointing unit normal at each point on the boundary surface, and utilizing the divergence theorem, we find that the continuity equation becomes:
:$(2): \quad \ds \int_{q^1_{\paren {\ell_1 - 1} } }^{q^1_{\ell_1} } \cdots \int_{q^N_{\paren {\ell_N - 1} } }^{q^N_{\ell_N} } \int_{p_1^{\paren {m_1 - 1} } }^{p_1^{m_1} } \cdots \int_{p_N^{\paren {m_N - 1} } }^{p_N^{m_N} } \left[ \frac {\partial f} {\partial t} + \boldsymbol \nabla \cdot \paren {\map f {\dot {\mathbf r} } } \right] \rd^{2 N} \mathbf r = 0$
In the equation immediately above, the gradient operator is defined as:
:$\boldsymbol \nabla \equiv \paren {\dfrac \partial {\partial q^1}, \cdots, \dfrac \partial {\partial q^N}, \dfrac \partial {\partial p_1} \cdots \dfrac  \partial {\partial p_N} }$
Now, note that the divergence $\boldsymbol \nabla \cdot \mathbf j$ may be rewritten by the following steps:
{{begin-eqn}}
{{eqn | l = \map {\boldsymbol \nabla \cdot} {f \, \dot {\mathbf r} }
      | r = \sum_{k \mathop = 1}^N \frac \partial {\partial q^k} \sqbrk {\map f {t, \map {\mathbf r} t} \, \map {\dot q^k} t} + \frac \partial {\partial p_k} \sqbrk {\map f {t, \map {\mathbf r} t} \, \map {\dot p_k} t}
      | c = {{Defof|Divergence}}
}}
{{eqn | r = \sum_{k \mathop = 1}^N \sqbrk {\frac {\map {\partial f} {t, \map {\mathbf r} t} } {\partial q^k} \map {\dot q^k} t + \frac {\map {\partial f} {t, \map {\mathbf r} t} } {\partial p_k} \map {\dot p_k} t} + \sum_{k \mathop = 1}^N \map f {t, \map {\mathbf r} t} \paren {\frac {\map {\partial \dot q^k} t} {\partial q^k} + \frac {\map {\partial \dot p_k} t} {\partial p_k} }
      | c = Product Rule
}}
{{eqn | r = \sum_{k \mathop = 1}^N \sqbrk {\frac {\map {\partial f} t} {\partial q^k} \map {\dot q^k} t + \frac {\map {\partial f} t} {\partial p_k} \map {\dot p_k} t } + \sum_{k \mathop = 1}^N \map f t \paren {\frac {\map {\partial^2 \HH} t} {\partial q^k p^k} - \frac {\map {\partial^2 \HH} t} {\partial p^k q^k} }
      | c = Hamilton's equations
}}
{{eqn | r = \sum_{k \mathop = 1}^N \frac {\map {\partial f} t} {\partial q^k} \map {\dot q^k} t + \frac {\map {\partial f} t} {\partial p_k} \map {\dot p_k} t
      | c = $\HH$ has continuous second partial derivatives
}}
{{end-eqn}}
With this, we we see that the entire integrand (that is, $\dfrac {\partial f} {\partial t} + \boldsymbol \nabla \cdot \mathbf j$) is nothing less and nothing more than the total derivative.
In other words:
:$\ds \dfrac {\map {\partial f} {t, \map {\mathbf r} t} } {\partial t} + \boldsymbol \nabla \cdot \mathbf j = \frac {\map {\partial f} {t, \map {\mathbf r} t} } {\partial t} + \sum_{k \mathop = 1}^N \frac {\map {\partial f} {t, \map {\mathbf r} t} } {\partial q^k} \map {\dot q^k} t + \frac {\map {\partial f} {t, \map {\mathbf r} t} } {\partial p_k} \map {\dot p_k} t = \frac {\map {\d f} {t, \map {\mathbf r} t} } {\d t}$
Thus, the continuity equation is restated as:
:$\ds \int_{q^1_{\paren {\ell_1 - 1} } }^{q^1_{\ell_1} } \cdots \int_{q^N_{\paren {\ell_N - 1} } }^{q^N_{\ell_N} } \int_{p_1^{\paren {m_1 - 1} } }^{p_1^{m_1} } \cdots \int_{p_N^{\paren {m_N - 1} } }^{p_N^{m_N} } \frac {\map {\d f} {t, \map {\mathbf r} t} } {\d t} \rd^{2 N} \mathbf r = 0$
Finally, consider that $\dfrac {\map {\d f} {t, \map {\mathbf r} t} } {\d t}$ is not a pathological function.
Therefore, it is possible to construct a sample space $\Omega$ such that:
:(i) the phase-space volume of each and every sample point is non-zero
:(ii) either $\dfrac {\map {\d f} {t, \map {\mathbf r} t} } {\d t}$ is strictly greater than or equal to zero throughout each and every sample point's phase-space volume or $\dfrac {\map {\d f} {t, \map {\mathbf r} t} } {\d t}$ is strictly less than or equal to zero throughout each and every sample point's phase-space volume.
This being so, it must be that  the total time derivative of  the joint distribution function is identically zero.
In other words:
:$\dfrac {\map {\d f} {t, \map {\mathbf r} t} } {\d t} = 0$
{{qed}}
{{Namedfor|Joseph Liouville|cat = Liouville}}