[proofplan]
We solve the first-order homological equation mode by mode. The nonresonance lower bound allows division by $k \cdot \omega(I)$ for the selected modes, producing a smooth real-valued generating Hamiltonian $S_1$. The time-$\varepsilon$ Hamiltonian flow of $S_1$ is canonical and near the identity, and Taylor expansion of $H \circ \Phi_\varepsilon$ shows that the selected Fourier modes cancel at order $\varepsilon$.
[/proofplan]
[step:Choose a neighbourhood where the nonresonance denominators stay bounded away from zero]
Since $K_0 \subset D$ is compact and each map
\begin{align*}
I \mapsto k \cdot \omega(I)
\end{align*}
is continuous, the hypothesis gives
\begin{align*}
|k \cdot \omega(I)| \geq \gamma
\end{align*}
for every $I \in K_0$ and every $k \in A$. Because $A$ is finite, there exists an open neighbourhood $U \subset D$ of $K_0$ such that
\begin{align*}
|k \cdot \omega(I)| \geq \frac{\gamma}{2}
\end{align*}
for every $I \in U$ and every $k \in A$.
This ensures that all denominators used below are smooth and uniformly bounded on $U$.
[/step]
[step:Define the generating Hamiltonian by solving the homological equation mode by mode]
Define
\begin{align*}
S_1: U \times \mathbb{T}^n &\to \mathbb{R}
\end{align*}
by
\begin{align*}
S_1(I,\theta)
=
\sum_{k \in A}
-\frac{H_{1,k}(I)}{i\, k \cdot \omega(I)}
e^{i k \cdot \theta}.
\end{align*}
The denominator is nonzero on $U$ by the previous step, so each coefficient is smooth on $U$. Since $A$ is finite, $S_1$ is smooth.
Because $H_1$ is real-valued, its Fourier coefficients satisfy
\begin{align*}
H_{1,-k}(I)=\overline{H_{1,k}(I)}.
\end{align*}
Since $A$ is symmetric and
\begin{align*}
\left(-\frac{H_{1,-k}(I)}{i\,(-k)\cdot\omega(I)}\right)e^{-i k\cdot\theta}
=
\overline{
\left(-\frac{H_{1,k}(I)}{i\,k\cdot\omega(I)}\right)e^{i k\cdot\theta}
},
\end{align*}
the sum defining $S_1$ is real-valued.
[guided]
The obstruction to removing a Fourier mode is the factor $k \cdot \omega(I)$. To remove the mode with index $k$, we want a generating Hamiltonian whose Poisson bracket with $H_0$ reproduces precisely that mode with the correct sign. This leads to division by $i\,k \cdot \omega(I)$.
We define the map
\begin{align*}
S_1: U \times \mathbb{T}^n &\to \mathbb{R}
\end{align*}
by
\begin{align*}
S_1(I,\theta)
=
\sum_{k \in A}
-\frac{H_{1,k}(I)}{i\, k \cdot \omega(I)}
e^{i k \cdot \theta}.
\end{align*}
This is legitimate because the previous step gives
\begin{align*}
|k \cdot \omega(I)| \geq \frac{\gamma}{2}
\end{align*}
for every $I \in U$ and every $k \in A$. Thus no denominator vanishes, and each coefficient
\begin{align*}
I \mapsto -\frac{H_{1,k}(I)}{i\, k \cdot \omega(I)}
\end{align*}
is a smooth complex-valued function on $U$. Since $A$ is finite, the Fourier sum is finite, so $S_1$ is smooth.
We also need $S_1$ to be real-valued, because it will generate a real Hamiltonian flow. The function $H_1$ is real-valued, so its Fourier coefficients satisfy
\begin{align*}
H_{1,-k}(I)=\overline{H_{1,k}(I)}.
\end{align*}
For the paired mode $-k$, the corresponding summand is
\begin{align*}
\left(-\frac{H_{1,-k}(I)}{i\,(-k)\cdot\omega(I)}\right)e^{-i k\cdot\theta}
=
\overline{
\left(-\frac{H_{1,k}(I)}{i\,k\cdot\omega(I)}\right)e^{i k\cdot\theta}
}.
\end{align*}
Because $A$ is symmetric, the terms occur in conjugate pairs. Therefore their sum is real-valued.
[/guided]
[/step]
[step:Verify that the Poisson bracket cancels exactly the selected modes]
Use the canonical action-angle Poisson bracket
\begin{align*}
\{F,G\}
=
\nabla_I F \cdot \nabla_\theta G
-
\nabla_\theta F \cdot \nabla_I G.
\end{align*}
Since $H_0$ depends only on $I$, we have $\nabla_\theta H_0=0$ and $\nabla_I H_0=\omega$. Hence
\begin{align*}
\{H_0,S_1\}
=
\omega(I)\cdot\nabla_\theta S_1.
\end{align*}
Differentiating the finite Fourier sum in $\theta$ gives
\begin{align*}
\nabla_\theta S_1(I,\theta)
=
\sum_{k \in A}
-\frac{H_{1,k}(I)}{i\, k \cdot \omega(I)}
i k e^{i k \cdot \theta}.
\end{align*}
Taking the dot product with $\omega(I)$ yields
\begin{align*}
\{H_0,S_1\}(I,\theta)
=
-\sum_{k \in A} H_{1,k}(I)e^{i k \cdot \theta}.
\end{align*}
Thus $S_1$ solves the homological equation for exactly the selected modes.
[/step]
[step:Expand the transformed Hamiltonian to first order]
Shrink $U$ if necessary and choose a compact neighbourhood $K_1\subset U$ of $K_0$. Let $\Phi_t:K_1 \times \mathbb{T}^n \to U \times \mathbb{T}^n$ denote the Hamiltonian flow generated by $S_1$, defined for $|t|$ sufficiently small. Since Hamiltonian flows preserve the canonical symplectic form, $\Phi_t$ is canonical. Since $\Phi_0=\operatorname{id}$ and $S_1$ is smooth, $\Phi_\varepsilon$ is near the identity as $\varepsilon \to 0$.
We use the Hamiltonian-flow convention for which differentiation along the flow gives
\begin{align*}
\frac{d}{dt}F\circ\Phi_t
=
\{F,S_1\}\circ\Phi_t.
\end{align*}
Applying Taylor expansion at $t=0$ to $F=H_0+\varepsilon H_1$ gives, uniformly on $K_0\times\mathbb T^n$,
\begin{align*}
H\circ\Phi_\varepsilon
=
H
+
\varepsilon\{H,S_1\}
+
O(\varepsilon^2).
\end{align*}
Since $H=H_0+\varepsilon H_1$, this becomes
\begin{align*}
H\circ\Phi_\varepsilon
=
H_0
+
\varepsilon H_1
+
\varepsilon\{H_0,S_1\}
+
O(\varepsilon^2),
\end{align*}
because the contribution $\varepsilon^2\{H_1,S_1\}$ is second order. The $O(\varepsilon^2)$ remainder is uniform on $K_0\times\mathbb T^n$ because $K_1\times\mathbb T^n$ is compact, the flow remains in $U\times\mathbb T^n$ for sufficiently small $|\varepsilon|$, and the required derivatives of $H_0$, $H_1$, and $S_1$ are bounded there.
[/step]
[step:Separate the Fourier modes and identify the averaged normal form]
Decompose $H_1$ as
\begin{align*}
H_1(I,\theta)
=
H_{1,0}(I)
+
\sum_{k \in A}H_{1,k}(I)e^{i k \cdot \theta}
+
\sum_{k \notin A \cup \{0\}}H_{1,k}(I)e^{i k \cdot \theta}.
\end{align*}
Using the homological equation from the previous step,
\begin{align*}
\{H_0,S_1\}
=
-
\sum_{k \in A}H_{1,k}(I)e^{i k \cdot \theta}.
\end{align*}
Therefore
\begin{align*}
H\circ\Phi_\varepsilon
=
H_0
+
\varepsilon
\left(
H_{1,0}(I)
+
\sum_{k \notin A \cup \{0\}}H_{1,k}(I)e^{i k \cdot \theta}
\right)
+
O(\varepsilon^2).
\end{align*}
The selected nonresonant angle-dependent modes have disappeared from the first-order part. The zeroth Fourier mode $H_{1,0}$ is the angle average, while modes outside $A$ and modes for which the denominator $k\cdot\omega(I)$ is not controlled remain. This proves the theorem.
[/step]