Solution (source code)

= Solution

For $S_t(z)=\Phi(t,0,z)$, differentiability with respect to initial data gives the variational equation
$$
Y'(t)=D_zb(t,S_t(z))Y(t),\qquad Y(0)=I_{2d},\qquad Y(t)=D_zS_t(z).
$$
The <Jacobi determinant derivative formula> gives
$$
J'(t)=\operatorname{tr}(D_zb(t,S_t(z)))J(t)
=(\operatorname{div}_zb)(t,S_t(z))J(t).
$$
The mixed derivatives in the <Hamiltonian vector field> cancel:
$$
\operatorname{div}_zb
=\sum_{i=1}^d\bigl(\partial_{x_i}\partial_{v_i}H
-\partial_{v_i}\partial_{x_i}H\bigr)=0.
$$
Consequently $J'(t)=0$, and $J(0)=1$ gives \b[preservation of <phase space> volume]:
$$
\boxed{\det D_zS_t(z)=1.}
$$
The same argument applies to $\Phi(t,s)$ for every starting time $s$. This is the <Liouville theorem in Hamiltonian mechanics>. Explicit time dependence of $H$ does not affect the cancellation.