Suppose we have the differential equation
\begin{align}
\dot{\mathbf{x}}&=v(t,\mathbf{x})\tag{0}\label{zero}
\end{align}
where $v:I\times\Omega\subset \mathbb{R}\times\mathbb{R}^n
\longrightarrow\mathbb{R}^n$ satisfies the smoothness conditions of the uniqueness and existence theorem.
An important property of the solutions to the differential
equation \eqref{zero} is the flow property, namely, if $\phi_{t,s}(y)=\phi(t;s,y)$ denotes the solution
in $t$ to \eqref{zero} with
initial conditions $\phi(s;s,y)=y$, then
- $\phi_{t,s}\circ\phi_{s,r}(y)=\phi_{t,r}(y)$ for all
$r,\,s,\,t\in I$ and $y\in \Omega$.
- $\phi_{t,t}(y)=y$ for all $t\in I$ and $y\in\Omega$.
Each $\phi_{s,t}(\cdot)$ defines a diffeomorpism on $\Omega$ to itself, its inverse being $\phi_{t,s}$
Suppose that $\mu$ is a measure on $(\Omega,\mathcal{B}(\Omega))$ with smooth density $g$ with respect to
Lebesgue measure on $\Omega$. The flow $\phi$ induces a family of measures measure
$\mu_t$ on phase space $(\Omega, \mathcal{B}(\Omega))$ given by
\begin{align}
\mu_t(B):=\mu(\phi^{-1}_{t,0}(B))=\mu(\phi_{0,t}(B))\tag{1}\label{push}
\end{align}
Furthermore, $\mu_t$ also has a density with respect to the Lebesgue
measure. This follows from
\begin{align}
\mu_t(B)=\int_{\phi_{0,t}(B)}g({\bf x})d{\bf x}=
\int_{B} g(\phi_{0,t}({\bf x}))
\det\left[\tfrac{\partial\phi_{0,t}}{\partial x}({\bf x})\right]d{\bf x}\tag{2}\label{push_density}.
\end{align}
So, $p(t;\mathbf{x})=g(\phi_{0,t}({\bf x}))
\det\left[\tfrac{\partial\phi_{0,t}}{\partial x}({\bf x})\right]$
is a density for $\mu_t$.
The following result, a generalization of Liouville's conservation law, is the content of the OP:
Therorem Let $\phi$ be the flow associated to equation
\eqref{zero}, and $\mu(d{\bf x})=g({\bf x})d{\bf x}$ be a smooth measure on phase space $(\Omega,\mathcal{B}(\Omega))$. Then, the density
$p:=p(t,{\bf x})$ of
the induced measure $\mu_t(d\mathbf{x})$ is smooth and satisfies the equation
\begin{align}
\frac{\partial p}{\partial t}+\nabla_{\bf x}\cdot(p\,v)=0\tag{3}\label{three}
\end{align}
For the proof of this fact, first I state the following elementary auxiliary results from Calculus and from ordinary differential equations:
Lemma D:
Let $\Delta:\mathbb{R}^{n^2}\longrightarrow\mathbb{R}$ be the
determinant function, i.e.
$$\Delta(\alpha_{11},\ldots,\alpha_{n1},\ldots,\alpha_{1n},
\ldots,\alpha_{nn})^{\top}
= \det[(\alpha_{ij})]$$
where $(\alpha_{ij})$ is the $n\times n$-matrix whose $ij$-th component is
$\alpha_{ij}$. Then,
$$\Delta_\alpha=
\frac{\partial \Delta}{\partial\alpha}=
(W_{11}\ldots,W_{n1},\ldots,W_{1n},\ldots,W_{nn})$$
where $W_{ij}$ is the $ij$-th cofactor of the matrix $(\alpha_{ij})$.
Proposition W:
Let $\phi(t;{\bf x})$ be a solution to the
equation $\dot{{\bf x}}= v(t,{\bf x})$,
with $\phi(0;{\bf x})={\bf x}$.
Define the function $W$ by
\begin{align}
W(t,{\bf x})&=\det\left[\frac{\partial \phi}
{\partial {\bf x}}(t;{\bf x})\right].
\end{align}
Then, $W$ satisfies the differential equation
\begin{align}
\dot{W}(t)=W(t)\, (\nabla_{\bf x}\cdot v)(t,\phi(t;\mathbf{x})); \qquad W(0)=1\tag{4}\label{four}
\end{align}
where $\left(\nabla_{\bf x}\cdot v\right)(t,\phi(t;{\bf x})) =\sum_{j=1}^n \frac{\partial v}{\partial x_j}(t,\phi(t;{\bf x}))$.
Remark: Proof of these results are discussed in this old posting here. I present them again at the end of this post for completeness.
Proof of \eqref{three}:
Notice that $\mu(B)=\mu_t(\phi_{t,0}(B))$ for all $t$ and any
$B\in\mathcal{B}(\Omega)$. Hence, $0=\tfrac{d}{dt}\mu_t(\phi_{t,0}(B))$.
From Proposition W and the multivariable change of variables formula
\begin{align}
0 &= \frac{d}{dt}\int_{\phi_{t,0}(B)} p(t,\mathbf{x})\,d\mathbf{x} =
\frac{d}{dt}\int_{B}p(t,\phi_{t,0}(\mathbf{x}))\det\left[\frac{\partial\phi_{t,0}}{\partial\mathbf{x}}(\mathbf{x})\right]\,d\mathbf{x}\\
&=\int_B \left(\tfrac{\partial p}{\partial t}(t,\phi_{t,0}(\mathbf{x})) + \big(\nabla_{\mathbf{x}}p\cdot v\big)(t,\phi_{t,0}(\mathbf{x}))+ \big(p\,\nabla_{\mathbf{x}}\cdot v\big)(t,\phi_{t,0}(\mathbf{x})\right) \operatorname{det}\Big[\frac{\partial \phi_{t,0}}{\partial x}(\mathbf{x})\Big]\,d\mathbf{x} \\
&=\int_B \left(\frac{\partial p}{\partial t}(t,\phi_{t,0}({\bf x}))
+ \big(\nabla_{\bf x}\cdot(p\,v)\big)(t,\phi_{t,0}({\bf x}))\right)
\det\left[\frac{\partial\phi_{t,0}}{\partial x}({\bf x})\right]\,d{\bf x}\\
&=\int_{\phi_{t,0}(B)}\tfrac{\partial p}{\partial t}(t,{\bf x})
+ \big(\nabla_{\bf x}\cdot(p\,v)\big)(t,{\bf x})\,d{\bf x}\tag{5}\label{five}
\end{align}
for any $B\in \mathcal{B}(\Omega)$. Since for any $t$ fixed, $\phi_t$ is a diffeomorphism, we conclude that
\begin{align}
\tfrac{\partial p}{\partial t}
+ \nabla_{\bf x}\cdot(p\,v)=0
\end{align}
Here are short proof to the main results used in the proof of the main result \eqref{three}.
Proof of Lemma D: Compute determinants using the cofactor formula.
Proof of Proposition W:
If $\phi(t;{\bf x})=(\phi^1(t;{\bf x}),\ldots,\phi^n(t;{\bf x}))^\top$, denote by
$\phi^{i}_{x_j}(t;{\bf x})= \frac{\partial\phi^i}{\partial x_j}(t;{\bf x})$.
Lemma D along with the chain rule yields
\begin{align}
\dot{W}&= \sum_i W_{i1}\dot{\phi}^{i}_{x_1} +\cdots+ \sum_i
W_{in}\dot{\phi}^{i}_{x_n}\\
&=\sum_{ij} W_{ij}\dot{\phi}^{i}_{x_j}
\tag{6}\label{chain}
\end{align}
where $W_{ij}$ is the $ij$--th cofactor of the matrix
$\left(\phi^i_{x_j}\right)$.
It is easy to check that
$\phi_{\bf x}(t;{\bf x}):=\frac{\partial\phi}{\partial
{\bf x}}(t;{\bf x})$ satisfies the variational equation
\begin{align}
\begin{matrix}
\dot{\phi}_{\bf x}(t;{\bf
x})&=&v_{\bf{x}}(t,\phi(t;{\bf{x}}))\phi_{\bf x}(t;{\bf x})\\
\phi_{\bf x}(0;{\bf x})&=&I
\end{matrix}
\tag{7}\label{vareq}
\end{align}
substituting \eqref{vareq} on \eqref{chain} and recalling the fact that
the determinant of a matrix that has two identical columns is zero, we obtain
\begin{align}
\dot{W}(t)&=\sum_{ijk} W_{ij}(t)
v^i_{x_k}(t,\phi(t;{\bf x}))\phi^k_{x_j}(t;{\bf x})\\
&= \sum_{ki} \left(v^i_{x_k}(t,\phi(t;{\bf x})\right)
\sum_j W_{ij}(t)\phi^k_{x_j}(t;{\bf x})\\
&=\sum_i v^i_{x_i}(t,\phi(t;{\bf x})
\sum_j W_{ij}(t)\phi^i_{x_j}(t;{\bf x})\\
&= \sum_i v^i_{x_i}(t,\phi(t;{\bf x})) W(t) =
W(t)\,\left(\nabla_{\bf x}\cdot
v\right)(t,\phi(t;{\bf x}))
\end{align}
