Using the divergence theorem to prove that $\frac{1}{|B_R(0)|} \int_{B_R(0)} M \textbf{y} . \textbf{y} dy = \frac{R^2}{ n + 2} \text{trace}(M)$ Here is the question I am trying to solve letter $(b)$ of it:



Here is a solution to it:









Which is a very very long solution. Does anyone have a more elegant and succinct solution please?
 A: 
Some few notations/properties I use here:

*

*$\partial B_1(0)$ denote the unit $n-$sphere of $\mathbb R^n$.


*For $R>0$, I denote $B_R(0):=\{\sigma _r:=r\sigma \mid r\in (0,R), \sigma \in \partial B_1(0)\}$.


*For $r>0$, I denote $\partial B_r(0):=\{\sigma _r:=r\sigma \mid \sigma \in \partial B_1(0)\}$


*If $\sigma \in \partial B_1(0)$, in particular, $\sigma $ is a unit normal vector of $\partial B_1(0)$. An infinitesimal area of $\partial B_1(0)$ is denoted by $\mathrm d \sigma$.


*Similarly, if $r>0$ is fixed, then $\mathrm d \sigma _r$ denote an infinitesimal area of $\partial B_r(0)$.


*Finally, one can prove that if $r>0$ is fixed, then $$\mathrm d \sigma _r= \mathrm d (r\sigma )=r^{n-1}\,\mathrm d \sigma .$$

Therefore
\begin{align*}
\frac{1}{|B_R(0)|}\int_{B_R(0)}My\cdot y\,\mathrm d y
&\underset{\text{Fubini}}{=}\frac{1}{|B_R(0)|}\int_0^R\int_{\partial B_r(0)}M(r\sigma )\cdot (r\sigma )\,\mathrm d \sigma _r\,\mathrm d r\\
&\underset{\mathrm d \sigma _r=r^{n-1}\mathrm d \sigma }{=}\frac{1}{|B_R(0)|}\int_0^R\int_{\partial B_1(0)}r^{n+1}M\sigma \cdot \sigma  \,\mathrm d \sigma \,\mathrm d r\\
&=\frac{1}{|B_R(0)|}\left(\int_0^Rr^{n+1}\,\mathrm d r\right)\left(\int_{\partial B_1(0)}M\sigma \cdot \sigma \,\mathrm d \sigma\right)\\
&=\frac{1}{|B_R(0)|}\cdot \frac{R^{n+2}}{n+2}\int_{\partial B_1(0)}M\sigma \cdot \sigma  \,\mathrm d \sigma\\
&\underset{\text{div. Thm.}}{=}\frac{1}{|B_R(0)|}\cdot \frac{R^{n+2}}{n+2}\int_{B_1(0)}\text{Tr}(M)\,\mathrm d x\\
&=\frac{1}{|B_R(0)|}\cdot \frac{R^{n+2}}{n+2}\cdot \text{Tr}(M)\cdot |B_1(0)|\\
&\underset{|B_R(0)|=R^n|B_1(0)|}{=}\frac{R^2}{n+2}\cdot \text{Tr}(M),
\end{align*}
as wished.
A: $
\def\b{\beta}\def\J{{\cal J}}
\def\LR#1{\left(#1\right)}
\def\vecc#1{\operatorname{vec}\LR{#1}}
\def\diag#1{\operatorname{diag}\LR{#1}}
\def\Diag#1{\operatorname{Diag}\LR{#1}}
\def\trace#1{\operatorname{Tr}\LR{#1}}
$Write the problem as
$$\eqalign{
\J = \int y^TMy\:dy = \int M:yy^T\:dy = M:\left[\int yy^T dy\right]   \\
}$$
where the constant $M$ matrix has been pulled out of the integral and the matrix inner product $(:)$ has been introduced
$$\eqalign{
A:B &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij} \;=\; \trace{A^TB} \\
A:A &= \|A\|^2_F \\
}$$
Applying the methods of this post, the remaining integral must yield an isotropic $2^{nd}$ order tensor, i.e. it must be a multiple of the identity matrix
$$\eqalign{
\J = M:\Big[\b I\Big] = \b\,\trace{M} \\
}$$
The normalization factor $\b$ is a function of the surface area of the $n$-sphere as shown in the linked post
$$\eqalign{
\b = \frac{S_n}{n^2+2n} \\
}$$
