I'm looking for a way to take a stationary phase approximation of an integral of the following form:
$$ \int_{-\infty}^\infty d\vec{q} \exp\left(2 \pi i N \left(S(q_{n+1}, \vec{q}, q_1) - \vec{K}^T\vec{q}\right)\right)$$
Here, $\vec{K}$ is just an integer valued constant vector, $\vec{q}$ is the vector $(q_2, q_3, \dotsc, q_n)$ and $S(q_{n+1}, \vec{q}, q_1)$ is just some scalar function that corresponds to the action of the system i'm looking at.
It has been suggested to me that this integral is highly oscillatory and thus the method of stationary phase can be utilised. However, I have no idea how to do stationary phase approximations for a multi-dimensional integral. Taking derivatives w.r.t $\,\vec{q}$ of the exponent suggests that stationary points occur when $$\frac{\partial S}{\partial \vec{q}} = \vec{K}^T.$$ Given I can find these, what form would the stationary phase approximation take? Googling i've only been able to find references for 1-dimensional integrals, and at most a passing reference that the method can be extended to multi-dimensional integrals.