Let $f:\mathbb R \to \mathbb R$ be a continuous function and define $F:S_{d-1} \times S_{d-1} \to \mathbb R$ by $F(u,v) := \mathbb E[f(X^\top u)f(X^\top v)]$, where $S_{d-1}$ is the unit-sphere in $\mathbb R^d$ and $X \sim N(0,I_d)$. Because of rotational invariance of the distribution of $X$, it is clear that $F(u,v) \equiv h(u^\top v)$ for some continuous function $h:[-1,1] \to \mathbb R$.
Questions:
- (1) How smooth can the function $h$ be on the open interval $(-1,1)$, especially around the point $0$ ?
- (2) In case $h$ is differentiable at $0$, what is the value of $h'(0)$ in terms of $f$ ?
Observation
For computing the value of $h(0)$, one can take any $a,b \in S_{d-1}$ such that $a^\top b = 0$. Then $X^\top a$ and $X^\top b$ are independent random variables and so $$ h(0) = h(a^\top b) = \mathbb E[f(X^\top a)f(X^\top b)] = \mathbb E[f(X^\top a)]\mathbb E[f(X^\top b)] = \mathbb E[f(X_1)]^2 = \sigma_0(f)^2, $$ where $\sigma_k(f)$ is the $k$th Hermite coefficient of $f$, defined by $\sigma_k(f) = \mathbb E_{G \sim N(0,1)}[f(G)\operatorname{He}_k(G)]$ and $\operatorname{He}_k$ is the probabilist's $k$th Hermite polynomial; in particular, $\mathrm{He}_0(z) = 1$, $\mathrm{He}_1(z) := z$, $\mathrm{He}_2(z) := z^2-1$, $\mathrm{He}_3(z) := z^3-3z^2$, etc.
Similarly, for computing $h(1)$ one can take $a \in S_{d-1}$ and note that $$ h(1)=h(a^\top a) = \mathbb E[f(X^\top a)^2] = \mathbb E[f(X^\top a)^2] = \mathbb E[f(X_1)^2]=\sigma_0(f^2). $$
Let $Z_1,Z_2\sim\text{i.i.d.}\operatorname N(0,1)$ and let $u\in[-1,1].$ Let $Z_3= uZ_1 + \sqrt{1-u^2}\cdot Z_2.$
Then the joint distribution of $(Z_1,Z_3)$ is the same as the distribution of the thing that in your notation would be called $(X^\top u, X^\top v)$ when the number I have called $u$ is the same as the thing that you called $u^\top v.$ Then, $(Z_1,Z_3)$ is jointly Gaussian with zero mean, covariance matrix $C=\begin{bmatrix}1 & u\\u & 1\end{bmatrix}$, and joint pdf $$ \begin{split} p_{Z_1,Z_3}(x,y) &= \frac{1}{(2\pi)^{2/2} \det(C)^{1/2}} \exp\left(-\frac{1}{2}(x,y) C^{-1} (x,y)^\top \right) \\ &= \frac 1 {2\pi\sqrt{1-u^2}} \cdot \exp\left( -\frac 1 {2(1-u^2)} (x^2 -2uxy+y^2) \right). \end{split} $$ Therefore \begin{align} & h(u) = \operatorname E\big(f(Z_1) f(Z_3)\big) \\[8pt] = {} & \int\limits_{\mathbb R^2} \frac{f(x)f(y)}{2\pi\sqrt{1-u^2}}\cdot \exp\left( -\frac 1 {2(1-u^2)} (x^2-2uxy+y^2) \right) \, d(x,y). \end{align} Here one should recall the circumstances under which one can differentiate under the integral sign.
Or one can find $\dfrac{h(u) - h(0)} u$ and let $u\to0,$ and remember what you found $h(0)$ to be.