[proofplan]
We differentiate the defining identity $\iota_{X_g}\omega=dg$ in the direction of the Hamiltonian vector field $X_f$. The Lie derivative commutator identity converts this into a contraction of $\omega$ with the Lie bracket $[X_f,X_g]$, while the Hamiltonian nature of $X_f$ gives $\mathcal L_{X_f}\omega=0$. The right-hand side becomes the differential of $X_f(g)=\{g,f\}$, which vanishes because $\{f,g\}=0$. Nondegeneracy of $\omega$ then forces the vector field $[X_f,X_g]$ to vanish.
[/proofplan]
[step:Show that the Hamiltonian vector field $X_f$ preserves $\omega$]
Since $\omega$ is symplectic, it is closed:
\begin{align*}
d\omega=0.
\end{align*}
By the defining equation for $X_f$,
\begin{align*}
\iota_{X_f}\omega=df.
\end{align*}
Cartan's formula for the Lie derivative of a differential form gives
\begin{align*}
\mathcal L_{X_f}\omega=d(\iota_{X_f}\omega)+\iota_{X_f}(d\omega).
\end{align*}
Substituting the two displayed identities,
\begin{align*}
\mathcal L_{X_f}\omega=d(df)+\iota_{X_f}(0)=0.
\end{align*}
Here $d(df)=d^2f=0$ by [nilpotence of the exterior derivative](/theorems/3913). Thus
\begin{align*}
\mathcal L_{X_f}\omega=0.
\end{align*}
[/step]
[step:Differentiate the identity $\iota_{X_g}\omega=dg$ along $X_f$]
Apply $\mathcal L_{X_f}$ to both sides of
\begin{align*}
\iota_{X_g}\omega=dg.
\end{align*}
The right-hand side satisfies
\begin{align*}
\mathcal L_{X_f}(dg)=d(\mathcal L_{X_f}g)=d(X_f(g)),
\end{align*}
because Lie derivative of a smooth function along $X_f$ is evaluation by $X_f$.
For the left-hand side, use the commutator identity
\begin{align*}
\mathcal L_X(\iota_Y\alpha)=\iota_{[X,Y]}\alpha+\iota_Y(\mathcal L_X\alpha)
\end{align*}
with $X=X_f$, $Y=X_g$, and $\alpha=\omega$. Since $\mathcal L_{X_f}\omega=0$, this gives
\begin{align*}
\mathcal L_{X_f}(\iota_{X_g}\omega)=\iota_{[X_f,X_g]}\omega.
\end{align*}
Therefore
\begin{align*}
\iota_{[X_f,X_g]}\omega=d(X_f(g)).
\end{align*}
[guided]
We now differentiate the equation defining $X_g$:
\begin{align*}
\iota_{X_g}\omega=dg.
\end{align*}
The operation used is the Lie derivative $\mathcal L_{X_f}$, because we want to measure how the one-form $\iota_{X_g}\omega$ changes along the flow direction generated by $X_f$.
On the right-hand side, Lie derivative commutes with exterior differentiation on functions:
\begin{align*}
\mathcal L_{X_f}(dg)=d(\mathcal L_{X_f}g).
\end{align*}
Since $g$ is a smooth real-valued function on $M$, its Lie derivative along $X_f$ is the smooth function $X_f(g):M\to\mathbb R$. Hence
\begin{align*}
\mathcal L_{X_f}(dg)=d(X_f(g)).
\end{align*}
On the left-hand side, the relevant identity is the Lie derivative-contraction commutator formula:
\begin{align*}
\mathcal L_X(\iota_Y\alpha)=\iota_{[X,Y]}\alpha+\iota_Y(\mathcal L_X\alpha)
\end{align*}
for vector fields $X,Y$ and a differential form $\alpha$. We apply it with $X=X_f$, $Y=X_g$, and $\alpha=\omega$. This yields
\begin{align*}
\mathcal L_{X_f}(\iota_{X_g}\omega)=\iota_{[X_f,X_g]}\omega+\iota_{X_g}(\mathcal L_{X_f}\omega).
\end{align*}
The previous step proved $\mathcal L_{X_f}\omega=0$, so the second term vanishes:
\begin{align*}
\mathcal L_{X_f}(\iota_{X_g}\omega)=\iota_{[X_f,X_g]}\omega.
\end{align*}
Since applying $\mathcal L_{X_f}$ to equal one-forms gives equal one-forms, the two computations combine to
\begin{align*}
\iota_{[X_f,X_g]}\omega=d(X_f(g)).
\end{align*}
This is the central identity: it turns the desired vector-field commutator into the differential of a scalar function.
[/guided]
[/step]
[step:Use the vanishing Poisson bracket to make the contraction zero]
By the Hamiltonian convention,
\begin{align*}
dg=\iota_{X_g}\omega.
\end{align*}
Evaluating both sides on $X_f$ gives
\begin{align*}
X_f(g)=dg(X_f)=\omega(X_g,X_f).
\end{align*}
By skew-symmetry of $\omega$,
\begin{align*}
\omega(X_g,X_f)=-\omega(X_f,X_g).
\end{align*}
Using the Poisson bracket convention $\{f,g\}=\omega(X_f,X_g)$, we obtain
\begin{align*}
X_f(g)=-\{f,g\}.
\end{align*}
The hypothesis says $\{f,g\}=0$ as a smooth function on $M$, so
\begin{align*}
X_f(g)=0.
\end{align*}
Taking exterior derivatives,
\begin{align*}
d(X_f(g))=0.
\end{align*}
From the identity established above,
\begin{align*}
\iota_{[X_f,X_g]}\omega=0.
\end{align*}
[/step]
[step:Apply nondegeneracy of $\omega$ pointwise to conclude the bracket vanishes]
Let $p\in M$ be arbitrary. Evaluating
\begin{align*}
\iota_{[X_f,X_g]}\omega=0
\end{align*}
at $p$ gives
\begin{align*}
\omega_p([X_f,X_g]_p,v)=0
\end{align*}
for every $v\in T_pM$. Since $\omega_p:T_pM\times T_pM\to\mathbb R$ is nondegenerate, the only vector in $T_pM$ whose contraction with $\omega_p$ is the zero covector is $0$. Hence
\begin{align*}
[X_f,X_g]_p=0.
\end{align*}
Because $p\in M$ was arbitrary,
\begin{align*}
[X_f,X_g]=0
\end{align*}
as a smooth vector field on $M$.
[/step]
