The method of stationary phase is very well-known and employed in many areas of physics and mathematics, and, of course, included in various versions as theorem in textbooks, especially on pseudors and microlocal analysis.

However, it always is somewhat dependent on local coordinates and the Fourier transform, despite being a quite invariant problem. To be precise, the question would be the following.

Let $M$ be a manifold and $\phi: M \longrightarrow \mathbb{C}$ be a smooth function with values in the closed right half plane. Let $u$ be an volume density on $M$ with compact support in $M$. Determine an asymptotic expansion as $t \rightarrow \infty$ of the integral $$ I(\phi, u, t) = \int_M e^{-\phi t} u$$ under some nondegeneracy conditions on $\phi$.

(For example, one could require $\phi$ to be Morse or, more general, require that the set where it vanishes is a submanifold $C$ of $M$ and that at a point $p \in C$, the Hessian of $\phi$ is non-degenerate on the space $T_pM/T_pC$.)

It is well-known that in these cases $I(\phi, u, t)$ has an asymptotic expansion of the form
$$ I(\phi, u, t) = (t/\pi)^{-(n-k)/2}\sum_{j=0}^\infty t^{-j} \int_C s_j,$$
where $k$ is the dimension of $C$ and the $s_j$ are certain volume densities on $C$. **In fact, they have to be certain universal terms, depending only on the $2j$-th jets of $\phi$ and $u$ at $C$.** This is not stated in most textbooks.

I wonder if it is possible to find these terms $s_j$ using Invariance theory alone. I would like if someone ever thought about this and knows a reference to this more invariant, geometric approach.

/Edit: To clarify my question: I was wondering if it is possible to determine the constants by **invariance theory**, i.e. some argument like "there is only one polynomial on the $2j$-jets of $u$ and $\phi$ that is invariant under coordinate transformation" or so. For the first term, this goes like this, supposed that $\phi$ is purely real:

Define the $n-k$-density $\mathrm{H}\phi$ on $C$ by setting $$\mathrm{H}\phi[X_1, \dots, X_{n-k}] := \sqrt{\left|\det \bigl( D^2\phi[X_i, X_j] \bigr)_{ij}\right|},$$ where $D^2\phi$ is the (on $C$ well-defined) Hessian of $\phi$. Now $u/\mathrm{H}\phi$ is a $k$-density on $C$ -- this is $s_0$.

Now there should be similar characterizations of the higher $s_j$ (which obviously can get arbitrarily complicated).