I would like to evaluate the integral asymptotically over the unit sphere surface $$ Z =\int e^{a \cos^2 \theta + b \sin^2\theta\cos2\phi + c\cos\theta} d\Omega = \int\limits_{0}^{\pi}\int\limits_{0}^{2\pi} e^{a \cos^2 \theta + b \sin^2\theta\cos2\phi + c\cos\theta} \sin\theta d\phi d\theta $$ for $a\rightarrow \pm\infty$ and $b\rightarrow\infty$, if $c$ is set implictly by the constraint $$ L = \frac{\partial \ln Z}{\partial c}$$ for a fixed $0<L<1$. With the variables change $s=\cos\theta$, \begin{align} Z &= 2\pi\int_{-1}^1 e^{a s^2+c s}\, I_0[b (1-s^2)] ds\,, \\ L &= \frac{\int_{-1}^1 s\, e^{a s^2+c s} I_0[b (1-s^2)]\, ds }{\int_{-1}^1 e^{a s^2+c s} I_0[b (1-s^2)]\, ds }\,. \end{align}
where $I_0(x)$ is the Bessel function of the first kind. (Alternatively, the $\theta$ integral may be evaluated for any $\phi$ using the error function.)
Is there a general method to look for the $a\rightarrow \pm\infty$ asymptotics of this integral?
Note that this is the integral of a Gaussian on the sphere in the sense that $$Z_1=\int_{\|\mathbf{x}\|=1} e^{-{\mathbf{x}\cdot \mathbf{M x} + \mathbf{v}\cdot\mathbf{x}}} dS$$ where the integration is over the unit sphere's surface, $\mathbf{M}$ is a $3\times 3$ symmetric traceless matrix, and $\mathbf{v}$ is a vector parallel to one of the eigenvectors. The integral in spherical coordinates $\mathbf{x}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$ and $dS = \sin\theta d\theta d\phi$ aligned with the eigenvectors, reduces to the integral given in the question as $Z_1=e^{-a/3} Z$, where $a=\frac32\lambda_3$, $b=\frac12(\lambda_1-\lambda_2)$, $\mathbf{v}=c\, \mathbf{u}_3$, and $\lambda_{1,2,3}$ and $\mathbf{u}_{1,2,3}$ are the three eigenvalues and normalized eigenvectors of $\mathbf{M}$.
The following is just the sketch of an idea. With some luck, somebody will be able to complete it (or show that it is wrong, of course!)
Differentiate with respect to $c$ to obtain \begin{align} \frac{\partial Z}{\partial c} = & \int_0^\pi\int_0^{2\pi}\left(e^{a\cos^2\theta}\cos\theta\sin\theta\right)e^{b\sin^2\theta\cos2\phi+c\cos\theta}d\phi d\theta\\ = & \int_0^\pi\int_0^{2\pi}\left(-\frac{1}{2a}\frac{\partial e^{a\cos^2\theta}}{\partial\theta}\right)e^{b\sin^2\theta\cos2\phi+c\cos\theta}d\phi d\theta. \end{align} Integrating by parts, we obtain \begin{align} \frac{\partial Z}{\partial c} = &\left[-\frac{1}{2a}\int_0^{2\pi}e^{a\cos^2\theta+b\sin^2\theta\cos2\phi+c\cos\theta}d\phi\right]_0^\pi\\ & +\frac{1}{2a}\int_0^\pi\int_0^{2\pi}e^{a\cos^2\theta+b\sin^2\theta\cos2\phi+c\cos\theta}\left(2b\sin\theta\cos\theta\cos2\phi - c\sin\theta\right)d\phi d\theta. \end{align} Now hopefully the term on the second line will go to zero as $a\to\infty$, and the first term can be evaluated using Bessel functions and the error function.
Edit: As remarked by @tired, the second term will not go to zero, as the exponential is increasing too fast as $a\to\infty$. I will leave this here anyway, so if someone wants to try something passing from partial integration as I did above they won't have to do the computation over again.