An integral about unit vector over the whole space [$\int n_in_jn_ln_m\mathrm{d}\Omega$ where $\textbf{n}$ is the unit vector]

Question

I stuck at this integral while studying electrodynamics $$\int n_in_jn_ln_m\mathrm{d}\Omega$$ where $\textbf{n}$ is the unit vector with $$n_i=\frac{x_i}{r},\quad r=\sqrt{x^2+y^2+z^2},\quad x_1=x,x_2=y,x_3=z$$ i,j,l,m are indices 1,2 or 3,and with the Kronecker $\delta$ symbol defined as $$\delta_{ij}=\begin{cases} 1, &i=j\\ 0, &i\neq j \end{cases}$$ $\Omega$ is the solid angle and the integral is over the whole space. The answer is $$\frac{4\pi}{15}(\delta_{ij}\delta_{lm}+\delta_{il}\delta_{jm}+\delta_{im}\delta_{jl})$$ or for example, $$\int n_1n_1n_2n_2\mathrm{d}\Omega =\int_\theta\int_\varphi\frac{x^2y^2}{r^4}\sin\theta\mathrm{d}\theta\mathrm{d}\varphi=\frac{4\pi}{15}$$ (using spherical coordinates)
But how can I get this answer (without too much discussion)?
If needed, this integral originates from this problem

Answer 1

(Too long for a comment)

The choice of the system of coordinate is not prescribed regarding rotation - any new system obtained by means of rotation of the old one is also acceptable. On the other rhand, the result of the integration cannot depend on the choice of the system and should be invariant vs. any rotation. Therefore, the answer should be expressed in terms of some tensor, invariant over rotation. There are only two such tensors: $\epsilon_{ijk}$ (antisymmetric over permutation of any pair of its indices), and $\delta_{ik}$ (symmetric over permutation). $n_in_jn_kn_l$ is evidently symmetric over permutation, so we should construct the answer of the combination of $\delta_{ik}$.

The general form is $\,a(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$ - evidently, symmetric over permutation of any pair of its indices. We just have to define the constant $a$. $$\int n_in_jn_kn_l \,d\Omega=a(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$$ Making the convolution of indices, for example, by means of $\delta_{ij}\delta_{kl}$ $$\int n_in_jn_kn_l \,\delta_{ij}\delta_{kl}\,d\Omega=a(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\delta_{ij}\delta_{kl}$$ (summation over repeated indices is implied). Given that $\,\delta_{ik}\delta_{kl}=\delta_{il}$; $\,\delta_{ik}\delta_{ik}=\delta_{ii}=1+1+1=3\,$, and $n_in_j\delta_{ij}=n_1^2+n_2^2+n_3^2=(\vec n)^2=1$ $$\int 1\cdot d\Omega=4\pi=a(3\cdot3+3+3)=15\,a$$ $$a=\frac{4\pi}{15}$$