[proofplan]
We first localize the integral to a small neighbourhood of $K=\operatorname{supp}a\cap C$; the complementary compact region has no critical points and therefore contributes a rapidly decreasing remainder by nonstationary phase. Near each clean component $C_\ell$, we use Euclidean tubular coordinates, writing points as a tangential parameter $y \in C_\ell$ plus a normal variable $v$. The clean condition says exactly that the Hessian in the normal variable is nondegenerate at $v=0$, so stationary phase with parameters applies uniformly in $y$. Integrating the resulting normal stationary phase expansion over $C_\ell$ gives the full expansion, and the leading coefficient follows from the normal-coordinate Jacobian being $1$ on the zero section.
[/proofplan]
[step:Localize the integral near the clean critical set]
Define the oscillatory integral
\begin{align*}
I(\lambda) := \int_U e^{i\lambda\phi(x)}a(x)\,d\mathcal L^N(x)
\end{align*}
for $\lambda > 0$. Since $K$ is compact and $V$ is an open neighbourhood of $K$, choose a function $\chi \in C_c^\infty(V;\mathbb R)$ such that $\chi(x)=1$ on an open neighbourhood of $K$. Then
\begin{align*}
I(\lambda)
=
\int_U e^{i\lambda\phi(x)}\chi(x)a(x)\,d\mathcal L^N(x)
+
\int_U e^{i\lambda\phi(x)}(1-\chi(x))a(x)\,d\mathcal L^N(x).
\end{align*}
The support of $(1-\chi)a$ is compact and disjoint from $C$, hence $|\nabla\phi|$ has a positive lower bound on that support. By the nonstationary phase estimate, cited here as a standard result not yet resolved to a wiki theorem because the theorem search service was unavailable, for every $R \in \mathbb N$ there is a constant $B_R>0$ such that
\begin{align*}
\left|
\int_U e^{i\lambda\phi(x)}(1-\chi(x))a(x)\,d\mathcal L^N(x)
\right|
\le B_R\lambda^{-R}
\end{align*}
for all $\lambda \ge 1$. Thus this term is rapidly decreasing and does not affect any coefficient in the asserted asymptotic expansion.
[guided]
The point of this first reduction is to remove every part of the amplitude where the phase has no critical point. Define
\begin{align*}
I(\lambda) := \int_U e^{i\lambda\phi(x)}a(x)\,d\mathcal L^N(x)
\end{align*}
for $\lambda>0$. The only possible stationary contributions can come from
\begin{align*}
K := \operatorname{supp}a \cap C,
\end{align*}
where $C=\{x\in U:\nabla\phi(x)=0\}$. Since $K$ is compact and $V$ is an open neighbourhood of $K$, we may choose $\chi \in C_c^\infty(V;\mathbb R)$ with $\chi=1$ on an open neighbourhood of $K$. Splitting the amplitude gives
\begin{align*}
I(\lambda)=\int_U e^{i\lambda\phi(x)}\chi(x)a(x)\,d\mathcal L^N(x)+\int_U e^{i\lambda\phi(x)}(1-\chi(x))a(x)\,d\mathcal L^N(x).
\end{align*}
We now verify why the second integral is negligible. The function $(1-\chi)a$ has compact support, and this support is disjoint from $C$ because $\chi=1$ near $\operatorname{supp}a\cap C$. Therefore $\nabla\phi(x)\neq 0$ at every point of $\operatorname{supp}((1-\chi)a)$. Since this support is compact and $x\mapsto |\nabla\phi(x)|$ is continuous, there is a number $m>0$ such that
\begin{align*}
|\nabla\phi(x)| \ge m
\end{align*}
for all $x\in \operatorname{supp}((1-\chi)a)$. This is precisely the hypothesis needed for the nonstationary phase estimate. Citing that standard estimate as a result not yet resolved to a wiki theorem because theorem search was unavailable, for every $R\in\mathbb N$ there exists $B_R>0$ such that
\begin{align*}
\left|
\int_U e^{i\lambda\phi(x)}(1-\chi(x))a(x)\,d\mathcal L^N(x)
\right|
\le B_R\lambda^{-R}
\end{align*}
for all $\lambda\ge1$. Hence the second integral is smaller than every negative power of $\lambda$, so all asymptotic coefficients come from the localized amplitude $\chi a$.
[/guided]
[/step]
[step:Isolate the finitely many clean components near the amplitude]
Let $C_1,\dots,C_L$ be the connected components of $C\cap V$ that meet $K$. For each $\ell$, define the compact set
\begin{align*}
K_\ell := \operatorname{supp}a \cap C_\ell.
\end{align*}
We first justify that each $K_\ell$ has a neighbourhood in which no other critical component occurs. Fix $p\in K_\ell$. Since $\ker D^2\phi_p=T_pC_\ell$, the differential of the map $\nabla\phi:U\to\mathbb R^N$ has rank $k_\ell$ transverse to $T_pC_\ell$. By the constant-rank form of the [implicit function theorem](/theorems/52) applied to $\nabla\phi$ at $p$, after shrinking to an open neighbourhood $W_p\subset V$ of $p$, the critical set $C\cap W_p$ is exactly $C_\ell\cap W_p$. Since $K_\ell$ is compact, finitely many such sets cover $K_\ell$; let $W_\ell\subset V$ be their union, then $W_\ell\cap C=W_\ell\cap C_\ell$.
Because the compact sets $K_1,\dots,K_L$ are pairwise disjoint and finite in number, shrink the neighbourhoods $W_\ell$ so that they are pairwise disjoint and still satisfy $W_\ell\cap C=W_\ell\cap C_\ell$. Choose functions $\rho_1,\dots,\rho_L\in C_c^\infty(V;\mathbb R)$ such that $\operatorname{supp}\rho_\ell\subset W_\ell$ and $\rho_\ell=1$ on an open neighbourhood of $K_\ell$. Shrink the original cutoff $\chi$ if necessary so that
\begin{align*}
\operatorname{supp}(\chi a)\subset \bigcup_{\ell=1}^L \{x\in W_\ell : \rho_\ell(x)=1\}.
\end{align*}
Then $\sum_{\ell=1}^L\rho_\ell(x)=1$ on a neighbourhood of $\operatorname{supp}(\chi a)$, and $\operatorname{supp}\rho_\ell$ meets $C\cap V$ only in $C_\ell$.
Define the smooth compactly supported function $a_\ell:U\to\mathbb C$ by
\begin{align*}
a_\ell(x) := \rho_\ell(x)\chi(x)a(x).
\end{align*}
Then, using the rapidly decreasing remainder from the localization step,
\begin{align*}
I(\lambda)=\sum_{\ell=1}^L\int_U e^{i\lambda\phi(x)}a_\ell(x)\,d\mathcal L^N(x)+O(\lambda^{-R})
\end{align*}
for every $R\in\mathbb N$. It remains to analyze one fixed component $C_\ell$.
For $x\in C_\ell$ and $w\in T_xC_\ell$, choose a smooth curve $\gamma:(-\varepsilon,\varepsilon)\to C_\ell$ with $\gamma(0)=x$ and $\gamma'(0)=w$. Since $\nabla\phi=0$ on $C_\ell$,
\begin{align*}
\frac{d}{dt}\phi(\gamma(t))\bigg|_{t=0}=\nabla\phi(x)\cdot w=0.
\end{align*}
Thus $\phi|_{C_\ell}$ has zero differential. Since $C_\ell$ is connected, $\phi$ is constant on $C_\ell$; denote this value by $\phi_\ell$.
[/step]
[step:Write the localized integral in Euclidean normal coordinates]
Fix $\ell\in\{1,\dots,L\}$ and write $S:=C_\ell$, $e:=e_\ell$, and $k:=k_\ell$. Let $NS\to S$ be the Euclidean normal bundle, whose fiber at $y\in S$ is
\begin{align*}
N_yS := (T_yS)^\perp \subset \mathbb R^N.
\end{align*}
By the [tubular neighbourhood theorem](/theorems/2276) for embedded submanifolds of Euclidean space, cited here as a standard result not yet resolved to a wiki theorem because theorem search was unavailable, after restricting to a relatively compact open subset of $S$ containing $\operatorname{supp}a_\ell\cap S$, there is an open neighbourhood $\Omega_\ell\subset NS$ of the zero section and a diffeomorphism
\begin{align*}
\Psi_\ell:\Omega_\ell &\to \Psi_\ell(\Omega_\ell)\subset U
\end{align*}
given by
\begin{align*}
\Psi_\ell(y,v) := y+v.
\end{align*}
Choose $\Omega_\ell$ so that $\operatorname{supp}a_\ell\subset \Psi_\ell(\Omega_\ell)$ and so that the projection of $\Psi_\ell^{-1}(\operatorname{supp}a_\ell)$ to $S$ is contained in a compact set $P_\ell\subset S$.
Let $d\mathcal L^k_y(v)$ denote [Lebesgue measure](/page/Lebesgue%20Measure) on the Euclidean [vector space](/page/Vector%20Space) $N_yS$, and let $J_\ell:\Omega_\ell\to(0,\infty)$ be the smooth density satisfying the change-of-variables formula
\begin{align*}
\int_{\Psi_\ell(\Omega_\ell)} f(x)\,d\mathcal L^N(x)=\int_S\int_{\Omega_{\ell,y}} f(\Psi_\ell(y,v))J_\ell(y,v)\,d\mathcal L^k_y(v)\,d\mathcal H^e(y)
\end{align*}
for every $f\in C_c(\Psi_\ell(\Omega_\ell);\mathbb C)$, where $\Omega_{\ell,y}:=\Omega_\ell\cap N_yS$. For $f(x)=e^{i\lambda\phi(x)}a_\ell(x)$ the absolute value is bounded by $|a_\ell(x)|$, which is compactly supported; hence the pulled-back integrand is absolutely integrable on the [measure space](/page/Measure%20Space) with measure $d\mathcal L^k_y(v)\,d\mathcal H^e(y)$. [Fubini's theorem](/theorems/2961) therefore justifies the iterated integral over $S$ and the normal fibers. Since $d\Psi_\ell$ maps tangent and normal directions orthogonally onto $\mathbb R^N$ at $v=0$, the density satisfies
\begin{align*}
J_\ell(y,0)=1
\end{align*}
for every $y\in S$.
Define the parameter-dependent amplitude
\begin{align*}
b_\ell:\Omega_\ell &\to \mathbb C
\end{align*}
by
\begin{align*}
b_\ell(y,v):=a_\ell(\Psi_\ell(y,v))J_\ell(y,v).
\end{align*}
Then the localized integral is
\begin{align*}
I_\ell(\lambda)
:=
\int_U e^{i\lambda\phi(x)}a_\ell(x)\,d\mathcal L^N(x)
=
\int_S\int_{\Omega_{\ell,y}} e^{i\lambda\phi(\Psi_\ell(y,v))}b_\ell(y,v)\,d\mathcal L^k_y(v)\,d\mathcal H^e(y).
\end{align*}
[/step]
[step:Apply stationary phase in the normal variable uniformly over the component]
For $y\in S$, define the normal phase
\begin{align*}
\Phi_{\ell,y}:\Omega_{\ell,y} &\to \mathbb R
\end{align*}
by
\begin{align*}
\Phi_{\ell,y}(v):=\phi(\Psi_\ell(y,v)).
\end{align*}
Because $y\in C$, the first derivative of $\Phi_{\ell,y}$ at $v=0$ vanishes. Its Hessian at $v=0$ is the [bilinear form](/page/Bilinear%20Form)
\begin{align*}
D^2\Phi_{\ell,y}(0)(v,w)=D^2\phi_y(v,w)
\end{align*}
for $v,w\in N_yS$. The clean condition $\ker D^2\phi_y=T_yS$ implies that this Hessian is nondegenerate on $N_yS$. Hence the signature of this normal Hessian is locally constant in $y$; since $S$ is connected, it is constant, and we denote it by $\sigma_\ell$.
For $y\in P_\ell$, define the normal-gradient map $G_y:\Omega_{\ell,y}\to N_yS$ by requiring
\begin{align*}
G_y(v)\cdot w=D_v\Phi_{\ell,y}(v)(w)
\end{align*}
for every $w\in N_yS$, where the dot product is the Euclidean [inner product](/page/Inner%20Product) on $N_yS$. At $v=0$, the derivative $D_vG_y(0):N_yS\to N_yS$ is represented by the nondegenerate bilinear form $H_\ell(y)$. By the [inverse function theorem](/theorems/51) with parameters and compactness of $P_\ell$, after shrinking $\Omega_\ell$ around the zero section and shrinking $\operatorname{supp}a_\ell$ inside that neighbourhood, $v=0$ is the only critical point of $\Phi_{\ell,y}$ on $\Omega_{\ell,y}\cap\operatorname{supp}_v b_\ell$ for every $y\in P_\ell$. Here $\operatorname{supp}_v b_\ell$ denotes the fiberwise support of $v\mapsto b_\ell(y,v)$.
Cover $P_\ell$ by finitely many coordinate patches on $S$ and finitely many smooth orthonormal frames of the normal bundle over those patches. Choose a subordinate smooth [partition of unity](/page/Partition%20of%20Unity) on $P_\ell$. In each patch, the [stationary phase theorem](/theorems/8198) with parameters applies in the fixed Euclidean variable representing $v$, with $y$ ranging over a compact parameter set; its remainder constants are uniform on that compact set. Summing over the finite partition and integrating the uniform remainders over $P_\ell$ gives smooth compactly supported coefficient densities
\begin{align*}
L_{j,\ell}:S\to\mathbb C
\end{align*}
such that, for every $M\in\mathbb N$,
\begin{align*}
I_\ell(\lambda)
=
e^{i\lambda\phi_\ell}\lambda^{-k/2}
\sum_{j=0}^{M-1}\lambda^{-j}
\int_S L_{j,\ell}(y)\,d\mathcal H^e(y)
+
O(\lambda^{-k/2-M}).
\end{align*}
The leading coefficient density is
\begin{align*}
L_{0,\ell}(y)
=
(2\pi)^{k/2}e^{i\pi\sigma_\ell/4}
b_\ell(y,0)|\det H_\ell(y)|^{-1/2}.
\end{align*}
[guided]
We now focus on the normal oscillation while treating the point $y\in S$ as a parameter. For each $y\in S$, define
\begin{align*}
\Phi_{\ell,y}:\Omega_{\ell,y} &\to \mathbb R
\end{align*}
by
\begin{align*}
\Phi_{\ell,y}(v):=\phi(\Psi_\ell(y,v)).
\end{align*}
The stationary point in the normal variable is $v=0$. Indeed, if $w\in N_yS$, then
\begin{align*}
D\Phi_{\ell,y}(0)(w)=\nabla\phi(y)\cdot w=0,
\end{align*}
because $y\in C$ and therefore $\nabla\phi(y)=0$.
Next we compute the Hessian in the normal variable. Since $\Psi_\ell(y,v)=y+v$ is affine in $v$, the second derivative in the $v$ directions has no extra coordinate term. Thus, for $v,w\in N_yS$,
\begin{align*}
D^2\Phi_{\ell,y}(0)(v,w)=D^2\phi_y(v,w).
\end{align*}
The clean critical set hypothesis says
\begin{align*}
\ker D^2\phi_y=T_yS.
\end{align*}
Therefore, after passing to the orthogonal normal space $N_yS=(T_yS)^\perp$, the induced Hessian has no kernel. This is exactly the nondegeneracy condition required for ordinary stationary phase in the normal variable.
The signature must also be stable in $y$. The map $y\mapsto H_\ell(y)$ is a smooth family of nondegenerate symmetric bilinear forms on the normal bundle. Eigenvalues of a symmetric form vary continuously in local orthonormal frames, and no eigenvalue can cross $0$ because the forms remain nondegenerate. Hence the number of positive eigenvalues and the number of negative eigenvalues are locally constant. Since $S=C_\ell$ is connected, the signature is constant on $S$; denote this constant by $\sigma_\ell$.
We must also ensure that $v=0$ is the only normal stationary point seen by the amplitude, uniformly in $y$. Let $P_\ell\subset S$ be the compact projection of $\Psi_\ell^{-1}(\operatorname{supp}a_\ell)$ to the base. For $y\in P_\ell$, define $G_y:\Omega_{\ell,y}\to N_yS$ by
\begin{align*}
G_y(v)\cdot w=D_v\Phi_{\ell,y}(v)(w)
\end{align*}
for every $w\in N_yS$. The derivative $D_vG_y(0)$ is the nondegenerate normal Hessian $H_\ell(y)$. The inverse function theorem with parameters gives a neighbourhood of the zero section in which $G_y(v)=0$ has the unique solution $v=0$; compactness of $P_\ell$ lets us choose this neighbourhood uniformly in $y$. We shrink the support of $a_\ell$ inside that neighbourhood, which does not change $a_\ell$ near $K_\ell$ and changes only a nonstationary contribution away from $C_\ell$.
Now cover $P_\ell$ by finitely many coordinate patches on $S$ and choose finitely many smooth orthonormal frames for the normal bundle over those patches. A subordinate partition of unity reduces the integral to finitely many parameter-dependent stationary phase problems in a fixed Euclidean normal variable. On each patch the parameter set is compact, so the stationary phase remainder is uniform in $y$; after integration over $P_\ell$, the same order $O(\lambda^{-k/2-M})$ remains. Therefore there are smooth compactly supported functions
\begin{align*}
L_{j,\ell}:S\to\mathbb C
\end{align*}
such that for every $M\in\mathbb N$,
\begin{align*}
I_\ell(\lambda)
=
e^{i\lambda\phi_\ell}\lambda^{-k/2}
\sum_{j=0}^{M-1}\lambda^{-j}
\int_S L_{j,\ell}(y)\,d\mathcal H^e(y)
+
O(\lambda^{-k/2-M}).
\end{align*}
The leading term in the stationary phase formula in dimension $k$ is the value of the amplitude at the critical point multiplied by the Gaussian factor. Hence
\begin{align*}
L_{0,\ell}(y)
=
(2\pi)^{k/2}e^{i\pi\sigma_\ell/4}
b_\ell(y,0)|\det H_\ell(y)|^{-1/2}.
\end{align*}
This is the place where the clean hypothesis is used quantitatively: it turns the transverse Hessian into the nondegenerate quadratic form whose determinant and signature enter the leading coefficient.
[/guided]
[/step]
[step:Identify the leading coefficient and assemble the expansion]
Since $\Psi_\ell(y,0)=y$ and $J_\ell(y,0)=1$, the definition of $b_\ell$ gives
\begin{align*}
b_\ell(y,0)=\rho_\ell(y)\chi(y)a(y)
\end{align*}
for $y\in S$. By construction, $\rho_\ell\chi=1$ on a neighbourhood of $K_\ell=\operatorname{supp}a\cap S$, while $a(y)=0$ for $y\in S\setminus K_\ell$. Hence
\begin{align*}
b_\ell(y,0)=a(y)
\end{align*}
for every $y\in S$. Therefore
\begin{align*}
\int_S L_{0,\ell}(y)\,d\mathcal H^e(y)
=
(2\pi)^{k/2}e^{i\pi\sigma_\ell/4}
\int_S a(y)|\det H_\ell(y)|^{-1/2}\,d\mathcal H^e(y).
\end{align*}
Returning to the original notation $S=C_\ell$, $e=e_\ell$, and $k=k_\ell$, define
\begin{align*}
A_{j,\ell}:=\int_{C_\ell} L_{j,\ell}(y)\,d\mathcal H^{e_\ell}(y)
\end{align*}
for $j\ge1$, and define $A_{0,\ell}$ by the displayed leading formula. Summing over the finitely many components and adding the rapidly decreasing nonstationary remainder gives, for every $M\in\mathbb N$,
\begin{align*}
I(\lambda)
-
\sum_{\ell=1}^{L}
e^{i\lambda\phi_\ell}\lambda^{-k_\ell/2}
\sum_{j=0}^{M-1}\lambda^{-j}A_{j,\ell}
=
\sum_{\ell=1}^{L}O(\lambda^{-k_\ell/2-M}).
\end{align*}
This is precisely the asserted asymptotic expansion, and the displayed expression for $A_{0,\ell}$ is the required explicit leading coefficient.
[/step]