Our students were tasked, using the divergence theorem $ \int_V \nabla \cdot \vec{F} \mathrm{d}^3 r = \int_{\partial V} \vec{F} \cdot \mathrm{d}\vec{a} $, to show that for any bounded, connected subset V of $\mathbb{R}^3$ with surface $\partial V$ the following identity holds:
$$\int_V \vec{r} \mathrm{d}^3r = \frac{1}{4} \int_{\partial V} \vec{r} \left( \vec{r} \cdot \mathrm{d}\vec{a} \right)$$
where $$ \vec{r} = \begin{pmatrix} x_1\\x_2\\x_3\\ \end{pmatrix}$$
Now, the sample solution considered the i-th component of the right hand side (using Einstein sum convention):
$$ \left( \frac{1}{4} \int_{\partial V} \vec{r} \left( \vec{r} \cdot \mathrm{d}\vec{a} \right) \right)_i = \frac{1}{4} \int_{\partial V} r_i r_j \mathrm{d}a_j = \frac{1}{4} \int_{\partial V} (r_i r_j) \mathrm{d}a_j \stackrel{Div. theorem}{=} \frac{1}{4} \int_V \partial_j (r_i r_j) \mathrm{d}^3r = \frac{1}{4} \int_V (\delta_{ij} r_j + 3r_i) \mathrm{d}^3r = \left( \int_V \vec{r} \mathrm{d}^3r \right)_i $$
On the other hand, many students used a corollary of the divergence theorem for gradient fields $ \int_V \nabla \phi \mathrm{d}^3 r = \int_{\partial V} \phi\mathrm{d}\vec{a} $ by realizing that for $\phi = \frac{1}{2} \left( x_1^2 + x_2^2 +x_3^2 \right) = \frac{1}{2} r^2 = \frac{1}{2} \left( \vec{r} \cdot \vec{r} \right) $ we find
$$ \int_V \vec{r} \mathrm{d}^3 r = \int_V \nabla \phi \mathrm{d}^3 r = \int_{\partial V} \phi\mathrm{d}\vec{a} = \frac{1}{2} \int_{\partial V} \left( \vec{r} \cdot \vec{r} \right) \mathrm{d}\vec{a} $$
Now, if both identities are correct, then $$\frac{1}{4} \int_{\partial V} \vec{r} \left( \vec{r} \cdot \mathrm{d}\vec{a} \right) = \frac{1}{2} \int_{\partial V} \left( \vec{r} \cdot \vec{r} \right) \mathrm{d}\vec{a}$$
which looks temptingly similar. However, nobody was able to prove this relation and it is (to me) far from trivial to see why this relation should hold.
Can anybody either show where a mistake is hidden or prove that the last identity is correct?
I realize that the sample solution effectively applies the divergence theorem to the outer product of $\vec{r}$ with itself: $$\frac{1}{4} \int_{\partial V} \vec{r} \left( \vec{r} \cdot \mathrm{d}\vec{a} \right) = \frac{1}{4} \int_{\partial V} \vec{r} \otimes \vec{r} \cdot \mathrm{d}\vec{a} = \frac{1}{4} \int_{\partial V} \vec{r} \, \vec{r}^T \cdot \mathrm{d}\vec{a} \stackrel{Div. theorem}{=} \frac{1}{4} \int_V \nabla \cdot \left( \vec{r} \, \vec{r}^T \right) \mathrm{d}^3r $$
which I find very curious. I did not know you could apply the divergence theorem to matrix-like objects, however, for every example we constructed we always found that the identitiy of the problem was correct.
Notice for any constant vector $\vec{k}$, we have
$$\begin{align} \nabla \times ( r^2 (\vec{k} \times \vec{r}) ) &= \nabla r^2 \times (\vec{k} \times \vec{r} ) + r^2 \nabla \times( \vec{k} \times \vec{r} )\\ &= 2 \vec{r} \times (\vec{k} \times \vec{r} ) + r^2 ((\nabla \cdot \vec{r}) \vec{k} - \color{red}{(\vec{k} \cdot \nabla )\vec{r}})\\ &= 2 (r^2 \vec{k} - (\vec{k}\cdot \vec{r})\vec{r}) + (3-\color{red}{1})r^2 \vec{k}\\ &= 4 r^2 \vec{k} - 2(\vec{k}\cdot\vec{r})\vec{r} \end{align} $$ Divide by $8$, integrate over any surface $S$ and apply Stoke's theorem, we obtain
$$\begin{align} \vec{k} \cdot \left[ \frac12 \int_S r^2 d\vec{a} - \frac14 \int_S \vec{r} (\vec{r}\cdot d\vec{a})\right] = & \frac18 \int_S( 4r^2 \vec{k} - 2(\vec{k}\cdot\vec{r})\vec{r}) \cdot d\vec{a}\\ = & \frac18 \int_S (\nabla \times ( r^2(\vec{k} \times \vec{r}))\cdot d\vec{a}\\ = & \frac18 \int_{\partial S} r^2(\vec{k} \times \vec{r})\cdot d\vec{r} \end{align} $$ When $S = \partial V$ is the boundary of some region $V$, the boundary of $S$ will be empty. This means the leftmost integral over $\partial S = \emptyset$ vanishes. As a result, $$\vec{k} \cdot \left[ \frac12 \int_{\partial V} r^2 d\vec{a} - \frac14 \int_{\partial V} \vec{r} (\vec{r}\cdot d\vec{a})\right] = 0$$ Since this is true for all constant vector $\vec{k}$, what's inside the square bracket vanishes and
$$\frac12 \int_{\partial V} r^2 d\vec{a} = \frac14 \int_{\partial V} \vec{r} (\vec{r}\cdot d\vec{a})$$