 Quadratically weighted inverse-radius Hessian integral

Away from the origin the Hessian matrix of $r^{-1}$ is $(3x_k{x}_l-r^2{\delta }_{kl})/r^5$. Multiplication by $x_i{x}_j$ makes the singularity locally integrable. Rotational invariance and the symmetries in each index pair reduce the integral to $A{\delta }_{ij}{\delta }_{kl}+B({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})$. The Laplacian trace gives $3A+2B=0$, while contracting the two pairs as $k=i,l=j$ gives $3A+12B={\int }_{r<a}2/rdV=4\pi a^2$. Solving these two equations yields the displayed expression.