Assume I have a diagonal matrix $L$ of size $n$. I want to compute the following integral: $$I_n(L) \equiv \int_{(\mathbb{S}^{n-1})^2} \mathrm{d}\sigma(x) \mathrm{d}\sigma(x') \exp[n x^\top L \, x']$$ In which $\mathbb{S}^{n-1}$ is the unit sphere in $\mathbb{R}^n$, and $ \mathrm{d}\sigma$ is the usual measure on $\mathbb{S}^{n-1}$, with $\int_{\mathbb{S}^{n-1}} \mathrm{d}\sigma(x) = 2 \pi^{n/2}/\Gamma(n/2)$. I am looking for either a general form of $I_n(L)$ or simply an asymptotic of $\frac{1}{n}\log I_n(L)$ as $n \to \infty$.
I tried but could not find such results. Any help would be appreciated. Thanks !
There is a way to get an infinite series using the expansion of $\exp$, however, it may be more efficient to compute the result using quadrature. The steps below first expand $\exp$, and use the multinomial theorem to get a polynomial expansion for the integrand. Taking the integral inside the sum to the monomials we can evaluate the integral directly to obtain an analytic (closed form) expression for the integral.
For clarity, first denote the integral by $I$, i.e. \begin{align} I := \int_{S_n}\int_{S_n}\exp(n \textbf{x}^T L \textbf{y}) d\sigma_n d\sigma_n \end{align}
Power series for exp
If we first consider the expansion of the function argument in the integral, \begin{align} \exp(n \textbf{x}^T L \textbf{y}) =& \sum_{k=0}^\infty \frac{(n\textbf{x}^T L \textbf{y})^k}{k!} \\ =& \sum_{k=0}^\infty \frac{\left(n \sum_{r=1}^M x_r L_{rr}y_r \right)^k}{k!} \tag{1} \end{align}
It is well known that the series converges uniformly, and the spherical integral is linear (see [1]), which means we can take the integrals inside the summation,
\begin{align} I = \sum_{k=0}^\infty \frac{n^k}{k!} \int_{S_n}\int_{S_n}\left( \sum_{r=1}^M x_r L_{rr}y_r \right)^kd\sigma_n d\sigma_n \tag{2} \end{align}
Use of the multinomial theorem
Using the multinomial theorem on the inner bracket gives, \begin{align} \left( \sum_{r=1}^M x_r L_{rr}y_r \right)^k =& \sum_{j_{1} + \ldots + j_{M} = k} \frac{k!}{j_{1}!\ldots j_{M}!} \prod_{t=1}^M (L_{tt} y_t x_t)^{j_t} \tag{3} \end{align}
We can then put (3) into (2) and distribute the integrals further over the addition to give,
\begin{align} I = \sum_{k=0}^\infty \frac{n^k}{k!}\sum_{j_{1} + \ldots + j_{M} = k} \rho_k \int_{S_n}\int_{S_n} \prod_{t=1}^M (L_{tt} y_t x_t)^{j_{t}}d\sigma_n d\sigma_n \tag{4} \end{align} where $\rho_k = \frac{k!}{j_{1}!\ldots j_{M}!} $.
Integrating a monomial
The spherical integral of a monomial is discussed in detail in [2], but the main result gives, \begin{align} \int_{S_n}\int_{S_n} \prod_{t=1}^M (L_{tt} y_t x_t)^{j_{t}}d\sigma_n d\sigma_n = \left\{ \begin{array}{rl} \left( 2\frac{\prod_{t=1}^M\Gamma(q_{t})}{\Gamma(\sum_{t=1}^M q_t)} \right)^2 \prod_{t=1}^M L_{tt}^{j_{t}} & : j_{t} \mbox{ all even} \\ 0& : \mbox{otherwise} \end{array} \right .\tag{5} \end{align}
Where $q_{t} = \frac{1}{2} (j_{t} + 1)$. Putting (5) into (4) then gives, \begin{align} I =& 4 \sum_{k=0}^\infty \frac{n^k}{k!}\sum_{j_{1} + \ldots + j_{M} = k \\ \forall t \mod(j_{t},2) = 0 } \rho_k \left(\frac{\prod_{t=1}^M\Gamma(q_{t})}{\Gamma(\sum_{t=1}^M q_{t})} \right)^2 \prod_{t=1}^M L_{tt}^{j_{t}} \tag{6} \end{align}
Simplification
Since the terms are non-zero only when all of the the $j_t$, and hence $k$, are even; after some cancellations, the sum can be written, \begin{align} I =& 4 \sum_{k=0}^\infty \sum_{j_{1} + \ldots + j_{M} = k} \frac{n^{2k}}{\prod_{t=1}^M (2j_{t})!}\left(\frac{\prod_{t=1}^M\Gamma(j_{t} + \frac{1}{2})}{\Gamma(\frac{M}{2} + \sum_{t=1}^M j_{t})} \right)^2 \prod_{t=1}^M L_{tt}^{2j_{t}} \\ =& 4 \sum_{k=0}^\infty \frac{n^{2k}}{\Gamma(\frac{M}{2} + k)^2} \sum_{j_{1} + \ldots + j_{M} = k}\prod_{t=1}^M \frac{\left(\Gamma(j_{t} + \frac{1}{2}) L_{tt}^{j_{t}}\right)^2}{ (2j_{t})!}\tag{7} \end{align}
[1] Baker, John A., Integration over spheres and the divergence theorem for balls, Am. Math. Mon. 104, No. 1, 36-47 (1997). ZBL0877.26008.
[2] Folland, Gerald B., How to integrate a polynomial over a sphere, Am. Math. Mon. 108, No. 5, 446-448 (2001). ZBL1046.26503.