[proofplan] The proof has two ingredients. First, the Hamiltonian structure forces the divergence of $b$ to vanish, because the mixed second derivatives of $H$ cancel pairwise. Second, the determinant of the Jacobian matrix of the flow satisfies Liouville's determinant equation, whose coefficient is precisely the divergence of $b$ along the trajectory. The determinant therefore remains equal to its initial value $1$, and the change of variables formula for the $C^1$ diffeomorphism $\Phi_t$ converts this pointwise Jacobian identity into preservation of [Lebesgue measure](/page/Lebesgue%20Measure). [/proofplan]

[step:Compute the divergence of the Hamiltonian vector field]Write a point of $\mathbb{R}^{2d}$ as $z = (x,p)$ with $x = (x_1,\dots,x_d) \in \mathbb{R}^d$ and $p = (p_1,\dots,p_d) \in \mathbb{R}^d$. For each $i \in \{1,\dots,d\}$, define the component functions \begin{align*} b_i: \mathbb{R}^{2d} \to \mathbb{R}, \quad (x,p) \mapsto \frac{\partial H}{\partial p_i}(x,p) \end{align*} and \begin{align*} b_{d+i}: \mathbb{R}^{2d} \to \mathbb{R}, \quad (x,p) \mapsto -\frac{\partial H}{\partial x_i}(x,p). \end{align*} Since $H \in C^2(\mathbb{R}^{2d})$, the vector field $b$ is $C^1$. Its divergence is therefore the [continuous function](/page/Continuous%20Function) \begin{align*} \operatorname{div} b: \mathbb{R}^{2d} \to \mathbb{R}, \quad (x,p) \mapsto \sum_{i=1}^d \frac{\partial b_i}{\partial x_i}(x,p) + \sum_{i=1}^d \frac{\partial b_{d+i}}{\partial p_i}(x,p). \end{align*} Substituting the component formulas gives \begin{align*} \operatorname{div} b(x,p) = \sum_{i=1}^d \frac{\partial^2 H}{\partial x_i \partial p_i}(x,p) - \sum_{i=1}^d \frac{\partial^2 H}{\partial p_i \partial x_i}(x,p). \end{align*} By equality of mixed partial derivatives for $C^2$ functions (citing a result not yet in the wiki: Equality of Mixed Partial Derivatives), each summand cancels with the corresponding summand in the second sum. Hence \begin{align*} \operatorname{div} b(x,p) = 0 \end{align*} for every $(x,p) \in \mathbb{R}^{2d}$.[/step]

[guided]The point of this step is to use the special Hamiltonian sign pattern. The first $d$ components of $b$ are the $p$-derivatives of $H$, while the last $d$ components are the negatives of the $x$-derivatives of $H$. Thus, when we take the divergence, each $x_i$-derivative of a $p_i$-derivative is paired with the negative $p_i$-derivative of the corresponding $x_i$-derivative. Formally, write $z = (x,p) \in \mathbb{R}^{2d}$ with $x = (x_1,\dots,x_d)$ and $p = (p_1,\dots,p_d)$. For each $i \in \{1,\dots,d\}$, the components of $b$ are \begin{align*} b_i(x,p) = \frac{\partial H}{\partial p_i}(x,p) \end{align*} and \begin{align*} b_{d+i}(x,p) = -\frac{\partial H}{\partial x_i}(x,p). \end{align*} Because $H$ is $C^2$, these component functions are $C^1$, so the divergence of $b$ is defined by differentiating the first $d$ components with respect to the $x$-variables and the last $d$ components with respect to the $p$-variables: \begin{align*} \operatorname{div} b(x,p) = \sum_{i=1}^d \frac{\partial}{\partial x_i}\left(\frac{\partial H}{\partial p_i}\right)(x,p) + \sum_{i=1}^d \frac{\partial}{\partial p_i}\left(-\frac{\partial H}{\partial x_i}\right)(x,p). \end{align*} Therefore \begin{align*} \operatorname{div} b(x,p) = \sum_{i=1}^d \frac{\partial^2 H}{\partial x_i \partial p_i}(x,p) - \sum_{i=1}^d \frac{\partial^2 H}{\partial p_i \partial x_i}(x,p). \end{align*} The hypothesis $H \in C^2(\mathbb{R}^{2d})$ is exactly what permits the use of equality of mixed partial derivatives for $C^2$ functions (citing a result not yet in the wiki: Equality of Mixed Partial Derivatives). Thus \begin{align*} \frac{\partial^2 H}{\partial x_i \partial p_i}(x,p) = \frac{\partial^2 H}{\partial p_i \partial x_i}(x,p) \end{align*} for every $i \in \{1,\dots,d\}$ and every $(x,p) \in \mathbb{R}^{2d}$. Each term in the first sum cancels the corresponding term in the second sum, so \begin{align*} \operatorname{div} b(x,p) = 0. \end{align*} This is the infinitesimal volume-preservation property of Hamiltonian vector fields.[/guided]

[step:Apply Liouville's determinant equation to the flow Jacobian] Fix $z \in U_t$. For each $s$ between $0$ and $t$, let \begin{align*} M_s(z) := J\Phi_s(z) \in \mathbb{R}^{2d \times 2d} \end{align*} denote the Jacobian matrix of the flow map $\Phi_s$ at $z$, where the Jacobian matrix is the matrix representation of the total derivative $D\Phi_s{}_z: \mathbb{R}^{2d} \to \mathbb{R}^{2d}$ in the standard bases. Let $I_{2d} \in \mathbb{R}^{2d \times 2d}$ denote the identity matrix. Since $\Phi_0$ is the identity map on its domain, \begin{align*} M_0(z) = I_{2d} \end{align*} and hence \begin{align*} \det M_0(z) = 1. \end{align*} The ODE flow Jacobian determinant formula, also called Liouville's formula for the determinant of the variational equation (citing a result not yet in the wiki: ODE Flow Jacobian Determinant Formula), applies because $b$ is $C^1$ and the assumed flow regularity gives a $C^1$ flow. It gives \begin{align*} \frac{d}{ds}\det M_s(z) = \operatorname{div} b(\Phi_s(z)) \det M_s(z) \end{align*} for every $s$ between $0$ and $t$. By the previous step, $\operatorname{div} b(\Phi_s(z)) = 0$, so \begin{align*} \frac{d}{ds}\det M_s(z) = 0. \end{align*} Therefore $s \mapsto \det M_s(z)$ is constant on the interval between $0$ and $t$. Since its value at $s = 0$ is $1$, we obtain \begin{align*} \det J\Phi_t(z) = \det M_t(z) = 1. \end{align*} This holds for every $z \in U_t$. [/step]

[step:Use change of variables to convert the Jacobian identity into measure preservation] Let $A \subseteq U_t$ be a Borel set. Since $\Phi_t: U_t \to U_{-t}$ is a $C^1$ diffeomorphism, the change of variables formula for $C^1$ diffeomorphisms (citing a result not yet in the wiki: Change of Variables Formula for $C^1$ Diffeomorphisms) applies to the nonnegative Borel function \begin{align*} \mathbb{1}_{\Phi_t(A)}: U_{-t} \to \mathbb{R}, \quad y \mapsto \mathbb{1}_{\Phi_t(A)}(y). \end{align*} Because diffeomorphisms map Borel sets to Borel sets, $\Phi_t(A)$ is Borel in $U_{-t}$. The change of variables formula gives \begin{align*} \int_{U_{-t}} \mathbb{1}_{\Phi_t(A)}(y)\, d\mathcal{L}^{2d}(y) = \int_{U_t} \mathbb{1}_{\Phi_t(A)}(\Phi_t(z)) |\det J\Phi_t(z)|\, d\mathcal{L}^{2d}(z). \end{align*} For every $z \in U_t$, the previous step gives $\det J\Phi_t(z) = 1$, hence $|\det J\Phi_t(z)| = 1$. Also, \begin{align*} \mathbb{1}_{\Phi_t(A)}(\Phi_t(z)) = \mathbb{1}_A(z) \end{align*} because $\Phi_t$ is injective. Therefore \begin{align*} \mathcal{L}^{2d}(\Phi_t(A)) = \int_{U_t} \mathbb{1}_A(z)\, d\mathcal{L}^{2d}(z) = \mathcal{L}^{2d}(A). \end{align*} Since $A \subseteq U_t$ was an arbitrary Borel set, $\Phi_t$ preserves $2d$-dimensional Lebesgue measure on its domain. This proves the theorem. [/step]