[proofplan]
The proof constructs a first-order differential operator whose action on the oscillatory exponential is exactly the identity. Since the phase has no critical point on the support of the amplitude, the coefficients of this operator are smooth and uniformly bounded near the support. Repeatedly moving the operator from the exponential onto the amplitude by compact-support integration by parts produces one factor of $h$ at each step. After $N$ integrations by parts, the remaining integral is bounded by the $L^1$ norm of a compactly supported amplitude controlled by finitely many derivatives of $a$ and $\phi$.
[/proofplan]
[step:Choose a neighbourhood where the phase has no critical points]
Since $K \subset U$ is compact and $U$ is open, choose an open set $V \subset \mathbb{R}^n$ such that $K \subset V$, $\overline{V} \subset U$, and $\overline{V}$ is compact.
By continuity of the map $x \mapsto |\nabla \phi(x)|$ and the lower bound on $K$, after replacing $V$ by a smaller open neighbourhood of $K$ if necessary, we may also assume
\begin{align*}
|\nabla \phi(x)| \geq \frac{c}{2}
\end{align*}
for every $x \in V$.
Define the smooth vector field $A: V \to \mathbb{R}^n$ by
\begin{align*}
A(x) = \frac{\nabla \phi(x)}{|\nabla \phi(x)|^2}.
\end{align*}
This is well-defined on $V$ because $|\nabla \phi|$ is bounded below there.
[/step]
[step:Build the differential operator that fixes the oscillatory exponential]
For each $h \in (0,h_0]$, define the first-order differential operator $L_h: C^\infty(V;\mathbb{C}) \to C^\infty(V;\mathbb{C})$ by
\begin{align*}
L_h f = \frac{h}{i} A \cdot \nabla f.
\end{align*}
For the smooth function $E_h: V \to \mathbb{C}$ given by
\begin{align*}
E_h(x) = e^{i\phi(x)/h},
\end{align*}
we have
\begin{align*}
\nabla E_h(x) = \frac{i}{h} e^{i\phi(x)/h} \nabla \phi(x).
\end{align*}
Therefore
\begin{align*}
L_h E_h(x) = \frac{h}{i}\frac{\nabla \phi(x)}{|\nabla \phi(x)|^2} \cdot \frac{i}{h}e^{i\phi(x)/h}\nabla \phi(x) = e^{i\phi(x)/h}
\end{align*}
for every $x \in V$.
[/step]
[step:Move the operator from the exponential onto the amplitude]
Define the formal adjoint operator $L_h^*: C_c^\infty(V;\mathbb{C}) \to C_c^\infty(V;\mathbb{C})$ by
\begin{align*}
L_h^* b = -\frac{h}{i}\operatorname{div}(bA)
\end{align*}
for $b \in C_c^\infty(V;\mathbb{C})$. Since $b$ has compact support in $V$, no boundary term appears when integrating by parts in $\mathbb{R}^n$, and for every $f \in C^\infty(V;\mathbb{C})$ and $b \in C_c^\infty(V;\mathbb{C})$,
\begin{align*}
\int_{\mathbb{R}^n} (L_h f)(x)b(x)\,d\mathcal{L}^n(x) = \int_{\mathbb{R}^n} f(x)(L_h^*b)(x)\,d\mathcal{L}^n(x).
\end{align*}
Here both integrals are over $\operatorname{supp} b \subset V$.
Since $\operatorname{supp} a \subset K \subset V$, the amplitude $a$ lies in $C_c^\infty(V;\mathbb{C})$. Using $L_h E_h = E_h$ and the adjoint identity with $f = E_h$ and $b = a$ gives
\begin{align*}
\int_{\mathbb{R}^n} E_h(x)a(x)\,d\mathcal{L}^n(x) = \int_{\mathbb{R}^n} E_h(x)(L_h^*a)(x)\,d\mathcal{L}^n(x).
\end{align*}
Iterating this identity $N$ times yields
\begin{align*}
\int_{\mathbb{R}^n} E_h(x)a(x)\,d\mathcal{L}^n(x) = \int_{\mathbb{R}^n} E_h(x)(L_h^*)^N a(x)\,d\mathcal{L}^n(x).
\end{align*}
[guided]
The purpose of $L_h$ is that it differentiates the exponential and gives back exactly the same exponential. Since
\begin{align*}
\nabla E_h(x) = \frac{i}{h}E_h(x)\nabla \phi(x),
\end{align*}
the coefficient $A = \nabla\phi/|\nabla\phi|^2$ is chosen so that $A \cdot \nabla\phi = 1$. Thus $L_h E_h = E_h$.
Now we want to replace the oscillatory integral by one where derivatives fall on the compactly supported amplitude. For $b \in C_c^\infty(V;\mathbb{C})$, define
\begin{align*}
L_h^*b = -\frac{h}{i}\operatorname{div}(bA).
\end{align*}
This is the correct adjoint because compact support removes boundary terms. Indeed,
\begin{align*}
\int_{\mathbb{R}^n} (L_h f)(x)b(x)\,d\mathcal{L}^n(x) = \frac{h}{i}\int_{\mathbb{R}^n} A(x)\cdot \nabla f(x)b(x)\,d\mathcal{L}^n(x).
\end{align*}
Expanding the dot product and integrating each compactly supported derivative term by parts gives
\begin{align*}
\frac{h}{i}\int_{\mathbb{R}^n} A(x)\cdot \nabla f(x)b(x)\,d\mathcal{L}^n(x) = -\frac{h}{i}\int_{\mathbb{R}^n} f(x)\operatorname{div}(bA)(x)\,d\mathcal{L}^n(x).
\end{align*}
Therefore
\begin{align*}
\int_{\mathbb{R}^n} (L_h f)(x)b(x)\,d\mathcal{L}^n(x) = \int_{\mathbb{R}^n} f(x)(L_h^*b)(x)\,d\mathcal{L}^n(x).
\end{align*}
Applying this with $f = E_h$ and $b = a$ gives
\begin{align*}
\int_{\mathbb{R}^n} E_h(x)a(x)\,d\mathcal{L}^n(x) = \int_{\mathbb{R}^n} (L_hE_h)(x)a(x)\,d\mathcal{L}^n(x) = \int_{\mathbb{R}^n} E_h(x)(L_h^*a)(x)\,d\mathcal{L}^n(x).
\end{align*}
The same argument can be repeated because $L_h^*a$ is still smooth and compactly supported in $V$. Repeating $N$ times gives
\begin{align*}
\int_{\mathbb{R}^n} E_h(x)a(x)\,d\mathcal{L}^n(x) = \int_{\mathbb{R}^n} E_h(x)(L_h^*)^N a(x)\,d\mathcal{L}^n(x).
\end{align*}
[/guided]
[/step]
[step:Estimate the iterated adjoint by finitely many seminorms]
Define the first-order differential operator $P: C_c^\infty(V;\mathbb{C}) \to C_c^\infty(V;\mathbb{C})$ by
\begin{align*}
Pb = -\frac{1}{i}\operatorname{div}(bA).
\end{align*}
Then $L_h^* = hP$, and hence
\begin{align*}
(L_h^*)^N a = h^N P^N a.
\end{align*}
The operator $P^N$ is a differential operator of order at most $N$ whose coefficients are finite sums of products of derivatives of $A$ of order at most $N$. Since $A = \nabla\phi/|\nabla\phi|^2$ and $|\nabla\phi| \geq c/2$ on $V$, these coefficients are bounded on $\overline{V}$ by constants depending only on $N$, $c$, and finitely many derivatives of $\phi$ on $\overline{V}$.
Therefore there exists a constant $M_N > 0$, depending only on $N$, $c$, $\overline{V}$, and finitely many derivatives of $\phi$ on $\overline{V}$, such that
\begin{align*}
|P^N a(x)| \leq M_N \sum_{|\alpha|\leq N}|D^\alpha a(x)|
\end{align*}
for every $x \in V$. Since $\operatorname{supp} P^N a \subset \operatorname{supp} a \subset K$, it follows that
\begin{align*}
\int_{\mathbb{R}^n}|P^N a(x)|\,d\mathcal{L}^n(x) \leq M_N \mathcal{L}^n(K)\sum_{|\alpha|\leq N}\sup_{x\in K}|D^\alpha a(x)|.
\end{align*}
If $\mathcal{L}^n(K)=0$, then $a=0$: otherwise some point $x_0 \in U$ would satisfy $a(x_0) \neq 0$, and continuity would give an open ball on which $a$ is nonzero, forcing $\operatorname{supp}a$ and hence $K$ to have positive $\mathcal{L}^n$-measure. Thus the displayed estimate is valid in all cases. Consequently, there is a finite constant $B_N>0$ such that
\begin{align*}
\|P^N a\|_{L^1(\mathbb{R}^n)} \leq B_N.
\end{align*}
[/step]
[step:Conclude the $h^N$ decay estimate]
Using $|E_h(x)| = 1$ for every $x \in V$, the iterated integration-by-parts identity gives
\begin{align*}
\left|\int_{\mathbb{R}^n} E_h(x)a(x)\,d\mathcal{L}^n(x)\right| \leq \int_{\mathbb{R}^n}|(L_h^*)^N a(x)|\,d\mathcal{L}^n(x).
\end{align*}
Since $(L_h^*)^N a = h^N P^N a$, we obtain
\begin{align*}
\left|\int_{\mathbb{R}^n} e^{i\phi(x)/h}a(x)\,d\mathcal{L}^n(x)\right| \leq h^N\|P^N a\|_{L^1(\mathbb{R}^n)} \leq B_N h^N.
\end{align*}
Taking $C_N = B_N$ proves the desired estimate for every $h \in (0,h_0]$. The dependence of $C_N$ is exactly through the compact neighbourhood $V$, the lower bound $c$, the number of integrations by parts $N$, and the finitely many derivatives of $a$ and $\phi$ used in the bound for $P^N a$.
[/step]
