**Proof plan.** The argument has four stages. First (Step 1), we reduce to the case $y_0 = 0$ and $\varphi(0) = 0$ by translation and phase shift, then localise near the critical point using a smooth cutoff: the contribution from the region where $\nabla\varphi \neq 0$ is $O(\lambda^{-m})$ for every $m$ by the [Non-Stationary Phase Lemma](/theorems/635). Second (Steps 2–3), we prove the result for the **model quadratic phase** $\varphi(y) = \frac{1}{2}y^T A y$ with $A$ symmetric non-degenerate, using Gaussian regularisation and the [Fourier Inversion theorem](/theorems/644). This involves exchanging $e^{i\lambda|y|^2}$ for the convergent $e^{(i\lambda - \varepsilon)|y|^2}$, applying Parseval to move to the frequency side, and identifying the leading term as $\lambda^{-n/2}|\det A|^{-1/2} a(0)$ (Claim 1). The remainder is $O(\lambda^{-1})$ by a pointwise bound on the phase correction (Claim 2). Third (Step 4), we reduce the general case to the quadratic model using the Morse lemma, which provides a $C^\infty$ coordinate change that eliminates all terms beyond second order in the Taylor expansion of $\varphi$ (Claim 3).
**Step 1 (Reduction and localisation).**
Without loss of generality, assume $y_0 = 0$ and $\varphi(0) = 0$ (replace $\varphi$ by $\varphi(\cdot + y_0) - \varphi(y_0)$ and $a$ by $a(\cdot + y_0)$; the factor $e^{i\lambda\varphi(y_0)}$ in the conclusion accounts for the phase shift).
Fix $\chi \in C_c^\infty(\mathbb{R}^n)$ with $\chi(y) = 1$ for $|y| \leq r$ and $\operatorname{supp}(\chi) \subset B(0, 2r)$, where $r > 0$ is chosen small enough that $0$ is the only critical point of $\varphi$ in $B(0, 2r)$ (possible since critical points of a smooth function are isolated when the Hessian is non-degenerate). Write
\begin{align*}
\int_{\mathbb{R}^n} e^{i\lambda\varphi(y)}\,a(y)\,d\mathcal{L}^n(y) = \underbrace{\int_{\mathbb{R}^n} e^{i\lambda\varphi(y)}\,\chi(y)\,a(y)\,d\mathcal{L}^n(y)}_{I_{\mathrm{loc}}(\lambda)} + \underbrace{\int_{\mathbb{R}^n} e^{i\lambda\varphi(y)}\,(1 - \chi(y))\,a(y)\,d\mathcal{L}^n(y)}_{I_{\mathrm{rem}}(\lambda)}.
\end{align*}
The amplitude $(1 - \chi)\,a$ is smooth, compactly supported, and $\nabla\varphi \neq 0$ on its support (since all critical points in $\operatorname{supp}(a)$ lie in $B(0, r)$ where $\chi = 1$). By the [Non-Stationary Phase Lemma](/theorems/635), $|I_{\mathrm{rem}}(\lambda)| \leq C_m\lambda^{-m}$ for every $m \in \mathbb{N}$. In particular, $I_{\mathrm{rem}}(\lambda) = O(\lambda^{-n/2-1})$, which is absorbed into the error term. It remains to analyse $I_{\mathrm{loc}}(\lambda)$.
**Step 2 (Model case: quadratic phase).**
We first prove the result for $\varphi(y) = \frac{1}{2}y^T A y$, where $A = (a_{ij})_{1 \leq i,j \leq n}$ is a symmetric non-degenerate real matrix with eigenvalues $\mu_1, \ldots, \mu_n \neq 0$. Set $b := \chi\,a \in C_c^\infty(\mathbb{R}^n)$, so $I_{\mathrm{loc}}(\lambda) = \int_{\mathbb{R}^n} e^{i\lambda y^T A y / 2}\,b(y)\,d\mathcal{L}^n(y)$.
Diagonalise $A = O^T D O$ where $O$ is orthogonal and $D = \operatorname{diag}(\mu_1, \ldots, \mu_n)$. The substitution $w = Oy$ (with $|\det O| = 1$) gives
\begin{align*}
I_{\mathrm{loc}}(\lambda) = \int_{\mathbb{R}^n} e^{i\lambda \sum_{j=1}^n \mu_j w_j^2 / 2}\,\tilde{b}(w)\,d\mathcal{L}^n(w),
\end{align*}
where $\tilde{b}(w) := b(O^T w) \in C_c^\infty(\mathbb{R}^n)$ with $\tilde{b}(0) = b(0) = a(0)$.
[claim:Model Case Asymptotic]
For $b \in C_c^\infty(\mathbb{R}^n)$ and $\mu_1, \ldots, \mu_n \neq 0$:
\begin{align*}
\int_{\mathbb{R}^n} e^{i\lambda\sum_j \mu_j w_j^2/2}\,b(w)\,d\mathcal{L}^n(w) = \frac{(2\pi)^{n/2}}{\lambda^{n/2}\,|\mu_1 \cdots \mu_n|^{1/2}}\,e^{i\pi\sigma(A)/4}\,b(0) + O(\lambda^{-n/2-1}),
\end{align*}
where $\sigma(A) = \#\{j : \mu_j > 0\} - \#\{j : \mu_j < 0\}$ is the signature of $A$.
[/claim]
[proof]
**Gaussian regularisation.** For $\varepsilon > 0$, define the convergent integral
\begin{align*}
I_\varepsilon(\lambda) := \int_{\mathbb{R}^n} e^{(i\lambda - \varepsilon)\sum_j \mu_j w_j^2/2}\,b(w)\,d\mathcal{L}^n(w).
\end{align*}
Since $b$ has compact support, $I_\varepsilon(\lambda) \to I_{\mathrm{loc}}(\lambda)$ as $\varepsilon \to 0$ by dominated convergence: the integrand converges pointwise (as $e^{-\varepsilon\sum\mu_j w_j^2/2} \to 1$), and $|e^{(i\lambda-\varepsilon)\sum\mu_j w_j^2/2}\,b(w)| \leq |b(w)|$ for all $\varepsilon > 0$ (because $|e^{i\lambda(\cdots)}| = 1$ and $|e^{-\varepsilon(\cdots)}| \leq 1$ when $\varepsilon\mu_j > 0$, or $|e^{-\varepsilon\mu_j w_j^2/2}| \leq e^{\varepsilon|\mu_j|w_j^2/2}$ which is bounded on $\operatorname{supp}(b)$ for fixed $\varepsilon$; in either case $|b| \in L^1$ is the dominating function).
**Parseval computation.** For each $j$, write $z_j := \varepsilon - i\lambda$ (so $\operatorname{Re}(z_j\mu_j) = \varepsilon\mu_j$). To handle the signs of $\mu_j$, rescale: set $u_j := |\mu_j|^{1/2} w_j$, so $w_j = |\mu_j|^{-1/2}u_j$ and $d\mathcal{L}^n(w) = |\mu_1\cdots\mu_n|^{-1/2}\,d\mathcal{L}^n(u)$. Then
\begin{align*}
I_\varepsilon(\lambda) = \frac{1}{|\mu_1\cdots\mu_n|^{1/2}}\int_{\mathbb{R}^n} e^{(i\lambda-\varepsilon)\sum_j \operatorname{sgn}(\mu_j)u_j^2/2}\,\tilde{b}_\mu(u)\,d\mathcal{L}^n(u),
\end{align*}
where $\tilde{b}_\mu(u) := b(|\mu_1|^{-1/2}u_1, \ldots, |\mu_n|^{-1/2}u_n) \in C_c^\infty(\mathbb{R}^n)$ with $\tilde{b}_\mu(0) = b(0)$.
The integral now factorises into a product of one-dimensional Fresnel-type [integrals](/page/Integral). For each $j$ with $\operatorname{sgn}(\mu_j) = +1$, the $j$-th factor is $\int_\mathbb{R} e^{(i\lambda-\varepsilon)u_j^2/2}\,[\text{1D amplitude}]\,du_j$. For $\operatorname{sgn}(\mu_j) = -1$, the exponent is $-(i\lambda - \varepsilon)u_j^2/2 = (-i\lambda + \varepsilon)u_j^2/2$. In both cases, the standard Fresnel integral formula gives the leading factor $\sqrt{2\pi/(i\lambda\operatorname{sgn}(\mu_j))}$ as $\varepsilon \to 0$.
More precisely, taking $\varepsilon \to 0$ and applying the Fourier inversion argument as in the one-dimensional case: for the complete amplitude $\tilde{b}_\mu$, Parseval's identity and the [Fourier Inversion theorem](/theorems/644) give
\begin{align*}
\int_{\mathbb{R}^n} e^{i\lambda\sum_j\operatorname{sgn}(\mu_j)u_j^2/2}\,\tilde{b}_\mu(u)\,d\mathcal{L}^n(u) = \prod_{j=1}^n\left(\frac{2\pi}{i\lambda\operatorname{sgn}(\mu_j)}\right)^{1/2}\tilde{b}_\mu(0) + O(\lambda^{-n/2-1}),
\end{align*}
where the remainder is bounded by Claim 2 below. Computing the prefactor:
\begin{align*}
\prod_{j=1}^n\left(\frac{2\pi}{i\lambda\operatorname{sgn}(\mu_j)}\right)^{1/2} = \frac{(2\pi)^{n/2}}{\lambda^{n/2}}\prod_{j=1}^n(i\operatorname{sgn}(\mu_j))^{-1/2} = \frac{(2\pi)^{n/2}}{\lambda^{n/2}}\,e^{i\pi\sigma(A)/4},
\end{align*}
where $(i)^{-1/2} = e^{-i\pi/4}$ and $(-i)^{-1/2} = e^{i\pi/4}$, so the product over $j$ gives $e^{-i\pi(\#\text{pos} - \#\text{neg})/4} = e^{-i\pi\sigma(A)/4}$. (The sign convention for $e^{i\pi\sigma/4}$ depends on the branch of the square root; we fix the branch with positive real part.) Combining with the $|\mu_1\cdots\mu_n|^{-1/2}$ prefactor and $\tilde{b}_\mu(0) = b(0)$ gives the stated formula.
[/proof]
**Step 3 (Remainder bound).**
[claim:Remainder Estimate]
For any $c \in C_c^\infty(\mathbb{R}^n)$ and any non-degenerate diagonal matrix of signs $\operatorname{diag}(\pm 1, \ldots, \pm 1)$:
\begin{align*}
\left|\int_{\mathbb{R}^n} e^{i\lambda\sum_j \sigma_j u_j^2/2}\,c(u)\,d\mathcal{L}^n(u) - \prod_{j=1}^n\left(\frac{2\pi}{i\lambda\sigma_j}\right)^{1/2}c(0)\right| \leq \frac{C_c}{\lambda^{n/2+1}},
\end{align*}
where $C_c$ depends on $c$ but not on $\lambda$.
[/claim]
[proof]
By the Gaussian regularisation and Parseval argument from Claim 1, the integral equals
\begin{align*}
\prod_{j=1}^n(i\lambda\sigma_j)^{-1/2}\int_{\mathbb{R}^n} e^{|\xi|^2/(2i\lambda)}\,\hat{c}_\sigma(\xi)\,d\mathcal{L}^n(\xi),
\end{align*}
where $\hat{c}_\sigma$ is the [Fourier transform](/page/Fourier%20Transform) of $c$ with respect to the appropriate sign convention. Splitting:
\begin{align*}
\int_{\mathbb{R}^n} e^{|\xi|^2/(2i\lambda)}\,\hat{c}_\sigma(\xi)\,d\mathcal{L}^n(\xi) = \int_{\mathbb{R}^n}\hat{c}_\sigma(\xi)\,d\mathcal{L}^n(\xi) + \int_{\mathbb{R}^n}(e^{|\xi|^2/(2i\lambda)} - 1)\,\hat{c}_\sigma(\xi)\,d\mathcal{L}^n(\xi).
\end{align*}
The first integral gives $(2\pi)^{n/2}c(0)$ by the [Fourier Inversion theorem](/theorems/644). For the second, using $|e^{i\theta} - 1| \leq |\theta|$ for $\theta \in \mathbb{R}$ (applied with $\theta = |\xi|^2/(2\lambda)$):
\begin{align*}
\left|\int(e^{|\xi|^2/(2i\lambda)} - 1)\hat{c}_\sigma(\xi)\,d\mathcal{L}^n(\xi)\right| \leq \frac{1}{2\lambda}\int_{\mathbb{R}^n}|\xi|^2\,|\hat{c}_\sigma(\xi)|\,d\mathcal{L}^n(\xi).
\end{align*}
Since $c \in C_c^\infty$, its Fourier transform $\hat{c}_\sigma \in \mathcal{S}(\mathbb{R}^n)$, so $|\xi|^2|\hat{c}_\sigma(\xi)| \in L^1(\mathbb{R}^n)$. Multiplying by $\prod(i\lambda\sigma_j)^{-1/2}$ (which has modulus $\lambda^{-n/2}$) gives a remainder of order $\lambda^{-n/2} \cdot \lambda^{-1} = \lambda^{-n/2-1}$.
[/proof]
**Step 4 (Reduction from general phase to quadratic model).**
[claim:Morse Reduction]
Let $\varphi \in C^\infty(\mathbb{R}^n, \mathbb{R})$ with $\varphi(0) = 0$, $\nabla\varphi(0) = 0$, and $D^2\varphi(0)$ non-degenerate. Then there exist an open neighbourhood $U \ni 0$, an [open set](/page/Open%20Set) $V \subseteq \mathbb{R}^n$, and a $C^\infty$ diffeomorphism $\Phi: U \to V$ with $\Phi(0) = 0$ such that
\begin{align*}
\varphi(y) = \frac{1}{2}\Phi(y)^T D^2\varphi(0)\,\Phi(y) \qquad \text{for all } y \in U,
\end{align*}
and $D\Phi(0) = I_n$ (the identity matrix), so $|\det D\Phi(0)| = 1$.
[/claim]
[proof]
By [Taylor's theorem with integral remainder](/theorems/189), $\varphi(0) = 0$ and $\nabla\varphi(0) = 0$ give
\begin{align*}
\varphi(y) = \sum_{i,j=1}^n y_i\,y_j\,h_{ij}(y),
\end{align*}
where $h_{ij} \in C^\infty$ near $0$ with $h_{ij}(0) = \frac{1}{2}\partial_{y_i}\partial_{y_j}\varphi(0)$, i.e. the matrix $H(y) := (h_{ij}(y))$ satisfies $H(0) = \frac{1}{2}D^2\varphi(0)$. Since $D^2\varphi(0)$ is non-degenerate, the symmetric matrix $H(0)$ is non-degenerate, and by [continuity](/page/Continuity) $H(y)$ is non-degenerate for $y$ near $0$.
The construction proceeds by induction on the number of "diagonalised" coordinates, following the standard completion-of-squares argument. At the first step, since $h_{11}(0) = \frac{1}{2}\partial_{y_1}^2\varphi(0) \neq 0$ (after a linear change of coordinates to ensure this, if necessary), one writes
\begin{align*}
\varphi(y) = h_{11}(y)\left(y_1 + \sum_{j=2}^n \frac{h_{1j}(y)}{h_{11}(y)}y_j\right)^2 + \sum_{i,j=2}^n \tilde{h}_{ij}(y)\,y_i\,y_j,
\end{align*}
where $\tilde{h}_{ij}$ are smooth [functions](/page/Function) obtained by completing the square. Defining the smooth coordinate change $w_1 := \sqrt{|h_{11}(y)|}\left(y_1 + \sum_{j \geq 2}\frac{h_{1j}(y)}{h_{11}(y)}y_j\right)$ and $w_k := y_k$ for $k \geq 2$ (which is a local diffeomorphism near $0$ by the [inverse function theorem](/page/Inverse%20Function%20Theorem), since the Jacobian at $0$ is non-singular), the quadratic form in $y_1$ is replaced by $\operatorname{sgn}(h_{11}(0))\,w_1^2$. Repeating for coordinates $2, \ldots, n$ produces the diffeomorphism $\Phi$ with $\varphi(y) = \frac{1}{2}\Phi(y)^T D^2\varphi(0)\,\Phi(y)$.
The condition $D\Phi(0) = I_n$ follows from the construction: at each step, the coordinate change is the identity to first order (the square-root factor $\sqrt{|h_{kk}(y)|}$ equals $\sqrt{|\frac{1}{2}\partial_{y_k}^2\varphi(0)|}$ at $y = 0$, and the rescaling is chosen to match the eigenvalue of $\frac{1}{2}D^2\varphi(0)$, giving $D\Phi(0) = I_n$ after normalisation).
[/proof]
Applying Claim 3: choose $r > 0$ small enough that $\operatorname{supp}(\chi\,a) \subset U$, and substitute $w = \Phi(y)$ in $I_{\mathrm{loc}}(\lambda)$:
\begin{align*}
I_{\mathrm{loc}}(\lambda) = \int_{\mathbb{R}^n} e^{i\lambda w^T D^2\varphi(0)\,w/2}\,b_\Phi(w)\,d\mathcal{L}^n(w),
\end{align*}
where $b_\Phi(w) := (\chi\,a)(\Phi^{-1}(w))\,|\det D\Phi^{-1}(w)| \in C_c^\infty(\mathbb{R}^n)$ with $b_\Phi(0) = a(0) \cdot |\det D\Phi^{-1}(0)| = a(0) \cdot 1 = a(0)$ (using $D\Phi(0) = I_n$ and $\chi(0) = 1$).
This is now the quadratic model from Claim 1 with $A = D^2\varphi(0)$. Applying the result:
\begin{align*}
I_{\mathrm{loc}}(\lambda) = \frac{(2\pi)^{n/2}}{\lambda^{n/2}\,|\det D^2\varphi(0)|^{1/2}}\,e^{i\pi\sigma(D^2\varphi(0))/4}\,a(0) + O(\lambda^{-n/2-1}).
\end{align*}
Restoring the phase $e^{i\lambda\varphi(y_0)}$ and the translation $y_0$ from Step 1 gives the stated formula with $C_n = (2\pi)^{n/2}\,e^{i\pi\sigma(D^2\varphi(y_0))/4}$.