Solution (source code)

= Solution

Choose the symmetric <tensor field>
$$
\Sigma_{kj}=(\partial_k\psi)(\partial_j\psi)-\frac12\delta_{kj}|\nabla\psi|^2.
$$
The <product rule> and equality of mixed <partial derivatives> give
$$
\begin{aligned}
\partial_k\Sigma_{kj}
&=(\nabla^2\psi)\partial_j\psi+(\partial_k\psi)\partial_k\partial_j\psi-\frac12\partial_j|\nabla\psi|^2\\
&=(\nabla^2\psi)\partial_j\psi.
\end{aligned}
$$
This is the <Euclidean scalar-gradient stress identity>. Substitute it into part (e), and use $\Sigma_{ki}\hat n_k=(\partial_n\psi)\partial_i\psi-\tfrac12|\nabla\psi|^2\hat n_i$, where $\partial_n\psi=\hat{\mathbf n}\cdot\nabla\psi$. Reassembling the components gives the required pure <surface integral>:
$$
\boxed{\int_V\mathbf r\,\nabla\cdot[(\nabla^2\psi)\nabla\psi]\,dV=\int_S\left[\mathbf r(\nabla^2\psi)\partial_n\psi-(\partial_n\psi)\nabla\psi+\frac12|\nabla\psi|^2\hat{\mathbf n}\right]dS.}
$$
Sufficient smoothness of $\psi$ is understood for the displayed <derivatives> and applications of the <divergence theorem>.