[proofplan]
We first separate the integral into a neighbourhood of the unique critical point and its complement. The complement is rapidly decaying by the nonstationary phase estimate, because $\nabla \phi$ is bounded away from zero there. Near $y_0$, the Morse lemma changes variables so that the phase becomes a signed quadratic form; after absorbing the Jacobian into the amplitude, the expansion follows by Taylor expanding the transformed amplitude and evaluating oscillatory Gaussian moments. The leading constant comes from the product of one-dimensional Fresnel integrals, and the finite-order remainder is controlled after the scaling $z = \lambda^{-1/2}w$.
[/proofplan]
[step:Split the integral into a critical and noncritical contribution]
Choose an open neighbourhood $V \subset U$ of $y_0$ whose closure is compact in $U$ and such that $y_0$ is the only critical point of $\phi$ in $\overline{V} \cap \operatorname{supp} a$. Choose $V$ small enough that it is contained in the Morse-coordinate image constructed in the next step. Let $\chi \in C_c^\infty(V;[0,1])$ be a cutoff function satisfying $\chi = 1$ on some open neighbourhood $V_0 \subset V$ of $y_0$ and $\operatorname{supp}\chi \subset V$ compactly. Define the oscillatory integral map $I:[1,\infty)\to\mathbb{C}$ by
\begin{align*}
I(\lambda) := \int_U e^{i\lambda\phi(y)}a(y)\,d\mathcal{L}^N(y).
\end{align*}
We decompose $I(\lambda) = I_{\mathrm{near}}(\lambda) + I_{\mathrm{far}}(\lambda)$, where $I_{\mathrm{near}}:[1,\infty)\to\mathbb{C}$ and $I_{\mathrm{far}}:[1,\infty)\to\mathbb{C}$ are defined by
\begin{align*}
I_{\mathrm{near}}(\lambda) := \int_U e^{i\lambda\phi(y)}\chi(y)a(y)\,d\mathcal{L}^N(y)
\end{align*}
and
\begin{align*}
I_{\mathrm{far}}(\lambda) := \int_U e^{i\lambda\phi(y)}(1-\chi(y))a(y)\,d\mathcal{L}^N(y).
\end{align*}
The amplitude $(1-\chi)a$ is compactly supported in $U$ and vanishes near $y_0$. Since $y_0$ is the unique critical point of $\phi$ on $\operatorname{supp} a$, the compact set $\operatorname{supp}((1-\chi)a)$ contains no critical point of $\phi$. Hence there is a constant $c_0 > 0$ such that
\begin{align*}
|\nabla \phi(y)| \geq c_0
\end{align*}
for every $y \in \operatorname{supp}((1-\chi)a)$.
By the nonstationary phase estimate, applied to the smooth real phase $\phi$ and the compactly supported amplitude $(1-\chi)a$, for every integer $m \geq 0$ there is a constant $C_m > 0$ such that
\begin{align*}
|I_{\mathrm{far}}(\lambda)| \leq C_m\lambda^{-m}
\end{align*}
for all $\lambda \geq 1$. Thus $I_{\mathrm{far}}(\lambda)$ is rapidly decreasing and contributes no term to the asymptotic expansion.
[/step]
[step:Put the phase into Morse normal form near $y_0$]
By the Morse lemma, since $\nabla\phi(y_0)=0$ and $H=J(\nabla\phi)_{y_0}$ is invertible, there exist an open neighbourhood $W \subset \mathbb{R}^N$ of $0$ and a $C^\infty$ diffeomorphism $\kappa: W \to \kappa(W) \subset U$ such that $\kappa(0)=y_0$, $\operatorname{supp}\chi \subset \kappa(W)$, and
\begin{align*}
\phi(\kappa(z)) = \phi(y_0) + \frac{1}{2}\sum_{r=1}^N \epsilon_r z_r^2
\end{align*}
for every $z=(z_1,\dots,z_N)\in W$, where each $\epsilon_r \in \{1,-1\}$. The number of indices $r$ with $\epsilon_r=1$ equals the number of positive eigenvalues of $H$, and the number with $\epsilon_r=-1$ equals the number of negative eigenvalues of $H$.
Define the transformed amplitude $b \in C_c^\infty(W;\mathbb{C})$ by
\begin{align*}
b(z) := \chi(\kappa(z))a(\kappa(z))|\det J\kappa_z|.
\end{align*}
Here $J\kappa_z$ is the Jacobian matrix of $\kappa$ at $z$. Because $\operatorname{supp}\chi \subset \kappa(W)$ is compact and $\kappa^{-1}:\kappa(W)\to W$ is continuous, $\operatorname{supp}b \subset \kappa^{-1}(\operatorname{supp}\chi \cap \operatorname{supp}a)$ is a compact subset of $W$. The change-of-variables theorem gives
\begin{align*}
I_{\mathrm{near}}(\lambda)=e^{i\lambda\phi(y_0)}\int_W e^{i\lambda Q(z)/2}b(z)\,d\mathcal{L}^N(z),
\end{align*}
where $Q:\mathbb{R}^N\to\mathbb{R}$ is the quadratic form
\begin{align*}
Q(z):=\sum_{r=1}^N \epsilon_r z_r^2.
\end{align*}
Moreover the Morse coordinate map may be chosen so that
\begin{align*}
|\det J\kappa_0|=|\det H|^{-1/2}.
\end{align*}
Consequently
\begin{align*}
b(0)=a(y_0)|\det H|^{-1/2}.
\end{align*}
[guided]
The goal of this step is to replace the original phase by the only part of it that matters at leading order: its nondegenerate quadratic behaviour at the critical point. Since $\nabla\phi(y_0)=0$ and $D^2\phi(y_0)$ is invertible, the Morse lemma applies. It gives local coordinates $z=(z_1,\dots,z_N)$ centered at $y_0$, written as a diffeomorphism $\kappa:W\to V$ with $\kappa(0)=y_0$, in which the phase has the exact form
\begin{align*}
\phi(\kappa(z))=\phi(y_0)+\frac{1}{2}\sum_{r=1}^N \epsilon_r z_r^2.
\end{align*}
The signs $\epsilon_r$ record the inertia of the Hessian: positive directions give $\epsilon_r=1$, and negative directions give $\epsilon_r=-1$. This is why the final answer contains the signature $\operatorname{sgn}(H)$.
Now apply the change-of-variables theorem to the near integral. The substitution is $y=\kappa(z)$, the domain $V$ becomes $W$, and the [Lebesgue measure](/page/Lebesgue%20Measure) transforms as
\begin{align*}
d\mathcal{L}^N(y)=|\det J\kappa_z|\,d\mathcal{L}^N(z).
\end{align*}
After this substitution, the smooth Jacobian factor is part of the amplitude. Thus we define
\begin{align*}
b(z):=\chi(\kappa(z))a(\kappa(z))|\det J\kappa_z|.
\end{align*}
Then $b\in C_c^\infty(W;\mathbb{C})$, because $\chi a$ is smooth, $\operatorname{supp}(\chi a)$ is compactly contained in $\kappa(W)$, and $\kappa$ is smooth. The near integral becomes
\begin{align*}
I_{\mathrm{near}}(\lambda)=e^{i\lambda\phi(y_0)}\int_W e^{i\lambda Q(z)/2}b(z)\,d\mathcal{L}^N(z),
\end{align*}
with
\begin{align*}
Q(z):=\sum_{r=1}^N \epsilon_r z_r^2.
\end{align*}
Finally, the linear part of the Morse coordinate change can be chosen to diagonalize and normalize the Hessian. With this normalization,
\begin{align*}
|\det J\kappa_0|=|\det H|^{-1/2}.
\end{align*}
Since $\chi(y_0)=1$ and $\kappa(0)=y_0$, the transformed amplitude satisfies
\begin{align*}
b(0)=a(y_0)|\det H|^{-1/2}.
\end{align*}
This identity is exactly the source of the factor $|\det H|^{-1/2}$ in the leading term.
[/guided]
[/step]
[step:Extend the quadratic integral without changing its asymptotics]
Since $\operatorname{supp}b$ is a compact subset of $W$, extending $b$ by $0$ outside $W$ defines $b_{\mathrm{ext}}\in C_c^\infty(\mathbb{R}^N;\mathbb{C})$ by $b_{\mathrm{ext}}(z)=b(z)$ for $z\in W$ and $b_{\mathrm{ext}}(z)=0$ for $z\notin W$. Then
\begin{align*}
\int_W e^{i\lambda Q(z)/2}b(z)\,d\mathcal{L}^N(z)=\int_{\mathbb{R}^N}e^{i\lambda Q(z)/2}b_{\mathrm{ext}}(z)\,d\mathcal{L}^N(z).
\end{align*}
We write $b$ for this compactly supported extension. This change is only notational: the value of the integral is unchanged because the extended function is zero outside $W$.
[/step]
[step:Compute the asymptotic coefficients for the quadratic phase]
Define the constant
\begin{align*}
\Gamma_\epsilon := \prod_{r=1}^N \lim_{\delta\downarrow0}\int_{\mathbb{R}} e^{i\epsilon_r t^2/2}e^{-\delta t^2/2}\,d\mathcal{L}^1(t).
\end{align*}
Thus each one-dimensional oscillatory Fresnel integral is defined by Abel regularisation, not by integration against a separate measure. The one-dimensional Fresnel formula, with each integral interpreted by the preceding Abel limit, gives
\begin{align*}
\Gamma_\epsilon=(2\pi)^{N/2}\exp\left(\frac{i\pi}{4}\sum_{r=1}^N\epsilon_r\right).
\end{align*}
Since $\sum_{r=1}^N\epsilon_r=\operatorname{sgn}(H)$, this is
\begin{align*}
\Gamma_\epsilon=(2\pi)^{N/2}e^{i\pi\operatorname{sgn}(H)/4}.
\end{align*}
Define the constant-coefficient differential operator $P$ on $C^\infty(\mathbb{R}^N)$ by
\begin{align*}
P:=\sum_{r=1}^N \epsilon_r\partial_{z_r}^2.
\end{align*}
For every integer $M\geq 0$, the quadratic [stationary phase lemma](/theorems/636) for compactly supported smooth amplitudes applies to the real nondegenerate quadratic form $Q$ and the amplitude $b\in C_c^\infty(\mathbb{R}^N;\mathbb{C})$. It gives a constant $C_{M,b,Q}>0$ and a threshold $\lambda_{M,b,Q}\geq1$ such that, for all $\lambda\geq\lambda_{M,b,Q}$,
\begin{align*}
\left|\int_{\mathbb{R}^N}e^{i\lambda Q(z)/2}b(z)\,d\mathcal{L}^N(z)-\lambda^{-N/2}\Gamma_\epsilon\sum_{j=0}^{M-1}\lambda^{-j}\frac{i^j}{2^j j!}(P^jb)(0)\right|\leq C_{M,b,Q}\lambda^{-N/2-M}.
\end{align*}
This is the precise external quadratic stationary phase result used in the proof. Its hypotheses are satisfied here: $Q$ is a real nondegenerate quadratic form on $\mathbb{R}^N$, because each $\epsilon_r$ is either $1$ or $-1$, and $b\in C_c^\infty(\mathbb{R}^N;\mathbb{C})$ by the preceding extension step. In one proof of the lemma, the away-from-zero estimate uses [integration by parts](/theorems/210) with the first-order differential operator
\begin{align*}
\mathcal{A}:=\frac{1}{i|z|^2}\sum_{r=1}^N \epsilon_r z_r\partial_{z_r},
\end{align*}
which satisfies $\mathcal{A}(e^{iQ(z)/2})=e^{iQ(z)/2}$ on $\mathbb{R}^N\setminus\{0\}$. If $\mathcal{A}^*$ denotes the formal adjoint of $\mathcal{A}$ with respect to Lebesgue measure $\mathcal{L}^N$, then compact support gives
\begin{align*}
\int e^{iQ(z)/2}g(z)\,d\mathcal{L}^N(z)=\int e^{iQ(z)/2}(\mathcal{A}^*g)(z)\,d\mathcal{L}^N(z)
\end{align*}
for every $g\in C_c^\infty(\mathbb{R}^N\setminus\{0\};\mathbb{C})$. Iterating this identity is the integration-by-parts mechanism behind the stated remainder bound. The constants in the lemma depend on $M$, $Q$, a compact set containing $\operatorname{supp}b$, and finitely many derivatives of $b$ up to an order determined by $M$ and $N$.
[guided]
After the Morse change of variables, the phase is exactly quadratic, so the problem is reduced to understanding
\begin{align*}
\int_{\mathbb{R}^N}e^{i\lambda Q(z)/2}b(z)\,d\mathcal{L}^N(z).
\end{align*}
The natural scaling is $z=\lambda^{-1/2}w$. Under this substitution, the Lebesgue measure transforms as
\begin{align*}
d\mathcal{L}^N(z)=\lambda^{-N/2}\,d\mathcal{L}^N(w),
\end{align*}
and the quadratic form satisfies $Q(\lambda^{-1/2}w)=\lambda^{-1}Q(w)$. Therefore the integral becomes
\begin{align*}
\lambda^{-N/2}\int_{\mathbb{R}^N}e^{iQ(w)/2}b(\lambda^{-1/2}w)\,d\mathcal{L}^N(w).
\end{align*}
This scaling explains the universal prefactor $\lambda^{-N/2}$.
We now expand the amplitude at $0$. For a fixed integer $M\geq0$, [Taylor's theorem](/theorems/827) with remainder expresses $b(\lambda^{-1/2}w)$ as a finite sum of homogeneous Taylor polynomials of degrees less than $2M$, plus a remainder whose derivatives are controlled by finitely many derivatives of $b$. Odd monomials integrate to zero against the even oscillatory kernel $e^{iQ(w)/2}$, so only terms of even degree contribute to the expansion through order $\lambda^{-M+1}$.
The Gaussian moments are encoded by the differential operator
\begin{align*}
P:=\sum_{r=1}^N\epsilon_r\partial_{z_r}^2.
\end{align*}
The Fresnel moment identity says that, for each integer $j\geq0$, the total contribution of the degree $2j$ Taylor terms is
\begin{align*}
\Gamma_\epsilon\frac{i^j}{2^j j!}(P^jb)(0),
\end{align*}
where
\begin{align*}
\Gamma_\epsilon=\prod_{r=1}^N\lim_{\delta\downarrow0}\int_{\mathbb{R}}e^{i\epsilon_r t^2/2}e^{-\delta t^2/2}\,d\mathcal{L}^1(t).
\end{align*}
The Abel-regularised one-dimensional Fresnel formulas give
\begin{align*}
\lim_{\delta\downarrow0}\int_{\mathbb{R}}e^{it^2/2}e^{-\delta t^2/2}\,d\mathcal{L}^1(t)=(2\pi)^{1/2}e^{i\pi/4}
\end{align*}
and
\begin{align*}
\lim_{\delta\downarrow0}\int_{\mathbb{R}}e^{-it^2/2}e^{-\delta t^2/2}\,d\mathcal{L}^1(t)=(2\pi)^{1/2}e^{-i\pi/4}.
\end{align*}
Multiplying over all coordinates yields
\begin{align*}
\Gamma_\epsilon=(2\pi)^{N/2}\exp\left(\frac{i\pi}{4}\sum_{r=1}^N\epsilon_r\right)=(2\pi)^{N/2}e^{i\pi\operatorname{sgn}(H)/4}.
\end{align*}
We now invoke the quadratic [stationary phase lemma](/theorems/645) in its finite-remainder form. The lemma applies because $Q$ is a real nondegenerate quadratic form and $b\in C_c^\infty(\mathbb{R}^N;\mathbb{C})$. It states that for each integer $M\geq0$ there are constants $C_{M,b,Q}>0$ and $\lambda_{M,b,Q}\geq1$ such that, for all $\lambda\geq\lambda_{M,b,Q}$,
\begin{align*}
\left|\int_{\mathbb{R}^N}e^{i\lambda Q(z)/2}b(z)\,d\mathcal{L}^N(z)-\lambda^{-N/2}\Gamma_\epsilon\sum_{j=0}^{M-1}\lambda^{-j}\frac{i^j}{2^j j!}(P^jb)(0)\right|\leq C_{M,b,Q}\lambda^{-N/2-M}.
\end{align*}
The proof of that lemma is where the Taylor remainder and the expanding support after $z=\lambda^{-1/2}w$ are controlled. Its integration-by-parts step uses the operator
\begin{align*}
\mathcal{A}_w:=\frac{1}{i|w|^2}\sum_{r=1}^N\epsilon_r w_r\partial_{w_r},
\end{align*}
which satisfies $\mathcal{A}_w(e^{iQ(w)/2})=e^{iQ(w)/2}$ on $\mathbb{R}^N\setminus\{0\}$. If $\mathcal{A}_w^*$ is its formal adjoint with respect to $\mathcal{L}^N$, then compact support away from $0$ permits repeated use of
\begin{align*}
\int e^{iQ(w)/2}g(w)\,d\mathcal{L}^N(w)=\int e^{iQ(w)/2}(\mathcal{A}_w^*g)(w)\,d\mathcal{L}^N(w).
\end{align*}
Thus the guided argument does not hide an extra assumption: the cited lemma supplies exactly the coefficient formula and the $O(\lambda^{-N/2-M})$ remainder needed here.
[/guided]
[/step]
[step:Transfer the quadratic coefficients back to the original amplitude]
For each integer $j\geq0$, define $L_j a$ by
\begin{align*}
L_j a:=|\det H|^{1/2}\frac{i^j}{2^j j!}(P^jb)(0),
\end{align*}
where $b(z)=\chi(\kappa(z))a(\kappa(z))|\det J\kappa_z|$ in the Morse coordinates constructed above. Since $\chi=1$ on a neighbourhood of $y_0$, every derivative of $\chi\circ\kappa$ at $0$ of positive order is zero and $(\chi\circ\kappa)(0)=1$. Therefore $(P^jb)(0)$ depends only on the derivatives at $0$ up to order $2j$ of the map
\begin{align*}
z\mapsto a(\kappa(z))|\det J\kappa_z|.
\end{align*}
By the chain rule, this uses derivatives of $a$ at $y_0$ of order at most $2j$ and derivatives of the chosen Morse coordinate map $\kappa$ at $0$ up to finite order. Thus, for the fixed phase $\phi$ and fixed Morse coordinates, $L_j$ is a linear differential functional of order at most $2j$, evaluated at $y_0$. The theorem only requires existence of such functionals; the displayed construction provides them. If two cutoffs or Morse coordinate maps are used, applying the already proved expansion with the two choices and subtracting gives an asymptotic expansion of the zero function. Uniqueness of asymptotic coefficients then forces equality of the numerical coefficients for each fixed amplitude $a$.
For $j=0$,
\begin{align*}
L_0a=|\det H|^{1/2}b(0)=|\det H|^{1/2}a(y_0)|\det H|^{-1/2}=a(y_0).
\end{align*}
[/step]
[step:Assemble the expansion and absorb the rapidly decaying part]
Combining the quadratic expansion with the definition of $L_j$ gives, for every integer $M\geq0$,
\begin{align*}
I_{\mathrm{near}}(\lambda)=e^{i\lambda\phi(y_0)}e^{i\pi\operatorname{sgn}(H)/4}\left(\frac{2\pi}{\lambda}\right)^{N/2}|\det H|^{-1/2}\sum_{j=0}^{M-1}\lambda^{-j}L_j a+O(\lambda^{-N/2-M}).
\end{align*}
From the first step, for every integer $m\geq0$,
\begin{align*}
I_{\mathrm{far}}(\lambda)=O(\lambda^{-m}).
\end{align*}
Choosing $m>N/2+M$ shows that the far contribution is absorbed into the same remainder. Therefore
\begin{align*}
I(\lambda)=e^{i\lambda\phi(y_0)}e^{i\pi\operatorname{sgn}(H)/4}\left(\frac{2\pi}{\lambda}\right)^{N/2}|\det H|^{-1/2}\sum_{j=0}^{M-1}\lambda^{-j}L_j a+O(\lambda^{-N/2-M}).
\end{align*}
For the explicit remainder formulation, fix $M\geq0$. The quadratic lemma gives constants $C_{M,b,Q}>0$ and $\lambda_{M,b,Q}\geq1$ for the near contribution. Choose an integer $m_M>N/2+M$. The nonstationary phase estimate with this integer gives a constant $C_{m_M}>0$, depending on $M$, $\phi$, and the compactly supported amplitude $(1-\chi)a$, such that $|I_{\mathrm{far}}(\lambda)|\leq C_{m_M}\lambda^{-m_M}$ for all $\lambda\geq1$. Since $m_M>N/2+M$ and $\lambda\geq1$, this implies $|I_{\mathrm{far}}(\lambda)|\leq C_{m_M}\lambda^{-N/2-M}$. Therefore, with $\lambda_M:=\lambda_{M,b,Q}$ and with $C_M:=C_{M,b,Q}+C_{m_M}$ enlarged if necessary to absorb the fixed factors in the near expansion, the stated inequality holds for every $\lambda\geq\lambda_M$. Since this holds for every integer $M\geq0$, it is precisely the asserted full asymptotic expansion.
[/step]