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!
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...)