**Edit:** The below answer addresses a separate question than that in the updated question. @WillieWong's comments above are most relevant to the question when the function depends on the difference $\mathbf{u}-\mathbf{v}$ instead of the distance.

'The function depends only on the distance between the two vectors, and their dot product.'

In other words the function can be taken to only be a function of the dot product, or distance (since $\|x-y\|^2=2-2\langle x,y \rangle$, for $x,y$ on the sphere).

One can expand the function $f$ into its Gegenbauer expansion $f(\langle x,y \rangle)=\sum\limits_{k=0}^{\infty} a_k C_{k}(\langle x,y \rangle)$ where $C_{k}(t)$ are Gegenbauer polynomials. Assuming the convergence is good enough to allow exchanging the sum and integral, the addition formula gives $$\int_{\mathbb{S}^{d-1}}\int_{\mathbb{S}^{d-1}} \sum\limits_{k=0}^{\infty} a_k C_{k}(\langle x,y \rangle) dx dy=\sum\limits_{k=0}^{\infty}\sum\limits_{m=1}^{\dim V_k}a_k b_k \left(\int_{\mathbb{S}^{d-1}} Y_{m,k}(x)dx \right)^2,$$

where spherical harmonics are denoted $Y_{m,k}$, $V_k$ is the finite-dimensional space of spherical harmonics of degree $k$, and $b_k$ are constants (which depend on how one normalizes the spherical harmonics $Y_{m,k}$). This reduces the double integral to a sum of integrals. It is not clear if $dx$ denotes the surface measure on the sphere. If so, then all terms with $k\geq 1$ are zero, meaning the integral is just the constant term in the Gegenbauer expansion of $f$.

difference$\mathbf{u}-\mathbf{v}$, and not just theirdistance$|\mathbf{u} - \mathbf{v}|$, can you clarify? $\endgroup$ – Willie Wong Nov 18 '19 at 15:02anda unit vector $\omega \perp (\mathbf{u}-\mathbf{v}$, you can solve for $$ \mathbf{u} = \sqrt{1 - |\mathbf{u}-\mathbf{v}|^2 / 4} \omega + \frac12 (\mathbf{u}-\mathbf{v}) $$ and similarly $\mathbf{v}$. So you can reparametrize the the product of the two unit spheres (a $2(D-1)$ dimensional manifold) by the product of the... $\endgroup$ – Willie Wong Nov 18 '19 at 17:27