Averaged dot product
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...)
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}$