Reference Request: Integral of Gaussian over Unit Sphere I am looking for a reference for integrals of the form
\begin{equation}
\tag{1} \int_{S^{n-1}} \mathcal{N}_{\omega} ( \mu , \Sigma ) d \omega
\end{equation}
where $S^{n-1}$ is the sphere in $\mathbb{R}^n$ and
\begin{equation}
\mathcal{N}_x (\mu , \Sigma) = \frac{e^{- \frac{1}{2} ( x - \mu ) \cdot \Sigma^{-1} (x - \mu) }}{\sqrt{ ( 2 \pi )^n \det{\Sigma} }}
\end{equation}
is the usual multivariate Gaussian distribution evaluated at $x \in \mathbb{R}^n$ with mean $\mu \in \mathbb{R}^n$ and covariance matrix $\Sigma \in \mathbb{R}^{n \times n}$. I am primarily concerned with the $n = 2, 3$ cases, but I find the case for general $n$ to be interesting as well.
There is a substantial simplification that can be made. By the spectral theorem and positive-definite nature of $\Sigma$ we may write $\Sigma = V D^2 V^T$ where $V \in \mathrm{SO}_n (\mathbb{R})$ and $D = \mathrm{diag}(\sigma_1 , \cdots , \sigma_n)$ is a diagonal matrix with positive, real values along the diagonal. Therefore, by spherical symmetry, I am essentially asking for a reference to the integral
\begin{equation}
\tag{2} \int_{S^{n-1}} e^{ - \frac{1}{2} \| D^{-1} (\omega - \tilde{\mu}) \|^2 } d \omega ,
\end{equation}
where $\tilde{\mu} = V^T \mu$.
Question: Does anyone have a reference or suggestion for how to calculate the closed-form expression for (2) and, hence, (1)?
Additional Thoughts:
There are special cases that are simple enough to compute. For example, if $\sigma_1 = \cdots = \sigma_n \equiv \sigma$, then the integral in question becomes
\begin{equation}
\frac{e^{- \frac{1 + |\mu|^2}{2 \sigma^2} }}{(2 \pi)^{n/2} \sigma} \int_{S^{n-1}} e^{\omega \cdot \mu / \sigma^2 }  d \omega . 
\end{equation}
The previous integral may be evaluated using the formula
\begin{equation}
\int_{S^{n-1}} e^{ \omega \cdot y } d \omega = (2 \pi)^{\frac{n}{2}} | y |^{1 - \frac{2}{n}} I_{\frac{n}{2} - 1} ( |y| ) ,
\end{equation}
where $I_{\nu}$ is the hyperbolic Bessel function of order $\nu$. This formula may be found in, e.g., Loss and Lieb, Analysis, Section 7.11. Therefore,
\begin{equation}
\tag{3} \int_{S^{n-1}} \mathcal{N}_{\omega} ( \mu , \sigma^2 ) d \omega = \frac{e^{- (1 + \mu^2)/(2 \sigma^2) }}{\sigma} \left( \frac{|\mu|}{\sigma^2} \right)^{1 - \frac{2}{n}} I_{\frac{n}{2} - 1} \left( \frac{|\mu|}{\sigma^2}  \right) .
\end{equation}
Another special case that I mention in the comments below is $\mu = 0$. As I point out there, the $n = 3$ case of (1) for $\mu = 0$ should be similar to the calculation of the normalization constant of the Kent distribution. For $n = 2$, (2) reduces to
\begin{equation}
\tag{4} \int_0^{2 \pi} e^{- \frac{1}{2} \left( \frac{\cos^2{\theta}}{\sigma_1^2} + \frac{\sin^2{\theta}}{\sigma_2^2} \right) } d \theta = 2 \pi e^{- \frac{1}{4} \left( \frac{1}{\sigma_1^2} + \frac{1}{\sigma_2^2} \right) } I_0 \left( \frac{1}{4} \left( \frac{1}{\sigma_2^2} - \frac{1}{\sigma_1^2} \right) \right) .
\end{equation}
The previous result follows from the cosine and sine reduction formula $\cos^2{x} = (1 + \cos{(2x)})/2$ and likewise for $\sin^2{x}$.
Considering that (3) and (4) have relatively simple expressions, I am hoping this is also the case for (2) with $\tilde{\mu} \neq 0$.
 A: For the $n = 2$ case, (2) in the OP can be evaluated as follows. Denote $\tilde{\mu} = \mu$. The integral in question may be reduced to
\begin{equation}
\tag{1} e^{- \frac{1}{2} \left( \frac{\mu_1^2}{\sigma_1^2} + \frac{\mu_2^2}{\sigma_2^2} \right) } \int_0^{2 \pi} e^{ - \frac{1}{2} \left( \frac{\cos^2{\theta}}{\sigma_1^2} + \frac{\sin^2{\theta}}{\sigma_2^2} - 2 \frac{\mu_1}{\sigma_1^2} \cos{\theta} - 2 \frac{\mu_2}{\sigma_2^2} \sin{\theta} \right) } d \theta .
\end{equation}
Using the double angle reduction formulas $\cos^2{\theta} = ( 1 + \cos{(2 \theta)} ) / 2$ and $\sin^2{\theta} = ( 1 - \cos{(2 \theta)} ) / 2$, we can simplify (1) to read
\begin{equation}
\tag{2} e^{- \frac{1}{4} \left( \frac{2\mu_1^2 + 1}{\sigma_1^2} + \frac{2\mu_2^2 + 1}{\sigma_2^2} \right) } \int_0^{2 \pi} e^{ - \frac{1}{2} \left( 
\left( \frac{1}{2 \sigma_1^2} - \frac{1}{2 \sigma_2^2} \right) \cos{(2 \theta)} - 2 \frac{\mu_1}{\sigma_1^2} \cos{\theta} - 2 \frac{\mu_2}{\sigma_2^2} \sin{\theta} \right) } d \theta .
\end{equation}
By the Jacobi-Anger expansion, we have
\begin{equation}
e^{z \cos{2 \theta}} = I_0 (z) + 2 \sum_{k \geq 1} I_k (z) \cos{(2 k\theta)} ,
\end{equation}
where $z = - \frac{1}{4} \left( \frac{1}{\sigma_1^2} - \frac{1}{\sigma_2^2} \right)$. Plugging this into (2), we have
\begin{equation}
\tag{3} e^{- \frac{1}{4} \left( \frac{2\mu_1^2 + 1}{\sigma_1^2} + \frac{2\mu_2^2 + 1}{\sigma_2^2} \right) } \int_0^{2\pi} \left( I_0 (z) + 2 \sum_{k \geq 1} I_k (z) \cos{(2 k \theta)} \right) e^{\frac{\mu_1}{\sigma_1^2} \cos{\theta} + \frac{\mu_2}{\sigma_2^2} \sin{\theta} } d \theta ,
\end{equation}
where, again, $z = - \frac{1}{4} \left( \frac{1}{\sigma_1^2} - \frac{1}{\sigma_2^2} \right)$. To proceed, note that
\begin{equation}
y_1 \cos{\theta} + y_2 \sin{\theta} = \| y \| \cos{(\theta - \arctan{(y_2/y_1)})} ,
\end{equation}
where $y_1 = \mu_1 / \sigma_1^2$ and $y_2 = \mu_2 / \sigma_2^2$. Plugging this into (3), changing variables, and using the cosine addition formula, we arrive at
\begin{equation}
\int_{S^1} e^{- \frac{1}{2} \| D^{-1} ( \omega - \mu ) \|^2} d \omega = 2 \pi e^{- \frac{1}{4} \left( \frac{2\mu_1^2 + 1}{\sigma_1^2} + \frac{2\mu_2^2 + 1}{\sigma_2^2} \right) } \left( I_0 (z) I_0 (\|y\|) + 2 \sum_{k \geq 1} I_k (z) I_{2k} (\| y \|) \cos{\left(2 k \arctan{\left( \frac{y_2}{y_1} \right)}\right)} \right) , \tag{*} 
\end{equation}
where, again, $z = - \frac{1}{4} \left( \frac{1}{\sigma_1^2} - \frac{1}{\sigma_2^2} \right)$, $y_1 = \mu_1 / \sigma_1^2$ and $y_2 = \mu_2 / \sigma_2^2$.
We note how ($\ast$) agrees with some of the special cases outlined above. Indeed, for $\sigma_1 = \sigma_2$, $z = 0$ and all higher order Bessel functions in the sum in ($\ast$) are zero. When $\mu_1 = \mu_2 = 0$, the $\arctan$ isn't defined, so this special case isn't covered.
