asymptotic with a very degenerate stationary phase Suppose $f(x_1,x_2)\in C^\infty_c(\mathbb R^2)$. I wonder how one may derive the asymptotic expansion of the following integral when the real paramter $\lambda\rightarrow \infty$:
\begin{equation}
\int_{\mathbb R^2}e^{i\lambda x_1x_2(x_1+x_2)}f(x_1,x_2) dx_1 dx_2.
\end{equation}
Note that the origin is the only critical point of the phase function $\phi(x_1,x_2)=x_1x_2(x_1+x_2)$, but it is a very degenerate one as the Hessian matrix $\big(\frac{\partial^2 \phi}{\partial x_i\partial x_j}\big)$ is zero at the origin. Actually, the Milnor number that describes the singularity of $\phi$ at the origin is $4$.
[EDIT: I have found an answer for this question myself. When the phase function $\phi$ is as simple as stated above, the integral's asymptotic behavior can be written down in a clean manner. The real complexity shows up when there is some perturbation in the phase, that is, when the phase is of the form $\Phi(x_1,x_2)=x_1x_2(x_1+x_2)+ax_1x_2+bx_1+cx_2$, with $a,b,c$ denoting the control parameters. In the latter case, an asymptotic $\lambda$-expansion that is uniform in $a,b,c$ is very hard to get and seems not existent in the literature yet.]
 A: (I am not really an expert, so it should be taken with a pinch of salt.)  
One may take $f=1$ and simply compute the integral using polar coordinates. (There are minor problems with convergence, but they are irrelevant.) The result is $c|\lambda|^{-2/3}$, where $c\neq 0$ is an absolute constant which I was too lazy to compute. So, in the general case the asymptotic must be $c|\lambda|^{-2/3}f(0,0)$. 
(To actually prove this may be a piece of work, but it certainly can be done if necessary.) 
[EDIT] Now I think I can help. First, I would change coordinates to make the phase function more symmetric, 
$\phi(x,y)=y^3-3x^2 y$. Then, I would make use of complex variables.  Let  $\Sigma\subset\mathbb{C}^2$ 
be a surface defined by $x=(\cos t-iu\sin 2t)r,\,y=(\sin t -iu\cos 2t)r$, where $u>0$ is not very big. ($u=1/2$ will do.) Here $0\le t<2\pi$ and $r>0$. 
The idea is to replace the original integral with 
$$\int_{\Sigma}e^{i\lambda\phi(x,y)}f(x,y)dx\wedge dy,$$
which must have the same asymptotic. This integral is convergent 
even if $f$ is an entire function, because on this surface $\Im \phi(x,y)>0$. The asymptotic series  simply comes from the Taylor series. (To make the argument rigorous one may use a partition of unity and other standard techniques.)
