Solution (source code)

= Solution

Reflection symmetry makes the <integral> zero for $i\ne j$. Rotational symmetry makes all diagonal <integrals> equal, and their sum is $R^2A_{n-1}(R)$. Therefore
$$
\boxed{\int_{r=R}x_ix_j\,dS=\delta_{ij}\frac{R^2A_{n-1}(R)}n
=\delta_{ij}\frac{\pi^{n/2}R^{n+1}}{\Gamma(n/2+1)}.}
$$