[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]