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$ ?
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)$$