[proofplan]
The Hilbert form can be written in Hamiltonian coordinates as
\begin{align*}
\omega=-H(x,y,P)\,dx+P\cdot dy.
\end{align*}
If it is closed, exactness of closed $1$-forms on simply connected domains gives a potential $S$ with $dS=\omega$. Comparing the $dy$ and $dx$ components gives $\nabla_yS=P$ and the Hamilton-Jacobi equation. Conversely, any such $S$ satisfies $dS=\omega$, so $d\omega=d^2S=0$.
[/proofplan]
[step:Rewrite the Hilbert form in Hamiltonian coordinates]
By definition,
\begin{align*}
\omega=L(x,y,q(x,y))\,dx+P(x,y)\cdot(dy-q(x,y)\,dx).
\end{align*}
Expanding the right-hand side gives
\begin{align*}
\omega=\bigl(L(x,y,q(x,y))-P(x,y)\cdot q(x,y)\bigr)\,dx+P(x,y)\cdot dy.
\end{align*}
The Hamiltonian identity along the field is
\begin{align*}
H(x,y,P(x,y))=P(x,y)\cdot q(x,y)-L(x,y,q(x,y)).
\end{align*}
Substituting this identity gives
\begin{align*}
\omega=-H(x,y,P(x,y))\,dx+P(x,y)\cdot dy.
\end{align*}
[/step]
[step:Integrate a closed Hilbert form to a Hamilton-Jacobi potential]
Assume $d\omega=0$. The coefficient $P$ is $C^1$ by hypothesis, and the coefficient $-H(x,y,P(x,y))$ is equal to $L(x,y,q(x,y))-P(x,y)\cdot q(x,y)$, which is $C^1$ because $L$ is $C^2$ and $q$ and $P$ are $C^1$. Thus the coefficients of $\omega$ are $C^1$. Since $D$ is simply connected, [Closed 1-Forms on Simply Connected Domains are Exact](/theorems/3618) gives a $C^2$ function $S:D\to\mathbb{R}$ such that
\begin{align*}
dS=\omega.
\end{align*}
Writing
\begin{align*}
dS=\partial_xS\,dx+\nabla_yS\cdot dy
\end{align*}
and comparing with
\begin{align*}
\omega=-H(x,y,P(x,y))\,dx+P(x,y)\cdot dy
\end{align*}
gives
\begin{align*}
\nabla_yS(x,y)=P(x,y)
\end{align*}
and
\begin{align*}
\partial_xS(x,y)=-H(x,y,P(x,y)).
\end{align*}
Substituting $P=\nabla_yS$ into the second identity gives
\begin{align*}
\partial_xS(x,y)+H(x,y,\nabla_yS(x,y))=0.
\end{align*}
[guided]
First rewrite the Hilbert form directly. By definition,
\begin{align*}
\omega=L(x,y,q(x,y))\,dx+P(x,y)\cdot(dy-q(x,y)\,dx).
\end{align*}
Expanding gives
\begin{align*}
\omega=\bigl(L(x,y,q(x,y))-P(x,y)\cdot q(x,y)\bigr)\,dx+P(x,y)\cdot dy.
\end{align*}
Since
\begin{align*}
H(x,y,P(x,y))=P(x,y)\cdot q(x,y)-L(x,y,q(x,y)),
\end{align*}
this becomes
\begin{align*}
\omega=-H(x,y,P(x,y))\,dx+P(x,y)\cdot dy.
\end{align*}
The coefficient $P$ is $C^1$. The coefficient $-H(x,y,P(x,y))$ is the same as $L(x,y,q(x,y))-P(x,y)\cdot q(x,y)$, so it is $C^1$ because $L$ is $C^2$ and $q$ and $P$ are $C^1$. Thus $\omega$ is a $C^1$ one-form. The closedness condition $d\omega=0$ is exactly the integrability condition for finding a potential. Since $D$ is simply connected, the exactness theorem for closed $1$-forms applies and gives $S$ with
\begin{align*}
dS=\omega.
\end{align*}
Now the whole argument is coefficient comparison. The differential of $S$ is
\begin{align*}
dS=\partial_xS\,dx+\nabla_yS\cdot dy.
\end{align*}
The Hilbert form in Hamiltonian coordinates is
\begin{align*}
\omega=-H(x,y,P(x,y))\,dx+P(x,y)\cdot dy.
\end{align*}
Equality of these two $1$-forms forces the $dy$ coefficients to agree, so
\begin{align*}
\nabla_yS=P.
\end{align*}
It also forces the $dx$ coefficients to agree, so
\begin{align*}
\partial_xS=-H(x,y,P).
\end{align*}
Replacing $P$ by $\nabla_yS$ gives the Hamilton-Jacobi equation
\begin{align*}
\partial_xS+H(x,y,\nabla_yS)=0.
\end{align*}
[/guided]
[/step]
[step:Differentiate an existing Hamilton-Jacobi potential back to closedness]
Conversely, suppose $S\in C^2(D)$ satisfies
\begin{align*}
\nabla_yS(x,y)=P(x,y)
\end{align*}
and
\begin{align*}
\partial_xS(x,y)+H(x,y,\nabla_yS(x,y))=0.
\end{align*}
Then
\begin{align*}
\partial_xS(x,y)=-H(x,y,P(x,y)).
\end{align*}
Therefore
\begin{align*}
dS=-H(x,y,P(x,y))\,dx+P(x,y)\cdot dy.
\end{align*}
By the Hamiltonian-coordinate expression for the Hilbert form, the right-hand side is $\omega$. Hence $\omega=dS$. Since exterior differentiation satisfies $d^2S=0$ for every $C^2$ function $S$, it follows that
\begin{align*}
d\omega=d(dS)=0.
\end{align*}
[/step]
