How to rigorously show tensor identities using symmetry arguments? I am wondering how to rigorously show tensor identities such as the following.

Let $n$ denote the radial unit vector in $D$ dimensions. Then $\langle n_i n_j \rangle = \frac 1 D \delta_{ij}$ and $\langle n_i n_j n_k n_l\rangle = \frac{1}{D^2+2D} \left[\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} \right]$, where the brackets indicate averaging over all orientations.

A 'proof' I sometimes encounter is to guess that (e.g.) $\langle n_i n_j \rangle = A \delta_{ij}$ by symmetry, and then find the value of $A$ by taking the trace of both sides.
I am wondering how to make the symmetry argument rigorous. I have a vague hunch that this problem is related to representations of the rotation group, but I am not sure. I would also be very happy with references that deal with this type of problems.
 A: 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} )$$ 
