Preface
The title of the OP's question is "A proof of Liouville’s theorem". Yet, the precise question asked therein by the OP is very specific as it relates to a specific proof of Liouville's theorem. The proof referred to by the OP is not cited by the OP. Separately to these issues, I have found several proofs of Liouville's theorem. All of those are unsuitable for my purposes. For example, the proof in [1] is, in my opinion, too terse, and thus not sufficient for my purposes. There are proofs elsewhere on the web that vary widely both with regards to the actual statement of the theorem, the complexity of the approach, and comprehensiveness of the proof. Because of link rot, I am not posting any cross references to these. All of these issues leave me thinking that I should produce my own statement of the theorem and my own proof. I post here to ask for verification and to welcome suggestions for how to improve my answer. Dear reader, please be kind. Specifically, I entreat you to not down vote my answer unless you are willing to post a comment so that I may later use your comment to improve my answer.
Liouville's theorem
Consider a Hamiltonian system characterised by the Hamiltonian $ \mathcal{H}\!\left(t,\boldsymbol {q},\boldsymbol {p} \right)$, where the state of the system, $\mathbf{r}$, is described by generalized coordinates $\boldsymbol {q}$ and $\boldsymbol {p}$, which correspond to $N$ generalized postions and $N$ generalized momenta, respectively. Let the trajectory, $\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 $f\!\left(t,\mathbf{r}(t)\right)$. 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 $\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 and (iv) the Hamiltonian has continuous second partial derivatives, then the total time derivative of the joint distribution function is identically zero. In other words
$$
\frac{df\!\left(t,\mathbf{r}(t)\right)}{dt} =0.
$$
Proof
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 continuity equation based on the several premises of the theorem.
I now describe the measure space for the Hamiltonian system, which is a triple $\left(\Omega, \mathcal{F}, \mu\right)$.
By $\Omega_{p_k}$ and $\Omega_{q_k}$ I respectively denote the sets
\begin{align*}
\Omega_{q^k}
&=
\left\{
[q^k_{{0}},q^k_{ {1}}),
[q^k_{ {1}},q^k_{ {2}}),
\ldots,
[q^k_{i_{k}-2},q^k_{i_{k}-1}),
[q^k_{i_{k}-1},q^k_{i_{k} }]
\right\}~\text{and}
\\
\Omega_{p_k}
&=
\left\{
[p_k^{{0}},p_k^{ {1}}),
[p_k^{ {1}},p^k_{ {2}}),
\ldots,
[p_k^{{k-2}},p_k^{j_{k}-1 }),
[p_k^{j_{k}-1},p_k^{j_{k}}]
\right\},
\end{align*}
where $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} ; $$
and
$$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} . $$
By $\Omega$, I mean 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} .
$$
By $\mathcal {F}$, I mean the $\sigma$-algebra of $\Omega$, which is a set of subsets of $\Omega$. In this proof I set
$
\mathcal {F} = \mathcal{P}(\Omega)$, where $\mathcal{P}(\Omega)$ is the power set on $\Omega$.
By $\mu$, I mean a finite measure on $\left(\Omega,\mathcal{F}\right)$ with total measure $N^*$, where $N^*$ is the quantity of particles in the Hamiltonian system. In other words
$ \mu{(\Omega)} =N^*$. For example, if our Hamiltonian is an ideal gas, each particle has 6 degrees of freedom (i.e. 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 $2N = 6N^*$.
Now that we have defined the measure space we proceed to prove Liouville's theorem by utilizing the continuity equation.
By $\omega_o$ I mean an arbitrary but fixed sample point in $\Omega$, and
by $V_o$, where $V_o\subset \mathbb{R}^{2N}$, I mean the phase-space volume of $\omega_o$. In other words,
$$
V_o =
\int_{q^1_{ \left(\ell_{1}-1\right) }}^{ q^1_{\ell_{1} }}
\cdots
\int_{q^N_{ \left(\ell_{N}-1\right) }}^{ q^N_{\ell_{N} }}
\int_{p_1^{ \left(m_{1}-1\right) }}^{ p_1^{m_{1} }}
\cdots
\int_{p_N^{ \left(m_{N}-1\right) }}^{ p_N^{m_{N} }}
d^{2N}\mathbf{r}.
$$
By $S_o$ I mean the piecewise smooth boundary that encloses the volume $V_o$.
By $\oiint_{S_o} d\mathbf{S}$ I mean a surface integral over the closed surface $S_o$.
By $N_o(t)$ I mean the measurable quantity of particles within $V_o$ at time $t$. I can write this as
$\mu\left(\omega_o\right)$. Or, I can write this 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 infinitesmal interval,
\begin{equation}
N_o(t) =
\int_{q^1_{ \left(\ell_{1}-1\right) }}^{ q^1_{\ell_{1} }}
\cdots
\int_{q^N_{ \left(\ell_{N}-1\right) }}^{ q^N_{\ell_{N} }}
\int_{p_1^{ \left(m_{1}-1\right) }}^{ p_1^{m_{1} }}
\cdots
\int_{p_N^{ \left(m_{N}-1\right) }}^{ p_N^{m_{N} }}
f(t,\mathbf{r}) \,
d^{2N}\mathbf{r} .
\tag{1}
\end{equation}
By $\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:
$$
\frac{\partial N_o}{\partial dt}
+
\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 $\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} = f(t,\mathbf{r}(t))\dot{\mathbf{r}}$. Then, the continuity equation becomes
$$
\frac{\partial N_o}{\partial dt}
+
\oiint_S f\!\left(t,\mathbf{r}(t)\right)\,\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}} \mathrm {d} S$, where $\mathbf {\hat {n}}$ is the outward pointing unit normal at each point on the boundary surface, and utilizing the the divergence theorem, we find that the continuity equation becomes
\begin{equation}
\label{eq:retyuyiyuiyuy}
\int_{q^1_{ \left(\ell_{1}-1\right) }}^{ q^1_{\ell_{1} }}
\cdots
\int_{q^N_{ \left(\ell_{N}-1\right) }}^{ q^N_{\ell_{N} }}
\int_{p_1^{ \left(m_{1}-1\right) }}^{ p_1^{m_{1} }}
\cdots
\int_{p_N^{ \left(m_{N}-1\right) }}^{ p_N^{m_{N} }}
\left[
\frac{\partial f
}{\partial t}
+
\boldsymbol{\nabla} \cdot
\left(f \,\dot{\mathbf{r}}\right)
\right] d^{2N}\mathbf{r}
=0.
\end{equation}
In
the equation immediately above, the gradient operator is defined as
$$
\boldsymbol{\nabla} \equiv \left(
\frac{\partial}{\partial q^1},\cdots, \frac{\partial}{\partial q^N},\frac{\partial}{\partial p_1} \cdots \frac{\partial}{\partial p_N}
\right).
$$
Now, note that the divergence $ \boldsymbol{\nabla} \cdot \mathbf{j}
$ may be rewritten by the following steps
\begin{align*}
\boldsymbol{\nabla} \cdot
\left(f \,\dot{\mathbf{r}}\right)
&=
\sum_{k=1}^N
\frac{\partial}{\partial q^k} \left[f(t,\mathbf{r}(t))\,\dot{q}^k(t)\right]
+
\frac{\partial}{\partial p_k} \left[f(t,\mathbf{r}(t))\,\dot{p}_k(t)\right]
&&
\text{def. divergence}
\\
&=
\sum_{k=1}^N
\left[
\frac{\partial f(t,\mathbf{r}(t))}{\partial q^k} \dot{q}^k(t)
+
\frac{\partial f(t,\mathbf{r}(t))}{\partial p_k} \dot{p}_k(t)
\right]
+
\sum_{k=1}^N
f(t,\mathbf{r}(t))\left(
\frac{\partial \dot{q}^k(t)}{\partial q^k}
+
\frac{\partial \dot{p}_k(t)}{\partial p_k} \right)
&&
\text{product rule}
\\
&=
\sum_{k=1}^N
\left[
\frac{\partial f(t)}{\partial q^k} \dot{q}^k(t)
+
\frac{\partial f(t)}{\partial p_k} \dot{p}_k(t)
\right]
+
\sum_{k=1}^N
f(t) \left(
\frac{\partial^2 \mathcal{H}(t)}{\partial q^k p^k}
-
\frac{\partial^2 \mathcal{H}(t)}{\partial p^k q^k} \right)
&&
\text{Hamilton's equations}
\\
&=
\sum_{k=1}^N
\frac{\partial f(t)}{\partial q^k} \dot{q}^k(t)
+
\frac{\partial f(t)}{\partial p_k} \dot{p}_k(t)
&&
\mathcal{H}~\text{has continuous second partial derivatives}.
\end{align*}
With, this we we see that the entire integrand (i.e., $\frac{\partial f
}{\partial t}
+
\boldsymbol{\nabla} \cdot \mathbf{j}
$) is nothing less and nothing more than the total derivative. In other words,
$$
\frac{\partial f \!\left(t,\mathbf{r}(t) \right)
}{\partial t}
+
\boldsymbol{\nabla} \cdot \mathbf{j}
=
\frac{\partial f \!\left(t,\mathbf{r}(t) \right)
}{\partial t}
+
\sum_{k=1}^N
\frac{\partial f\!\left(t,\mathbf{r}(t) \right)}{\partial q^k} \dot{q}^k(t)
+
\frac{\partial f\!\left(t,\mathbf{r}(t) \right)}{\partial p_k} \dot{p}_k(t)
= \frac{df\!\left(t,\mathbf{r}(t) \right)}{dt} .
$$
Thus, the continuity equation is restated at
$$
\int_{q^1_{ \left(\ell_{1}-1\right) }}^{ q^1_{\ell_{1} }}
\cdots
\int_{q^N_{ \left(\ell_{N}-1\right) }}^{ q^N_{\ell_{N} }}
\int_{p_1^{ \left(m_{1}-1\right) }}^{ p_1^{m_{1} }}
\cdots
\int_{p_N^{ \left(m_{N}-1\right) }}^{ p_N^{m_{N} }}
\frac{df\!\left(t,\mathbf{r}(t) \right)}{dt} d^{2N}\mathbf{r}
=0.
$$
Finally, consider that $ \frac{df\!\left(t,\mathbf{r}(t) \right)}{dt} $ is not a pathological function. By this I mean that 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 and (ii) either $\frac{df\!\left(t,\mathbf{r}(t) \right)}{dt} $ is strictly greater than or equal to zero throughout each and every sample point's phase-space volume or $\frac{df\!\left(t,\mathbf{r}(t) \right)}{dt} $ 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,
$$
\frac{df\!\left(t,\mathbf{r}(t) \right)}{dt}
=0.
$$
Q.E.D.
Bibliography
[1] Huang, "Statistical Mechanics," John Wiley and Sons, 1987, p 64.
