Integral over the $\mathcal{S}^{n-1}$ sphere I have been running into the following integral again and again:
Let $S^{n-1}= \{x \in \mathbb{R}^{n} \: | \: ||x||=1 \}$ and let $\lambda_{S^{n-1}}$ denote the surface measure over $S^{n-1}$ as defined in Stroock (2000) page 86. 
Consider a fixed symmetric, positive definite matrix $Q$ of dimension $n \times n$, and a fixed scalar $a\in \mathbb{R}_{+}$
Question 1) Do you know if there is a closed form solution for the integral: 
$$\int_{S^{n-1}} \exp\Big(a \omega'Q\omega \Big) \lambda_{S^{n-1}} (d \omega) $$
When $n=2$, I can express this integral as a modified Bessel function of the first kind $I_{v}(x)$ with $v=0$ evaluated at the eigenvalues of $Q$.
Question 2) Any suggestion about good numerical method for solving this integral? 
Thanks!
*Stroock (2000) "A concise introduction to the theory of integration"
 A: Let $I(a)$, with $a\geqslant0$, denote the integral to be computed. For every $t$, the decomposition of $x$ in $\mathbb R^n$ into a radial part $r$ and a spherical part $\omega$ yields
$$
\int_0^{+\infty}r^{n-1}\mathrm e^{-\frac12tr^2}I\left(\tfrac12r^2\right)\mathrm dr\propto\int_{\mathbb R^n}\mathrm e^{-\frac12x^*(tI-Q)x}\mathrm dx,
$$
where the $\propto$ sign subsumes some proportionality factors depending on $n$.
The RHS is a gaussian integral hence, for every $t$ large enough, is proportional to $(\det(tI-Q))^{-1/2}$. Let $(q_k)_{1\leqslant k\leqslant n}$ denote the eigenvalues of $Q$. The change of variable $s=\frac12r^2$ in the LHS yields
$$
\int_0^{+\infty}\mathrm e^{-ts}I(s)s^{\frac{n-2}2}\mathrm ds\propto\prod_{k=1}^n\frac1{\sqrt{t-q_k}}.
$$
For every fixed $q$, define $g_q:s\mapsto\mathrm e^{qs}s^{-\frac{1}2}$ (hence,
the function $g_q$ is similar the density of a gamma distribution $\Gamma(\frac12,\beta)$ but for some negative $\beta=-q$). Then,
$$
\frac1{\sqrt{t-q}}\propto\int_0^{+\infty}\mathrm e^{-ts}g_{q}(s)\mathrm ds,
$$
hence, by identification,
$$
I(s)\propto s^{-\frac{n-2}2}\cdot h(s),\qquad h=g_{q_1}\ast g_{q_2}\ast\cdots\ast g_{q_n}.
$$
The convolution $h$ has no simple expression in the general case but, when $s\to0$, $h(s)$ is proportional to $s^{\frac{n-2}2}$, hence the prefactor $s^{-\frac{n-2}2}$ in the expression of $I(s)$ cancels out and there is no singularity at $s=0$, as was to be expected since $I(0)=1$. 
This remark and an estimation of $h(s)$ when $s\to0$ show that $I(s)=c_n\cdot s^{-\frac{n-2}2}\cdot h(s)$ with
$$
\frac1{c_n}=\iint_{0\leqslant s_1+\cdots+s_{n-1}\leqslant1}s_1^{-1/2}s_2^{-1/2}\cdots (1-s_1-\cdots-s_{n-1})^{-1/2}\mathrm ds_1\cdots\mathrm ds_{n-1},
$$
that is, $c_n=\Gamma(n/2)\cdot\pi^{-n/2}$. Note finally that, as a confirmation, $I(s)=\mathrm e^{sq}$ in every dimension $n\geqslant1$ when $q_1=\cdots=q_n=q$.
