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