Hello it has been many years since this comment was posted, however there was some renewed interest on pages that have lead to a resolution of this issue.
[1] Averaged dot product
[2] Surface integral of normal components summations on a sphere
I'll demonstrate how to generalize the proof given on page [2], in order to solve this problem. I have based my answer off Greg's work in [2], because I wanted to think about spherical integration in $D$ dimensional space. It seems easier to conceptualize. However, I will note that the proof given in [1] would require no logical adjustments to yield your desired result.
(Also my thanks goes to @JeanMarie for pointing out this page in his answer to [1])
----------------------------Relation 1-------------------------------------------
So to start, we have by definition that, the expected value is the normalized surface integral.
$$ \left< n_i n_j \right> = \frac{1}{S_D} \int \hat n \hat n dA \tag{*}$$
Where $S_D$ is the surface area of the sphere in $D$ dimensions.
Due to arguments of symmetry, which I have enumerated on [1] and Greg has enumerated on [2] we have that we have an isotropic second order tensor,
$$ \oint \hat n \hat n dA = \beta \delta_{ij} \tag{**}$$
The scale factor can be determined by contraction in $D$ dimensional space.
$$ \delta_{ij} \delta_{ij} = \delta_{ii} = \text{tr}\{ I_{D \times D}\} = D$$
Thus for the normalization condition reads
$$ D\beta = \oint \hat n \cdot \hat n dA = \oint dA = S_D $$
Which implies that,
$$ \beta = \frac{S_D}{D} \tag{***}$$
Inserting expressions ($***$) and ($**$) into ($*$) we arrive at the desired,
$$ \left< n_i n_j \right> = \frac{1}{S_D} \int \hat n \hat n dA = \frac{\delta_{ij}}{D}$$
----------------------------Relation 2-------------------------------------------
We can do something similar for the next relation. We have as before,
$$ \left< n_i n_j n_k n_l \right> = \frac{1}{S_D} \int \hat n \hat n \hat n \hat n dA \tag{x}$$
By arguments of symmetry provided in both [1] and [2],
$$\oint \hat{n}\,\hat{n}\,\hat{n}\,\hat{n}\,dA = \beta\,{\mathbb Y}$$
where ${\mathbb Y}$ is a 4th order isotropic tensor and given as,
$${\mathbb Y}_{ijkl} = \delta_{ij}\,\delta_{kl} + \delta_{ik}\,\delta_{jl} + \delta_{il}\,\delta_{jk} \tag{xx}$$
From here you can follow Greg's proof, so that the contraction $\big\{I:{\mathbb Y}:I\big\}$ yields
$$\eqalign{
\delta_{ij}\,(\delta_{ij}\,\delta_{kl} + \delta_{ik}\,\delta_{jl} + \delta_{il}\,\delta_{jk})\,\delta_{kl}
&= (D\,\delta_{kl} + \delta_{kl} + \delta_{kl})\,\delta_{kl} \cr
&= (D*D + D + D)\cr
&= D^2+2D \cr
}$$
The contraction $\big\{I:(\hat{n}\,\hat{n}\,\hat{n}\,\hat{n}):I\big\}$ yields
$$\eqalign{
(\hat{n}\cdot\hat{n})\,(\hat{n}\cdot\hat{n}) &= (1)\,(1) = 1
}$$
So
$$\eqalign{
(D^2+2D)\,\beta &= \oint dA = S_D \cr
\beta &= \frac{S_D}{D^2+2D} \tag{xxx}
}$$
Inserting expressions ($xxx$) and ($xx$) into ($x$) we arrive at the desired,
$$ \left< n_i n_j n_k n_l \right> = \frac{1}{S_D} \int \hat n \hat n \hat n \hat n dA = \frac{1}{D^2+2D} (\delta_{ij}\,\delta_{kl} + \delta_{ik}\,\delta_{jl} + \delta_{il}\,\delta_{jk} )$$
