Suppose we have an $n$-dimensional sphere $S^{n-1}$, defined as $S^{n-1}:=\{x\in \mathbb{R}^n: |x|=1\}$. Let $f:S^{n-1}\to \mathbb{R}{\geq 0}$ be the probability density function of $\frac{X}{|X|}$, where $X\sim N(0,\Sigma)$ is a multivariate Gaussian random variable. Since the hyperspherical harmonics form a basis for the functions on the $n$-sphere, we can expand $f$ in terms of the harmonics in higher dimensions $Y{k_1,\dots, k_{n-1}}(\theta_1,\dots,\theta_{n-1},\phi)$, where $\theta_1,\dots,\theta_{n-1}$ are the angular variables that parameterize the $(n-1)$-dimensional hypersphere $S^{n-1}.$
However, it seems that analytically calculating the coefficients of the expansion can be very challenging when $n>3$. I was wondering if the coefficients can be obtained by exploiting the rich properties of the Gaussian weight kernel?