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"

-
Just a quick, half baked idea: Due to rotational symmetry, you may as well assume that $Q$ is diagonal. Now it seems to me that you can do one dimension at a time, getting a single integral with an integrand involving the answer for $S^{n-2}$. Iterate. – Harald Hanche-Olsen Aug 3 '12 at 13:08
what is $\omega'$? – Mercy King Aug 3 '12 at 13:08
I guess $\omega'$ means $\omega^t$, the transposed vector of $\omega$. – Siminore Aug 3 '12 at 13:20
Yes, sorry: $\omega'$ is the transpose of $\omega$. – JLp Aug 3 '12 at 13:40
Thanks to HHO for the half baked idea. Let me work that out. – JLp Aug 3 '12 at 13:41

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$.