I'm trying to calculate $\mathcal{I} = \int_{S^2} e^{-(Ax^2 + By^2)} dS$ for real constants $A$ and $B$, where $S^2$ is the two-dimensional sphere..
I have external constraints on these constants which are $-\frac13 < A < \frac23$, $-\frac13 < B < \frac23$ and $-\frac13 < -(A + B) < \frac23$, though I don't believe that anything blows up with this particular function if those constraints are not met -- plotting the integral numerically for different values of $A$ and $B$ gives something that looks smooth.
So far I have attempted to put it in spherical coordinates, and done the common substitution $\mu = \cos\theta$ which yields the following (after using the Pythagorean identity): $$ \mathcal{I} = e^{-B} \int_{-1}^1 d\mu \: e^{-B \mu^2} \int_{0}^{2\pi} d\phi \: e^{-\alpha (1 - \mu^2) \sin^2\phi} $$ with $\alpha = (A - B)$ just another real constant. One may compute either of the two integrals to get: $$ \mathcal{I} = 2\pi e^{-B} \int_{-1}^1 d\mu \: e^{-B \mu^2} I_0\left(\alpha(1 - \mu^2)\right) $$ or $$ \mathcal{I} = \sqrt{\pi} e^{-B} \int_0^{2\pi} d\phi \: \frac{e^{-\alpha \sin\phi} \text{erf}\left(\sqrt{B - \alpha \sin\phi}\right)}{\sqrt{B - \alpha \sin\phi}} $$ where $\text{erf}$ is the error function, and $I_0$ is the zeroth modified Bessel function of the first kind. At this point I am somewhat stuck. Maybe complex integration could do something for the latter expression, but the square roots mean branch cuts which seem like bad news.
For context, I'm trying to write an efficient numerical solver that has to evaluate this function a lot, so any expression that increases the efficiency above numerical quadrature is helpful.