[proofplan]
After putting the quadratic part in oscillator form, we work only with the Taylor polynomial of $H$ through degree $N$. The Poisson bracket with the quadratic oscillator Hamiltonian acts diagonally on complex monomials, with eigenvalue determined by the integer frequency combination $(\alpha-\beta)\cdot\omega$. The nonresonance hypothesis makes this homological operator invertible on every non-action monomial of degree at most $N$, while the action monomials are exactly its kernel. We then remove the non-action part degree by degree using time-one Hamiltonian flows generated by homogeneous polynomials, and the accumulated error is absorbed into a remainder vanishing to order $N+1$.
[/proofplan]
[step:Put the quadratic Hamiltonian in oscillator coordinates]
By hypothesis, after a linear symplectic change of coordinates on $\mathbb{R}^{2n}$, the degree-two Taylor polynomial of $H$ at $0$ is
\begin{align*}
H_2(Q,P)=\sum_{j=1}^n \omega_j I_j,
\end{align*}
where
\begin{align*}
I_j=\frac{1}{2}(Q_j^2+P_j^2).
\end{align*}
Because the change of coordinates is linear and symplectic, it preserves the standard Poisson bracket and changes neither the vanishing order of the Taylor remainder nor the local nature of the desired symplectomorphism. Thus it suffices to prove the theorem in these coordinates.
Let $\mathcal{P}_r$ denote the finite-dimensional real [vector space](/page/Vector%20Space) of homogeneous real polynomials of phase-space degree $r$ in $(Q,P)\in\mathbb{R}^{2n}$. Since $H\in C^{N+1}(U)$, [Taylor's theorem](/theorems/827) at $0$ gives
\begin{align*}
H(Q,P)=H(0)+H_2(Q,P)+H_3(Q,P)+\cdots+H_N(Q,P)+\rho_{N+1}(Q,P),
\end{align*}
where $H_r\in\mathcal{P}_r$ for $3\leq r\leq N$, and $\rho_{N+1}$ vanishes to order $N+1$ at $0$.
[/step]
[step:Diagonalize the homological operator on complex monomials]
For each $j\in\{1,\dots,n\}$, introduce the complex coordinate $z_j:\mathbb{R}^{2n}\to\mathbb{C}$ by
\begin{align*}
z_j(Q,P)=\frac{1}{\sqrt{2}}(Q_j+iP_j).
\end{align*}
Also define $\bar z_j:\mathbb{R}^{2n}\to\mathbb{C}$ by
\begin{align*}
\bar z_j(Q,P)=\frac{1}{\sqrt{2}}(Q_j-iP_j).
\end{align*}
Then $I_j=z_j\bar z_j$. For multi-indices $\alpha,\beta\in\mathbb{N}_0^n$, define
\begin{align*}
z^\alpha\bar z^\beta=\prod_{j=1}^n z_j^{\alpha_j}\bar z_j^{\beta_j}.
\end{align*}
The real polynomial space $\mathcal{P}_r$ is obtained by imposing the real-valuedness condition on complex linear combinations of monomials $z^\alpha\bar z^\beta$ with $|\alpha|+|\beta|=r$.
Define the homological operator $\mathcal{L}_r:\mathcal{P}_r\to\mathcal{P}_r$ by
\begin{align*}
\mathcal{L}_r(G)=\{H_2,G\},
\end{align*}
where $G\in\mathcal{P}_r$ and $\{\cdot,\cdot\}$ is the standard Poisson bracket. Since
\begin{align*}
\{H_2,z_j\}=-i\omega_j z_j
\end{align*}
and
\begin{align*}
\{H_2,\bar z_j\}=i\omega_j \bar z_j,
\end{align*}
the Leibniz rule for the Poisson bracket gives
\begin{align*}
\{H_2,z^\alpha\bar z^\beta\}=i((\beta-\alpha)\cdot\omega)z^\alpha\bar z^\beta.
\end{align*}
[guided]
The reason for passing to $z_j$ and $\bar z_j$ is that the oscillator flow rotates each complex coordinate independently. For each $j\in\{1,\dots,n\}$, define $z_j:\mathbb{R}^{2n}\to\mathbb{C}$ by
\begin{align*}
z_j(Q,P)=\frac{1}{\sqrt{2}}(Q_j+iP_j).
\end{align*}
Also define $\bar z_j:\mathbb{R}^{2n}\to\mathbb{C}$ by
\begin{align*}
\bar z_j(Q,P)=\frac{1}{\sqrt{2}}(Q_j-iP_j).
\end{align*}
Then
\begin{align*}
I_j=z_j\bar z_j.
\end{align*}
The standard Poisson bracket on $\mathbb{R}^{2n}$ is
\begin{align*}
\{F,G\}=\sum_{j=1}^n\left(\frac{\partial F}{\partial Q_j}\frac{\partial G}{\partial P_j}-\frac{\partial F}{\partial P_j}\frac{\partial G}{\partial Q_j}\right).
\end{align*}
Using $H_2=\sum_{j=1}^n\omega_j z_j\bar z_j$, direct differentiation gives
\begin{align*}
\{H_2,z_j\}=-i\omega_j z_j
\end{align*}
and
\begin{align*}
\{H_2,\bar z_j\}=i\omega_j\bar z_j.
\end{align*}
Now let $\alpha,\beta\in\mathbb{N}_0^n$ and define
\begin{align*}
z^\alpha\bar z^\beta=\prod_{j=1}^n z_j^{\alpha_j}\bar z_j^{\beta_j}.
\end{align*}
Applying the Leibniz rule to each factor gives
\begin{align*}
\{H_2,z^\alpha\bar z^\beta\}=i((\beta-\alpha)\cdot\omega)z^\alpha\bar z^\beta.
\end{align*}
Thus the homological operator $G\mapsto\{H_2,G\}$ is diagonal in the complex monomial basis. This is the central algebraic point: solving the homological equation reduces to dividing each nonresonant monomial coefficient by the scalar $i((\beta-\alpha)\cdot\omega)$.
[/guided]
[/step]
[step:Identify the resonant monomials with action polynomials]
Let $r\leq N$ and let $z^\alpha\bar z^\beta$ be a monomial of phase-space degree $r$, so $|\alpha|+|\beta|=r$. If $\alpha\neq\beta$, then $k=\beta-\alpha$ is a nonzero element of $\mathbb{Z}^n$ and
\begin{align*}
|k_1|+\cdots+|k_n|\leq |\alpha|+|\beta|=r\leq N.
\end{align*}
The nonresonance hypothesis gives $k\cdot\omega\neq 0$, so this monomial is not in the kernel of $\mathcal{L}_r$.
Therefore the kernel of $\mathcal{L}_r$ consists exactly of complex linear combinations of monomials with $\alpha=\beta$. Such monomials have the form
\begin{align*}
z^\alpha\bar z^\alpha=I_1^{\alpha_1}\cdots I_n^{\alpha_n}.
\end{align*}
Consequently, if $r$ is odd, $\ker\mathcal{L}_r=\{0\}$; if $r=2m$ is even, $\ker\mathcal{L}_{2m}$ is precisely the space of homogeneous polynomials of degree $m$ in the actions $I_1,\dots,I_n$.
[/step]
[step:Solve the homological equation in each degree]
For each $r\in\{3,\dots,N\}$, decompose any polynomial $F_r\in\mathcal{P}_r$ uniquely as
\begin{align*}
F_r=Z_r^{\mathrm{ph}}+F_r^{\mathrm{nr}},
\end{align*}
where $Z_r^{\mathrm{ph}}\in\ker\mathcal{L}_r$ and $F_r^{\mathrm{nr}}$ is the sum of all nonresonant monomials. If $r$ is odd, $Z_r^{\mathrm{ph}}=0$. If $r=2m$, then $Z_{2m}^{\mathrm{ph}}=Z_m(I)$ for a homogeneous polynomial $Z_m:\mathbb{R}^n\to\mathbb{R}$ of degree $m$.
Since $\mathcal{L}_r$ is invertible on the nonresonant summand, there exists a real homogeneous polynomial $G_r\in\mathcal{P}_r$ satisfying
\begin{align*}
\{H_2,G_r\}=-F_r^{\mathrm{nr}}.
\end{align*}
Indeed, if
\begin{align*}
F_r^{\mathrm{nr}}=\sum_{\{(\alpha,\beta): |\alpha|+|\beta|=r,\ \alpha\neq\beta\}} c_{\alpha\beta}z^\alpha\bar z^\beta,
\end{align*}
then one may take
\begin{align*}
G_r=\sum_{\{(\alpha,\beta): |\alpha|+|\beta|=r,\ \alpha\neq\beta\}} \frac{-c_{\alpha\beta}}{i((\beta-\alpha)\cdot\omega)}z^\alpha\bar z^\beta,
\end{align*}
and the real-valuedness of $G_r$ follows from the conjugate symmetry of the coefficients of the real polynomial $F_r^{\mathrm{nr}}$.
[/step]
[step:Remove nonresonant terms by successive Hamiltonian flows]
We prove by induction on $r\in\{3,\dots,N\}$ that there is a local symplectomorphism $\Phi_r$ fixing $0$ and tangent to the identity such that
\begin{align*}
H\circ\Phi_r=H(0)+H_2+Z_3^{\mathrm{ph}}+\cdots+Z_r^{\mathrm{ph}}+F_{r+1}+F_{r+2}+\cdots+F_N+\rho_{N+1,r},
\end{align*}
where $F_s\in\mathcal{P}_s$ for $s>r$ and $\rho_{N+1,r}$ vanishes to order $N+1$ at $0$.
For the initial stage $r=2$, take $\Phi_2=\operatorname{id}_{\mathbb{R}^{2n}}$. Suppose the assertion holds through degree $r-1$, where $3\leq r\leq N$. Decompose the current degree-$r$ term as
\begin{align*}
F_r=Z_r^{\mathrm{ph}}+F_r^{\mathrm{nr}},
\end{align*}
and choose $G_r\in\mathcal{P}_r$ with $\{H_2,G_r\}=-F_r^{\mathrm{nr}}$.
Let $\varphi_r^t$ denote the local Hamiltonian flow generated by $G_r$. Since $G_r$ is a smooth polynomial and its Hamiltonian vector field vanishes to order $r-1\geq2$ at $0$, the time-one map $\varphi_r^1$ is defined on a sufficiently small neighbourhood of $0$, is symplectic, fixes $0$, and is tangent to the identity. Taylor expansion of $H\circ\Phi_{r-1}\circ\varphi_r^1$ gives
\begin{align*}
H\circ\Phi_{r-1}\circ\varphi_r^1=H\circ\Phi_{r-1}+\{H\circ\Phi_{r-1},G_r\}+\text{terms of phase-space degree at least } r+1.
\end{align*}
Only $H_2$ contributes to degree $r$ in the Poisson bracket, because $\{F_s,G_r\}$ has phase-space degree $s+r-2\geq r+1$ whenever $s\geq3$. Hence the new degree-$r$ term is
\begin{align*}
F_r+\{H_2,G_r\}=Z_r^{\mathrm{ph}}.
\end{align*}
Setting $\Phi_r=\Phi_{r-1}\circ\varphi_r^1$ proves the induction step.
[/step]
[step:Collect the action terms and absorb the higher-order remainder]
After completing the induction at $r=N$, define
\begin{align*}
\Phi=\Phi_N.
\end{align*}
This map is a finite composition of local symplectomorphisms, hence is a local symplectomorphism on some neighbourhood $W$ of $0$ onto a neighbourhood $U_0\subset U$ of $0$. It satisfies $\Phi(0)=0$ and $D\Phi_0=\operatorname{id}_{\mathbb{R}^{2n}}$ because each factor has these properties.
The transformed Hamiltonian has the form
\begin{align*}
H\circ\Phi=H(0)+H_2+Z_3^{\mathrm{ph}}+\cdots+Z_N^{\mathrm{ph}}+R_{N+1}.
\end{align*}
By the kernel identification above, $Z_r^{\mathrm{ph}}=0$ for odd $r$, while for $r=2m$ one has
\begin{align*}
Z_{2m}^{\mathrm{ph}}=Z_m(I)
\end{align*}
with $Z_m$ homogeneous of degree $m$ in the actions. Since $2m\leq N$ is equivalent to $m\leq \lfloor N/2\rfloor$, this becomes
\begin{align*}
H\circ\Phi(Q,P)=H(0)+\sum_{j=1}^n\omega_j I_j+Z_2(I)+Z_3(I)+\cdots+Z_{\lfloor N/2\rfloor}(I)+R_{N+1}(Q,P).
\end{align*}
The remainder $R_{N+1}$ vanishes to order $N+1$ at $0$ because the construction has eliminated or retained explicitly every Taylor term of phase-space degree at most $N$. The Taylor coefficients of $\Phi$ through degree $N$ are determined recursively by the uniquely specified nonresonant solutions $G_r$, with the usual freedom of composing with action-dependent resonant flows, which does not change the displayed normal form through degree $N$. This proves the theorem.
[/step]