Analytic expression for $\mathbb E_X[(X^\top u)^p (X^\top v)^p]$, where $u$ and $v$ are fixed vectors in $\mathbb R^d$ and $X$ is uniform on sphere Let $u$ and $v$ be fixed vectors in $\mathbb R^d$. Let $X$ be a random vector uniformly distributed on the unit-sphere in $\mathbb R^d$, and let $f_p$ be a real polynomial of degree $p \ge 1$ (for simplicity, we might simply consider $f_p(t) \equiv t^p$).

Question. What is an analytic expression for the correlatetion $c(u,v) := \mathbb E_X[f(X^\top u)f(X^\top v)]$ ?

Special case $p=1$
In this case, a simple computation gives
$$
c(u,v) = \mathbb E[X^\top uv^\top X] = \mbox{trace}(\mbox{cov}(X)uv^\top) = \mbox{trace}((1/d) I_d uv^\top) = \frac{u^\top v}{d}
$$

Is it too crazy to conjecture that in general, $c(u,v) = K_{d,p} \cdot (u^\top v)^p$, for some constant $K_{d,p}$ which only depends on $d$ and $p$ ?

 A: Here's a strategy for computing $c_{p,q}(u,v)=\mathbb{E}[(X^T u)^p(X^T v)^q]$. First compute the generating function $G(\lambda)=\mathbb{E}(e^{i\lambda^TX})$, then set $\lambda=au+bv$. The result is that we can express the quantities above as derivatives of the generating function:
$$c_{p,q}(u,v)=(-i)^{p+q}\frac{\partial^{p+q}}{\partial a^p \partial b^q}G(au+bv)\Bigg|_{a=b=0}$$
Calculating the generating function is fairly simple for $X$ uniformly distributed on the sphere by performing a rotation so that the vector $\lambda$ is in the $x_1$ direction:
$$G(\lambda)=\int e^{i\lambda\cdot x}\frac{d\Omega_{d-1}}{\mu(S_{d-1})}=\frac{1}{\mu(S_{d-1})}\int_{R^{d}}d^dx \delta(|x|-1)e^{i\lambda\cdot x}=\frac{1}{\mu(S_{d-1})}\int_{R^{d}}d^dx \delta(|x|-1)e^{i|\lambda| x_1}$$
Reverting back to spherical coordinates
$$G(\lambda)=\int e^{i|\lambda|\cos\theta}\frac{d\Omega_{d-1}}{\mu(S_{d-1})}=\frac{\mu(S_{d-2})}{\mu(S_{d-1})}\int_{0}^{\pi}\sin^{d-2}\theta~ e^{i|\lambda|\cos\theta}d\theta$$
We recognize the last integral as the Poisson integral for Bessel functions, and hence $G(\lambda)$ can be succinctly written in closed form
$$G(\lambda)=2^{d-2/2}\Gamma\left(\frac{d}{2}\right)\frac{J_{\frac{d-2}{2}}(|\lambda|)}{|\lambda|^{\frac{d-2}{2}}}$$
Using the series representation of the Bessel function one can show that
$$G(au+bv)=\sum_{n=0}^{\infty}\frac{(-1)^n \Gamma(d/2)}{n!\Gamma(n+d/2)2^{2n}}(a^2|u|^2+b^2|v|^2+2ab (u\cdot v))^n$$
We see immediately that the result is non-zero only when $p+q$ is even. It is also obvious that the conjecture cannot hold since
$$c_{2,2}(u,v)\propto 4(u\cdot v)^2+2 |u|^2 |v|^2$$
I have not found a satisfactory algebraic formula for arbitrary $p,q$, to be continued!
EDIT: I came up with a good way to expand this sum.
Obviously when one choose particular values $p+q=2m$, only the term with $n=m$ in the sum will contribute. All we need to do is expand the multivariable polynomial in parentheses above for any value of $n$. The trinomial theorem is not very illuminating so we pick a different route. First rewrite with $z=b/a$
$$(a^2u^2+b^2v^2+2ab(u\cdot v))^m=a^{2m}|v|^{2m}(z+\lambda e^{i\theta})^m(z+\lambda e^{-i\theta})^m$$
where $\lambda=|u|/|v|$, and, $|u||v|\cos\theta=u\cdot v$, $|u||v|\sin\theta=\sqrt{u^2 v^2-(u\cdot v)^2}$
Now we can use the binomial expansion on each of the individual terms
$$(a^2u^2+b^2v^2+2ab(u\cdot v))^m=a^{2m}|u|^{2m}\sum_{kl}{m\choose k}{m\choose l}\left(\frac{z}{\lambda}\right)^{k+l}(e^{i\theta})^{l-k}$$
All we need to do now is reindex the sum so that everything is written as a coefficient of $z^s$:
$$(a^2u^2+b^2v^2+2ab(u\cdot v))^m=a^{2m}|u|^{2m}\sum_{s=0}^{2m}z^s(\lambda e^{i\theta})^{-s}\sum_{l=\max(s-m,0)}^{\min(s,m)}{m\choose s-l}{m\choose l}e^{i2l\theta}$$
One can explicitly show that this expression is real. We take the real part of the right hand side and since $\cos(n\theta)=T_n(\cos\theta)$ are the Chebyshev polynomials
$$(a^2u^2+b^2v^2+2ab(u\cdot v))^m=a^{2m}|u|^{2m}\sum_{s=0}^{2m}z^s\lambda^{-s}\sum_{l=\max(s-m,0)}^{\min(s,m)}{m\choose s-l}{m\choose l}T_{2l-s}\left(\frac{u\cdot v }{|u||v|}\right)$$
with the understanding that $T_{-n}(x)=T_n(x), n>0$. Finally, we conclude that
$$c_{2m-p,p}(u,v)=\frac{\Gamma(d/2)\Gamma(p+1)\Gamma(2m-p+1)}{2^{2m}\Gamma(m+d/2)\Gamma(m+1)}|u|^{2m-p}|v|^p\sum_{l=\max(p-m,0)}^{\min(p,m)}{m\choose p-l}{m\choose l}T_{2l-p}\left(\frac{u\cdot v }{|u||v|}\right)$$
I haven't attempted to try and simplify the Chebyshev polynomials much further. Further simplification occurs when $m=p$ but not to the point of a compact expression better than
$$c_{p,p}=\frac{\Gamma(d/2)\Gamma(p+1)}{2^{2p}\Gamma(p+d/2)}|u|^p |v|^p\sum_{l=0}^{p}{p \choose l}^2T_{2l-p}\left(\frac{u\cdot v }{|u||v|}\right)$$
A: Let $Y=(Y_1,\ldots,Y_d)$ be a collection of iid $N(0,1)$'s. Then by using (hyper)spherical coordinates
$$
E[(Y\cdot u)^p(Y\cdot v)^q]=(2\pi)^{-\frac{d}{2}}\int_{0}^{\infty}dr\ {\rm vol}(S^{d-1}) \ r^{d-1+p+q}\ e^{-\frac{r^2}{2}}
\ E[(X\cdot u)^p(X\cdot v)^q]
$$
$$
=\frac{2^{\frac{p+q}{2}}\Gamma\left(\frac{p+q+d}{2}\right)}{\Gamma\left(
\frac{d}{2}\right)}\times
\ E[(X\cdot u)^p(X\cdot v)^q]
$$
We need $p+q$ even for a nonzero result, so we assume that. The Gaussian expectation
$E[(Y\cdot u)^p(Y\cdot v)^q]$ can be computed with the Isserlis-Wick Theorem. It amounts to a sum over complete matchings of a set with $p+q$ elements subdivided into two compartments, one of size $p$ and one of size $q$. Let me organize the count according to $k$, the number of matched pairs joining the two compartments. We then get
$$
E[(Y\cdot u)^p(Y\cdot v)^q]=
$$
$$
\sum_{k}\binom{p}{k}\binom{q}{k}
\times k!\times\frac{(p-k)!}{2^{\frac{p-k}{2}}\left(\frac{p-k}{2}\right)!}
\times\frac{(q-k)!}{2^{\frac{q-k}{2}}\left(\frac{q-k}{2}\right)!}
\times (u\cdot v)^{k}(u\cdot u)^{p-k}(v\cdot v)^{q-k}
$$
where the sum is over $0\le k\le\min(p,q)$ of same parity as $p$ and $q$.
Finally, after some cleanup,
$$
E[(X\cdot u)^p(X\cdot v)^q]=
$$
$$
\frac{p!q!\ \Gamma\left(
\frac{d}{2}\right)}{2^{p+q}\ \Gamma\left(
\frac{p+q+d}{2}\right)}\times \sum_k
\frac{2^k}{k!\left(\frac{p-k}{2}\right)!\left(\frac{q-k}{2}\right)!}\times
(u\cdot v)^{k}(u\cdot u)^{p-k}(v\cdot v)^{q-k}
$$
with the same range of summation for $k$.
(this is a quick computation, didn't yet check for errors typos etc.)
A: Here's an alternate solution I came up with. So, let $Z$ be a random vector in $\mathbb R^d$ with iid components distributed according to $N(0,1)$. Set $R := \|X\|$ and note that

*

*$X$ and $Z/R$ have the same distribution.

*$R$ and $Z$ are independent.

Let $A := (uv^\top + vu^\top)/2$, a psd matrix, and note that $E[(X^\top u)^p(X^\top v)^p] = E[(X^\top A X)^p]$. One computes
$$
E[(Z^\top A Z)^p] = E[R^{2p}]\cdot E[(X^\top AX)] = E[R^{2p}]\cdot E[(X^\top u)^p(X^\top v)^p]. \tag{1}
$$
Now, $E[R^{2p}] = E[(\|Z\|^2)^p] = 2^p \cdot \dfrac{\Gamma(p + d/2)}{\Gamma(d/2)}$, the $p$th moment of a chi-squared random variable with $d$ degrees of freedom It remains to compute. Combining with (1) then gives
$$
E((X^\top u)^p(X^\top v)^p] = \frac{\Gamma(d/2)}{2^p\Gamma(p+d/2)}\cdot\sigma_p(A).\tag{2}
$$
It remains to compute $\alpha_p(A) := E[(Z^\top A Z)^p]$.
Now, since as early as this paper https://www.janmagnus.nl/papers/JRM003.pdf, it is known that the $\sigma_p(A)$ is a polynomial in the first $p$ cumulants $\kappa_1,\ldots,\kappa_p$ of $N(0,1)$, and in fact for the first $2$ values of $p$, evaluates to
$$
\begin{split}
\sigma_1(A) &= \mbox{tr} A = u^\top v\\
\sigma_2(A) &= (\mbox{tr} A)^2 + 2\mbox{tr} A^2 = 2(u^\top v)^2 + \|u\|^2 \|v\|^2\\
&\;\;\vdots
\end{split}
$$
In particular, using (2) gives $E[(X^\top u)(X^\top v)] = \dfrac{u^\top v}{d}$ (which was already computed directly) and
$$
E[(X^\top u)^2(X^\top v)^2] = \dfrac{2(u^\top v)^2 + \|u\|^2 \|v\|^2}{4(d/2 + 1)(d/2)} = \dfrac{2(u^\top v)^2 + \|u\|^2\|v\|^2}{d(d+2)}.
$$
