Given $\mathbb E[f(X^\top u)f(X^\top v)] \equiv h(u^\top v)$ for $X \sim N(0,I)$ and unit-vectors $u,v$, compute $h(1)$ and $h'(0)$ in terms of $f$ 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).
$$
 A: 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.
A: Let's prove the following explicit formula

Claim. If $f \in L^2(\mathbb R,N(0,1))$, then we have
$$
h(0) = \sigma_0(f)^2,\,h'(0) = \sigma_1(f)^2,
$$
where $\sigma_0(f):=\mathbb E_{G \sim N(0,1)}[f(G)]$ and $\sigma_1(f) := \mathbb E_{G \sim N(0,1)}[Gf(G)]$ are respectively the order-$0$ and order-$1$ Hermite coefficients of $f$.

Indeed, from Michael Hardy's post, we know that
$$
h(u) = \int_{-\infty}^\infty\int_{-\infty}^\infty f(x)f(y)R(x,y,u)dxdy, \tag{*}
$$
where $R(x,y,u):= \dfrac{1}{2\pi\sqrt{1-u^2}}\exp(-\dfrac{1}{2(1-u^2)}(x^2-2uxy+y^2))$. As a sanity-check, we use (*) to compute the value of $h(0)$ as
$$
h(0)=\int_{-\infty}^\infty\int_{-\infty}^\infty f(x)f(y)\exp(-(x^2+y^2)/2)/(2\pi)dxdy=\left(\int_{-\infty}^\infty f(x)\frac{\exp(-x^2/2)}{\sqrt{2\pi}}dx\right)^2 =\sigma_0(f)^2,
$$
which is the same value obtained as an "observation" in the question.
Let us continue to compute $h'(0)$, our main interest. Using (*), we deduce from Reynold's transport theorem that
$$
\begin{split}
h'(u) &= \int_{-\infty}^\infty\int_{-\infty}^\infty f(x)f(y)\frac{\partial R(x,y,u)}{\partial u}dxdy\\
&=\int_{-\infty}^\infty\int_{-\infty}^\infty f(x)f(y)\frac{\exp(-\frac{(x^2-2uxy+y^2)}{2(1-u^2)})(-u^3+u^2xy-u(x^2+y^2-1)+xy)}{2\pi(1-u^2)^{5/2}}dxdy,
\end{split}
$$
which makes sense for all $u \in (-1,1)$. Repeating this argument, we deduce that $h$ is $\mathcal C^\infty$ on $(-1,1)$. In particular, we compute
$$
\begin{split}
h'(0) = \int_{-\infty}^\infty\int_{-\infty}^\infty xyf(x)f(y)\frac{\exp(-(x^2+y^2)/2)}{2\pi}dxdy = \left(\int_{-\infty}^\infty xf(x)\frac{\exp(-x^2/2)}{\sqrt{2\pi}}dx\right)^2
=\sigma_1(f)^2,
\end{split}
$$
as claimed.
