Almost independence of $x^\top a$ and $x^\top b$ for $x$ uniform on the sphere in $\mathbb R^d$ and $a,b \in \mathbb R^d$ with $a^\top b = 0$ Let $d$ be a large positive integer. Let $x$ be uniformly distributed on the unit-sphere in $\mathbb R^d$ and let $a$ and $b$ be perpendicular vectors in $\mathbb R^d$, i.e such that $a^\top b=0$. Let $f:\mathbb R \to \mathbb R$ be a continuous function which is twice-continuously differentiable at $0$ (you may assume more regularity if that helps).

Question. How large can the absolute difference $$
\Delta(a,b)=|\mathbb E[f(x^\top a)f(x^\top b)] - \mathbb E[f(x^\top a)]\mathbb E[f(x^\top b)]|
$$
be as a function of the dimension $d$ ?

Observation 1. In the case where $x \sim N(0,I_d)$, or more generally, with any diagonal covariance matrix, we have $\Delta(a,b) = 0$. Indeed, in this case, the joint distribution of $(x^\top a,x^\top b)$ is $N(0, C)$, where $C$ is the $2 \times 2$ psd matrix with entries given by

*

*$c_{11}=\mathbb E[(x^\top a)^2] = \mbox{trace}(aa^\top)=\|a\|^2$,

*$c_{12}=c_{21}=\mathbb E[(x^\top a)(x^\top b)] = \mbox{trace}(ab^\top) = a^\top b = 0$,

*$c_{22}=\mathbb E[(x^\top b)^2] = \mbox{trace}(bb^\top)=\|b\|^2$.

Thus, $x^\top a$ and $x^\top b$ are jointly gaussian random variables with zero correlation, and so are independent. Thus $\Delta(a,b) = 0$.
Observation 2. Since $d$ is large both $x^\top a$ and $x^\top b$ are close to centered normal random variables with variance $\|a\|^2/d$ and $\|b\|^2/d$ respectively, so maybe an adaptation of the previous argument might be of use ?
 A: $\newcommand{\de}{\delta}\newcommand{\ep}{\varepsilon}$The covariance of two random variables (r.v.'s) does not change if one of them is shifted by a constant. So, without loss of generality $f(0)=0$. Let $n:=d$.
To compute the asymptotics, we need to assume that
\begin{equation*}
    f(x)=Ax+Bx^2+Cx^3+Dx^4+e_1x^5+O(x^6) \tag{1}
\end{equation*}
for some real $A,B,C,D,e_1$ and all real $x\in[-1,1]$.
The joint distribution of $x^\top a$ and $x^\top b$ is the same as that of $Y_1$ and $Y_2$, where
\begin{equation*}
    Y_j:=X_j/|X|,
\end{equation*}
$X_1,\dots,X_n$ are iid standard normal r.v.'s and $|X|:=\sqrt{X_1^2+\dots+X_n^2}$. So, we want to find the asymptotics of
\begin{equation*}
    \de:=Ef(Y_1)f(Y_2)-(Ef(Y_1))^2 \tag{2}
\end{equation*}
(as $n\to\infty$).
Note that
\begin{equation*}
    EY_1^2=\dots=EY_n^2=\frac1n, \tag{3}
\end{equation*}
since the $Y_i$'s are exchangeable and $Y_1^2+\dots+Y_n^2=1$. Also, using e.g. a Chernoff exponential concentration inequality for $|X|$, we see that $|X|$ is highly concentrated near $\sqrt n$, say in the sense that, for each real $\ep>0$,
$P(|\,|X|-\sqrt n|>\ep\sqrt n)$ goes to $0$ faster than any negative power of $n$. It follows that for each real $p>0$
\begin{equation*}
    E|Y_1|^p\sim E|X_1|^p n^{-p/2}. \tag{4}
\end{equation*}
So, by (1) and the symmetry of (the distribution of) $Y_1$,
\begin{equation*}
    Ef(Y_1)=\frac Bn+\frac{(3+o(1))D}{n^2}.  \tag{5}
\end{equation*}
Next, by (1), the symmetry of the $Y_i$'s, and (4)
\begin{equation*}
    Ef(Y_1)f(Y_1)\\
    =B^2\,EY_1^2Y_2^2+2BD\,EY_1^2Y_2^4+o(n^{-3}).  \tag{6}
\end{equation*}
Similarly to (4),
\begin{equation*}
    EY_1^2Y_2^4\sim EX_1^2X_2^4 n^{-3}=\frac3{n^3}. \tag{7}
\end{equation*}
The main difficulty here is to estimate $EY_1^2Y_2^2$. We have
\begin{equation}
    EY_1^2Y_2^2=Eh(G), \tag{8}
\end{equation}
where $G^2$ has the $\chi^2$ distribution with $n-2$ degrees of freedom and
\begin{align*}
    h(x)&:=E\frac{X_1^2X_2^2}{(X_1^2+X_2^2+x^2)^2} \\ 
    &=\frac1{2\pi}\int_0^{2\pi}dt\int_0^\infty r\,dr\,e^{-r^2/2}\frac{r^4\cos^2t\sin^2t}{(r^2+x^2)^2} \\ 
    &=\frac1{16}\int_0^\infty du\,\frac{u^2e^{-u/2}}{(u+x^2)^2} \\ 
    &=\frac{x^2}{16}\int_0^\infty dt\,\frac{t^2e^{-x^2 t/2}}{(1+t)^2}. 
\end{align*}
Writing $1/(1+t)^2=1-2t+O(t^2)$ for $t\in(0,1)$, we now get
\begin{equation}
    h(x)=\frac1{16}\,(x^{-4}-(12+o(1))x^{-6}) \tag{9}
\end{equation}
as $x\to\infty$.
Also,
\begin{equation}
    EG^{-4}=\frac1{n^2}+\frac{10+o(1)}{n^3},\quad EG^{-6}=\frac{1+o(1)}{n^3}.
\end{equation}
Hence, in view of (9),
\begin{equation}
    EY_1^2Y_2^2=Eh(G)=\frac1{n^2}-\frac{2+o(1)}{n^3}. 
\end{equation}
Finally, recalling (2), (5), and (7), we get
\begin{equation}
    \de=-\frac{2B^2+o(1)}{n^3}. 
\end{equation}
So, as should be expected, we have a small negative correlation between $x^\top a$ and $x^\top b$ (on the order of $1/n^2$, in view of (3)).
A: It turns out that one can get a stronger result than demanded in the question: compute $\Delta(a,b)$ for any $a,b \in S_{d-1}$, perpendicular or not. Indeed,

Claim. If $f(\rho)=a_0 + a_1 \rho + a_2 \rho^2 + a_3 \rho^3 + \mathcal O(\rho^4)$, then for every $u,v \in S_{d-1}$, and $x$ be uniform on the sphere, then we have the following approximation
\begin{eqnarray}
\begin{split}
\mathbb E_x[f(x^\top u)f(x^\top v)]-\mathbb E_x[f(x^\top u)]\mathbb E_x[f(x^\top v)] &= -\frac{a_2^2}{d^2}\\
&+(\dfrac{a_1^2}{d}+\frac{6a_1a_3}{d^2})u^\top v\\
& + \dfrac{2a_1^2}{d^2}(u^\top v)^2\\
&+\mathcal O(\dfrac{1}{d^3}).
\end{split}
\end{eqnarray}
In paricular, if $u^\top v = 0$, then
$$
\mathbb E_x[f(x^\top u)f(x^\top v)]-\mathbb E_x[f(x^\top u)]\mathbb E_x[f(x^\top v)] = -\frac{a_2^2}{d^2}+\mathcal O(\frac{1}{d^3}),
$$
as proved by Iosif Pinelis (in the accepted answer).

For the proof, we will need the following result

Lemma.
Let $u$ and $v$ be fixed and $x$ be uniformly random on the unit-sphere. Let $p$ and $q$ be nonnegative integers and define $c_{p,q}(u,v):=\mathbb E_x[(x^\top u)^p(x^\top v)^q]$. If $p$ and $q$ have different parities, then $c_{p,q}(u,v)=0$. Otherwise, we have the formula
\begin{eqnarray}
c_{p,q}(u,v) = \dfrac{p!q!\Gamma(\frac{d}{2})}{2^{p+q}\Gamma(\frac{d+p+q}{2})} \sum_t \dfrac{2^t}{t!(\frac{p-t}{2})!(\frac{q-t}{2})!}(u^\top v)^t,
\end{eqnarray}
where the sum is over all $t$ between $0$ and $p \land q$ inclusive, that have the same parity as $p$ and $q$. The formula is simplified in the table below for special values of $p$ and $q$.
The above lemma is proved (ME link) here https://math.stackexchange.com/a/4004804/168758.

Proof of the claim.
WLOG, let $a_0=0$. Thanks to the lemma, one may compute
$$
(\mathbb E_x[f(x^\top u)])^2 = (\frac{a_2}{d}+\mathcal O(\frac{1}{d^2}))^2 = \frac{a_2^2}{d^2} + \mathcal O(\frac{1}{d^3}),
$$
and similarly
\begin{eqnarray*}
\begin{split}
    \mathbb E_x[f(x^\top u)f(x^\top v)] &= \sum_{p=0}^3\sum_{q=0}^3 a_p a_q\mathbb E_x[(x^\top u)^p(x^\top v)^q]+\mathcal O(\dfrac{1}{d^3})\\
    &= \dfrac{a_1^2}{d}u^\top v + \dfrac{2a_2^2}{d(d+2)}(u^\top v)^2 + \dfrac{6 a_1a_3}{d(d+2)}u^\top v+\mathcal O(\dfrac{1}{d^3})
\\ &=(\dfrac{a_1^2}{d}+\frac{6a_1a_3}{d^2})u^\top v + \dfrac{2a_2^2}{d^2}(u^\top v)^2 +\mathcal O(\dfrac{1}{d^3}),
\end{split}
\end{eqnarray*}
and the claim follows after subtracting the previous display. $\quad\quad\quad\quad\Box$
