1. Transport of Differential Forms
Let ##M## be a smooth manifold with local coordinates ##x=(x^{1},\dots,x^{m})^{T}##. Consider a ##k##-form ##\omega## depending on ##t##:
$$\omega(t,x) = \sum_{i_{1}<\dots<i_{k}} \omega_{i_{1}\dots i_{k}}(t,x) dx^{i_{1}} \wedge \dots \wedge dx^{i_{k}}$$
Let ##G_{t_{0}}^{t}: M \to M## be the flow along the trajectories of the system defined by the vector field ##v##:
$$\dot{x} = v(t,x), \quad \frac{d}{dt} G_{t_{0}}^{t}(x) = v(t, G_{t_{0}}^{t}(x)), \quad G_{t_{0}}^{t_{0}}(x) = x$$
The following identity is of key importance for the Lie derivative ##L_v##:
$$\frac{d}{dt} (G_{t_{0}}^{t})^* \omega(t, \cdot) = (G_{t_{0}}^{t})^* \left( \frac{\partial \omega}{\partial t} + L_v \omega \right) \qquad (1)$$
2. The Continuity Equation
Recall that if we have a volume form ##\nu = \rho(t,x) dx^1 \wedge \dots \wedge dx^m##, its Lie derivative is given by:
$$L_v \nu = \text{div}(\rho v) dx^1 \wedge \dots \wedge dx^m, \quad \text{where } \text{div}(\rho v) = \frac{\partial(\rho v^i)}{\partial x^i} \qquad (2)$$
Formula (1) provides a method for solving the following Cauchy problem:
$$\frac{\partial \omega}{\partial t} + L_{v(t,x)} \omega = 0, \quad \omega|_{t=t_0} = \hat{\omega} \qquad (3)$$
Indeed, from (1) it follows that ##\frac{d}{dt}(G_{t_0}^t)^* \omega(t, \cdot) = 0##. Since the pullback does not depend on ##t##, we obtain:
$$(G_{t_0}^t)^* \omega(t, \cdot) = \hat{\omega} \implies \omega(t,x) = ((G_{t_0}^t)^{-1})^* \hat{\omega}$$
3. Derivation of Liouville's Formula
Now, let ##M = \mathbb{R}^m## and consider a linear system where ##v(t,x) = A(t)x,## and let ##A(t)## be an ##m\times m## matrix continuously depending on ##t##.
In this case, the flow is linear: ##G_0^t(x) = X(t)x##, where ##X(t)## is the fundamental matrix:
$$\dot X=AX,\quad X(0)=E.$$
Consider the constant volume form ##\hat{\omega} = \hat{\rho} dx^1 \wedge \dots \wedge dx^m## (with ##\hat{\rho} =\mathrm{const}\neq 0##). Then the form
$$\omega(t) = \rho(t) dx^1 \wedge \dots \wedge dx^m = ((G_0^t)^{-1})^* \hat{\omega} = \det(X^{-1}(t)) \hat{\rho} dx^1 \wedge \dots \wedge dx^m$$
is a solution to the Cauchy problem (3) with ##t_0 = 0##. By formula (2), this is equivalent to the continuity equation:
$$\frac{\partial \rho}{\partial t} + \text{div}(\rho v) = 0, \quad \rho|_{t=0} = \hat{\rho}$$
Since ##\rho(t)## depends only on ##t## and ##v## is linear, we have:
$$\text{div}(\rho v) = \rho(t) \text{tr}(A(t))$$
This leads to the differential equation ##\dot{\rho} + \rho \, \text{tr}(A) = 0##, the solution of which is:
$$\rho(t) = \hat{\rho} e^{-\int_{0}^{t} \text{tr} A(s) ds}$$
Comparing this with our expression ##\rho(t) = \det(X^{-1}(t)) \hat{\rho} = \frac{1}{\det X(t)} \hat{\rho}##, we arrive at the celebrated Liouville formula:
$$\det X(t) = e^{\int_{0}^{t} \text{tr} A(s) ds}$$