Averaged dot product In one of the textbooks I came across this:
"the average over the angles of the scalar product such as $(\vec{n} \cdot \vec{a}) (\vec{n} \cdot \vec{b})$ is equal to $\frac{1}{3}(\vec{a} \cdot \vec{b})$"
The author also gives the average of $(\vec{n} \cdot \vec{a}) (\vec{n} \cdot \vec{b})(\vec{n} \cdot \vec{c}) (\vec{n} \cdot \vec{d})$ which is equal to $\frac{1}{15}[(\vec{a} \cdot \vec{b})(\vec{c} \cdot \vec{d})+(\vec{a} \cdot \vec{d})(\vec{c} \cdot \vec{b})+ (\vec{a} \cdot \vec{c})(\vec{b} \cdot \vec{d})]$.
In all of these $\vec{n} = \vec{r}/r$ and $\vec{a},\vec{b},\vec{c},\vec{d}$ are constant vectors. How to prove it? What does "average over angles" mean, just $d\theta$ and $d\phi$ or the full spherical angle $d\phi \;d\theta \sin \theta$?
Thanks in advance!
 A: First formula :
Here are two different proofs :
Method 1 :
Let $u \cdot v$ denote the dot product of $u$ and $v$.
Using expansion :
$$(n \cdot (a-b))^2 = (n \cdot a)^2 + (n \cdot b)^2 - 2 (n \cdot a)(n \cdot b),$$
we can express :
$$(n \cdot a)(n \cdot b) = \frac12 \left((n \cdot a)^2 + (n \cdot b)^2 - (n \cdot (a-b))^2\right)$$
which gives the equivalent identity in terms of average values :
$$\overline{(n \cdot a)(n \cdot b)} = \frac12 \left(\overline{(n \cdot a)^2 }+ \overline{(n \cdot b)^2} - \overline{(n \cdot (a-b))^2}\right)\tag{*}$$
Thus, we are brought back to compute the average value $\overline{(c \cdot n)^2} $of $(c \cdot n)^2$ for a fixed vector $c$ and a variable vector $n \in \mathbb{S}$, the unit sphere.
We are going to prove that
$$\overline{(n \cdot c)^2} = \frac13\|c\|^2 \tag{**}$$
Indeed, taking into account (**) in (*), we will get :
$$\overline{(n \cdot a)(n \cdot b)} = \frac12 \frac13 \left(\|a\|^2 + \|b\|^2 - \|a-b\|^2\right)=$$
$$= \frac12 \frac13 \left(\|a\|^2 + \|b\|^2 - \|a\|^2 - \|b\|^2+2 a \cdot b\right)=\frac13 a \cdot b$$
as desired.
It remains to prove (**). We can consider WLOG $c$ to be with coordinates $(0,0,k)$ with $k:=\|c\|$. Let $(H)$ be the hemisphere of $\mathbb{S}$ with $z>0$. If $n=(x,y,z) \in (H)$, 
$$(n \cdot c)^2=k^2z^2.\tag{***}$$
(H) can be "sliced" into spherical segments of height d$z$ with area $2\pi dz$ (http://mathworld.wolfram.com/SphericalSegment.html
). Integrating all these contributions on hemisphere $(H)$, we obtain, using (***) ; 
$$\overline{(n \cdot c)^2}=\color{red}{\frac{1}{2 \pi}}\int_0^1 (k^2 z^2) (2 \pi dz) =  \frac13 \|c\|^2 $$
i.e., result (**). Why this front division by $\color{red}{2 \pi}$ ? Because the sum of all the spherical segments areas is equal to the value  $2 \pi$ which is the area of hemisphere $(H)$. Said otherwise,  dividing $2 \pi dz $ by $2 \pi$ provides a pdf (probability distribution function) whose total "mass" is 1. 

Method 2 :
Identifying $a,b$ and $n$ (all assumed unit norm vectors) with the associated column vectors of coordinates with respect to a certain basis, let us consider the following symmetric matrix :
$$M:=ab^T+ba^T$$
Let the dot product of $a$ and $b$ be denoted :
$$\delta:=a \cdot b \ \ \text{identified with} \ \ a^Tb.$$
Quadratic form :
$$Q(n):=n^TMn$$
is such that $Q(n)=(n^Ta)(b^Tn)+(n^Tb)(a^Tn)=2(n^Ta)(n^Tb),$
i.e., twice the expression for which we have to determine its mean.
With the hypothesis $\|n\|=1$, $Q$ is the Rayleigh quotient (https://en.wikipedia.org/wiki/Rayleigh_quotient) associated with matrix $M$. The set of values of a Rayleigh quotient is known to be expressed as the set of all weighted averages (with positive weights) of its eigenvalues.   
The 3 eigenvalues of $M$ being
$$\delta-1 \leq 0 \leq \delta +1\tag{3}$$
(see appendix below), we have to compute the mean of all expressions :
$$w_1(\delta-1)+w_20+w_3(\delta+1) \ \text{such that} \ w_1+w_2+w_3=1$$
which is :
$$\frac13(\delta-1)+\frac13(\delta+1)=\frac23\delta $$
($\delta-1$ and $\delta+1$ playing symmetrical rôles).
Dividing by 2 yields the result.
Appendix Explanation of (3) ;
The eigenvalues in (3) are associated with  eigenvectors $a-b$, $a \times b$, $a+b$ resp.
For example, for the first vector $a-b$ :
$$M(a-b)=(ab^T+ba^T)(a-b)=$$
$$=a \delta -a\underbrace{(b^Tb)}_{=1}+b\underbrace{(a^Ta)}_{=1}-b \delta=(\delta-1)(a-b)$$
establishing that $a-b$ is an eigenvector associated with eigenvalue $\delta-1$.
2nd formula :
I have no proof for it but @ThomasTuna has given a valuable one using a tensorial method ; another very interesting answer, using almost the same tools : 
Surface integral of normal components summations on a sphere
Another reference whose interest is in the question mentionning equivalent formulas in $N$ dimensions with $3$ replaced by $D$, and $15$ replaced by $(D^2+2D)$: 
How to rigorously show tensor identities using symmetry arguments?
(without answer, but a valuable comment by Ted Shifrin), and finally
Eigenvalues of a rank 2 tensor defined by an integral
(without any answer...)
A: FIRST IMPRESSION: This problem was written in "Some problems of nuclear theory" by Akchiezer, this is a physics textbook.
Classically trained physicists solve problems with tensors/indices (i.e Einstein Notation) and arguments of symmetry. Using such tools both asserted relations follow in a few lines of algebra. 
GAME PLAN:I will convert the problem into tensor contraction and then I will deduce the answer by symmetry. 
NOTE: I will use explicit summation, rather than Einstein's implicit summation convention.
DERIVATION: Consider $(\vec a \cdot \vec n )(\vec b \cdot \vec n)$.
$$(\vec a \cdot \vec n )(\vec b \cdot \vec n) = \sum_{ij} a_i n_i b_j n_j$$
$$ = \sum_{ij} a_i  b_j n_i n_j$$
Let $M_{ij} \equiv a_i b_j $ and $N_{ij} \equiv n_i n_j$.
$$ = \sum_{ij} M_{ij} N_{ij} \tag{X}$$
Eq.X is how a classically trained physicist sees this problem. As the contraction of 2 different 2-tensors, here $M$ and $N$ are obviously both rank 2 tensors.
By arguments of symmetry (see Appendix), we conclude that 
$$N_{ij}=\frac{1}{3} \delta_{ij} \tag{XX}$$
Combine equations (*), (**), and the definition of $M$ to arrive at:
$$(\vec a \cdot \vec n )(\vec b \cdot \vec n) = \sum_{ij} M_{ij} N_{ij} = \sum_{ij} (a_i b_j) (\frac{1}{3} \delta_{ij})= \frac{1}{3} \sum_{ij} a_i \delta_{ij} b_j  = \frac{a^Tb}{3}$$

REVIEW AND CONNECTIONS: This is exactly, the result already provided by Jean Marie's Method 1 and 2. So this is a third method to arrive at your first relation. It has an interesting connection with Jean's method 2. In this method Jean regroups the vectors as 
$$(\vec{n} \cdot \vec{a}) (\vec{n} \cdot \vec{b}) = \vec{n}^T\vec{a} \vec{b}^T\vec{n}^T$$
and distinguishes the the outer product of vectors $\vec{a} \vec{b}^T$ as a matrix.
Then solves the problem via properties of $M=\vec{a} \vec{b}^T +\vec{b} \vec{a}^T$.
Analogously, I have grouped $$(\vec{n} \cdot \vec{a}) (\vec{n} \cdot \vec{b}) = \vec{a}^T\vec{n} \vec{n}^T\vec{b}^T$$
and distinguished the the outer product of vectors $\vec{n} \vec{n}^T=N$ as a matrix. Then solved the problem via properties of $N$.
These approaches solve the same problem because of commutativity of the inner product.

Now I will prove the second statement.
GAME PLAN: Exactly as before, I will convert the problem into a tensor contraction and then I will deduce the answer by symmetry. 
DERIVATION: Consider 
$$(\vec{n} \cdot \vec{a}) (\vec{n} \cdot \vec{b})(\vec{n} \cdot \vec{c}) (\vec{n} \cdot \vec{d}) = \sum_{ijkl} a_i n_i b_j n_j c_k n_k d_l n_l $$
Let $M_{ijkl} = a_i b_j c_k d_l$ and $N_{ijkl}=n_i n_j n_k n_l$.
$$= \sum_{ijkl} M_{ijkl}N_{ijkl} \tag{*}$$
We now reduced the right hand side to a contraction of $4$-tensors 
As before, by arguments of symmetry (see appendix)
$$N_{ijkl}=n_i n_j n_k n_l = \frac{1}{15}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} ) \tag{**}$$.
Combining Eq.**, Eq.*, and the definition of $M$, we arrive at the solution
$$\sum_{ijkl} M_{ijkl}N_{ijkl} = \sum_{ijkl} ( a_i b_j c_k d_l ) \frac{1}{15}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} )$$
$$ =\frac{1}{15}[(\vec{a} \cdot \vec{b})(\vec{c} \cdot \vec{d})+(\vec{a} \cdot \vec{d})(\vec{c} \cdot \vec{b})+ (\vec{a} \cdot \vec{c})(\vec{b} \cdot \vec{d})]$$
-----------------ARGUMENTS OF SYMMETRY------------------
Consider a rotated vector $n'$ where $n'_i=\sum_{j} R_{ij} n_j = R \vec n$ (Here $R$ is defined as a rotation matrix - https://en.wikipedia.org/wiki/Rotation_matrix). Then the rotated 2-tensor $N'$ is defined as,
$$N' = R \vec{n} (R \vec n)^T = R \vec{n} \vec{n}^T R^T =R N R^{-1}$$
Thus the rotation of $\vec n$ results in a basis transformation of $N$ (https://en.wikipedia.org/wiki/Matrix_similarity).
NOW HERE'S THE SYMMETRY, Since the random vector $\vec n$ is sampled from spherically symmetric distribution on a sphere, then $N$ must invariant w.r.t. rotations, thus
$$N'=N$$
$$N'_{ij}= R_{ik}n_k R_{jl} n_l = R_{ik} R_{jl} n_k n_l = R_{ik} R_{jl} \delta_{kj} n_j \delta_{li} n_i =  R_{ij} R_{ji} n_i n_j = \delta_{ij} n_i n_j = N_{ij}$$
Thus, we conclude the only such two 2-tensors are proportional to the identity so that $N_{ij} = c_0 \delta_{ij}$. 
Furthermore, $\vec n$ is normalized. Thus $\sum_i n_i n_i =1$ which is equivalently, $tr\{ N\} = \sum_{ij} N_{ij} \delta_{ij} = \sum_i c_0  \delta_{ij} \delta_{ij} = 3*c_0$, it follows that,
$$N_{ij}=\frac{1}{3} \delta_{ij} \tag{XX}$$
This is result Eq.XX
Now consider the second problem, Eq.**, which is a generalization of the first. As before we still have spherical symmetry, so
$$N_{ijkl}=n_i n_j n_k n_l = N'_{ijkl}$$
$$N'_{ijkl}=\sum_{wxyz} R_{iw}n_w R_{jx}n_x R_{ky}n_y R_{lz} n_z = N_{ijkl}$$
$$N'_{ijkl}=\sum_{wxyz} R_{iw}R_{jx} R_{ky} R_{lz} n_w n_x n_y n_z = N_{ijkl}$$
Now unlike above we have freedom to mix indices contractions, I could insert $n_z = \delta_{zk} n_k$ or $n_z = \delta_{zj} n_j$ or $n_z = \delta_{zi} n_i$, (from here commutativity determines all other possible combinations) and arrive at 
$$N_{ijkl}= c_1 (\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} )$$
As before the normalization condition is,
$$1 =\sum_i n_i n_i$$
$$= \sum_i n_i n_i \sum_k n_k n_k $$
$$= \sum_{ijkl} c_1 (\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} ) \delta_{ij}\delta_{kl}$$
$$=9+3+3$$
Thus, we conclude that 
$c_1 = \frac{1}{15}$
