Steepest descent integration in several dimensions The method of steepest descent provides an asymptotic approximation for integrals of the form:
$$I = \int_C \exp(M f(z))\mathrm dz$$
for large positive $M$, where $f(z)$ is analytic in the region of interest, $C$ a contour and $f(z)$ goes to zero at the endpoints of the contour. The asymptotic approximation is:
$$I \sim \sqrt{\frac{2\pi}{f''(z_0)}}\exp(M f(z_0)+i\theta)$$
where $z_0$ is a saddle-point of $f$ ($f'(z_0)=0$) such that the original contour $C$ can be deformed (fixing the endpoints) to pass through $z_0$ in a steepest descent direction (which defines the angle $\theta$) and assuming that it doesn't go through any other saddle-points for simplicity.
Is there a generalization for higher-dimensional integrals? For example,
$$J = \int_{C_1}\mathrm dz_1\int_{C_2}\mathrm dz_2 \exp(M g(z_1,z_2))$$
I haven't studied much of multi-dimensional complex analysis, so I'm not sure what conditions should be imposed on $g(z_1,z_2)$. But as suggested in (https://en.wikipedia.org/wiki/Several_complex_variables), say that $g(z_1,z_2$ can be represented as a convergent power series in the region of interest.
 A: You can iterate the the integrals assuming there is a saddle point, using the implicit function theorem and a bit of work, as explained here.
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.661.8737&rep=rep1&type=pdf
The key is expressing the  the variables in terms of each other (using implicit function theorem), and having a non vanishing Hessian determinant. 
A: You can look for Multivariable Morse Lemma to get an extension of Steepest Descent into multiple complex variables $z=[z_1, z_2, ...,z_n]$. Higher dimensional asymptotics as $M\rightarrow\infty$ for this multiple integral with $f:\mathbb{C}^n\rightarrow\mathbb{C}$ and $C=C_1\times C_2\times ...\times C_n$ a multiple complex contour domain, $$J(M)=\int_C e^{M f(z)}dz$$ is obtained taking $z_0=[z_{10}, z_{20}, ...,z_{n0}]$, the point where $f'(z_0)=\nabla f(z)|_{z_0}=0$  then, under regularity conditions (single non-degenerate saddle point),  $$J(M)=\frac{(2\pi/M)^{\frac{n}{2}}e^{Mf(z_0)}}{\sqrt{\det{[-f''(z_0)]}}}[1+O(1/M)]$$ here $f''(z)$ is the Hessian matrix having eigenvalues $\lambda=[\lambda_1,\lambda_2,...,\lambda_n]$ with $|\arg(-\lambda_k)|<\frac{\pi}{2}$ $$\sqrt{\det{[-f''(z_0)]}}=e^{\frac{i}{2}\phi}\prod_{k=1}^n|\lambda_k|^{1/2}\ne0$$ where $\phi=\sum_{k=1}^n \arg(-\lambda_k)$. Note that if $z_0\in\mathbb{R}^n$ and $\Im[f(z)]=0\  \ \forall z\in\mathbb{R}^n\ $  then  $\ \phi=0$.
Also if $z_0\in\mathbb{R}^n,\  \Re[f(z)]=0\ \ \forall z\in\mathbb{R}^n\ $ and $|\arg\sqrt{-\lambda_k}|\le\frac{\pi}{4}$ then  $\ \phi=m\cdot\frac{\pi}{2}\ $ where $m$ is the number of negative eigenvalues minus the number of positive ones (stationary phase method).
The case of multiple saddle points $z_0^{(\ell)}$ s.t. $f'(z_0^{(\ell)})=0$ and $\det{[-f''(z_0^{(\ell)})]}\ne0,\ \ \ell=1,2,...L$ is worked the same way. This gives $$J(M)=\sum_{\ell=1}^L\frac{(2\pi/M)^{\frac{n}{2}}e^{Mf(z_0^{(\ell)})}}{\sqrt{\det{[-f''(z_0^{(\ell)})]}}}[1+O(1/M)]$$
Chapters VIII-IX in
Wong, R., Asymptotic approximations of integrals, Classics in Applied Mathematics 34. Philadelphia, PA: SIAM (ISBN 0-89871-497-4/pbk). xvii, 543 p. (2001). ZBL1078.41001.
provide asymptotic methods for multidimensional integrals. Cases of degenerate saddles where Hessian vanishes $\det[-f''(z_0^{(\ell)})]=0$ for some $\ell$ is more complex. They are found in Ch VIII sect. 5 pg. 435 and Ch IX sect. 4 pg 491.
I hope this helps.
A: You can find a referenced C^N saddle point method here:
Pinna, Francesco; Viola, Carlo, The saddle-point method in (\mathbb{C}^N) and the generalized Airy functions, Bull. Soc. Math. Fr. 147, No. 2, 221-257 (2019). ZBL1472.41018.
See also the link at the publisher. There are also some explained examples.
