An invariant method of stationary phase 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).
 A: Check Proposition 1.2.4 from the book Fourier Integral Operators by the late great  Duistermaat. This result applies  in the case when the phase $\phi$ is Morse.  If the phase is not Morse, but the critical points    are finitely determined (finite Milnor number) then things are a bit more complicated. The vol 1 book  of Arnold-Gusein-Zade Varchenko Singularities of differentiable maps is a good source.      You can also  have a look  at the senior thesis of a Zach Lamberty, a  former student of mine.  There he deals with the $2$-dimensional case ($\dim M =2$) and he essentially works out the toric resolution trick of Arnold and comp.  for a special and quite degenerate two variable phase.
A: $\phi=\Re \phi+i\Im \phi$. You have assumed $\Re \phi\ge 0$ and you deal with a complex phase function. Note that the standard notation is not yours, since what is usually called the real stationary phase method coincides here with the case $\phi$ purely imaginary. 
Never mind, let's follow your notations and note that $t\rightarrow+\infty$ (the $+$ is missing in your formulation and is quite important since $\Re \phi\ge 0$).
(1) Let us assume that $\Im \phi$ does not have a stationary point on the support of the amplitude $u$($d\Im \phi\not=0$ on $supp u\cap${$\Re \phi= 0$}: then $I(\phi, u,t)=O(t^{-N})$ for any $N>0$.
(2) Let us assume that $\Im \phi$ is such that
$$
d\Im \phi=0\text{ at}\quad  supp u\cap{\text{{$\Re \phi= 0$}}\text{}}\Longrightarrow
\text{Hessian}(\Im \phi) \text{ non-singular}
$$
then 
$I(\phi, u,t)\sim ct^{-n/2}$, where $n=dim M$.
The constant $c$ can be computed explicitly in terms of the indices of the Hessian at the stationary points, the value of the amplitude there and appears as a finite sum corresponding to the finite number of stationary points of $\Im \phi$ on the compact $supp u\cap${$\Re \phi= 0$}.
A  complete expansion is available in Hormander ALPDO first volume, Chapter 7, section devoted to the complex stationary phase method. To sum-up the simple case exposed here: the integral is largest when $\Re \phi$ vanishes at a critical point of $\Im \phi$, and if that critical point is non-degenerate, you find a behavior in $t^{-n/2}$.
No coordinate choice is involved here.
I should say that the real stationary phase method is easier to understand: it corresponds here to your case with $\phi$ purely imaginary (!). You may for instance assume that $i\phi$ is a real-valued Morse function and the Morse lemma is providing a normal form (No such thing exists for a complex valued function). You find a finite number of stationary points on the support of $u$, and you can take advantage of the normal form on a neighborhood of each critical point. Anyhow the contribution elsewhere is $O(t^{-\infty}).$
Morse lemma reduces the problem to an integral with an exactly quadratic phase for which you have a full expansion since you know explicitly the Fourier transform of a Gaussian function.
