$\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_1,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
$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)).

$\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)).