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.