I am quite confused on how to evaluate successive mixed gradients of functions of the form $f(r)/r$, appearing often as potential in physics. In particular, I was reading this paper, where Eq. (3.3) is
$$ H_{ijk}(uR) = (-\nabla^2\delta_{ij} + \nabla_i \nabla_j)\nabla_k \frac{e^{-uR}}{R} $$
While the paper gives the result of this expression, I tried to evaluate it myself, but I can't see how the components work for the successive gradients.
In particular, I understand that $\nabla \frac{1}{R}=-\frac{\hat{R}}{R^2}$, so I'm pretty sure that $\nabla_k \frac{1}{R}=-\frac{\hat{R}_k}{R^2}$. Then, in the next step, I should take $\nabla_j$ of this result, and since I assume $\nabla_j \hat{R}_k=0$, I would get
$\nabla_j \nabla_k \frac{1}{R}=\delta_{jk} \frac{2 \hat{R}_k}{R^3}$
I am pretty sure that already at this step my derivation is incorrect, and I assume it is because the mixed derivative does not give $\delta_{jk}$, but I can't see why.
Could anyone walk me through how to obtain the components of this successive gradient expression?
Denote $$f(r)=\frac{e^{-ur}}{r}$$ with $r = \| \mathbf{x} \|$. Using the concept of differentials, it holds $$ df= - \left[ u+\frac{1}{r} \right]f dr= - \left[ u+\frac{1}{r} \right]\frac{f}{r} \mathbf{x}^T d\mathbf{x} $$ since $dr = \frac{1}{r} \mathbf{x}^T d\mathbf{x}$.
From here, denote $$ \phi = \frac{\partial f}{\partial x_k} = - \left[ u+\frac{1}{r} \right] \frac{f}{r} x_k $$
The gradient writes $$ \mathbf{g} = \frac{\partial \phi}{\partial \mathbf{x}} = \left[3+3ur+u^2 r^2\right] \frac{f x_k}{r^4} \mathbf{x} - \left[ u+\frac{1}{r} \right] \frac{f}{r} \mathbf{e}_k $$ with $\mathbf{e}_k$ the canonical vector basis for $\mathbb{R}^3$.
The Hessian follows similarly and yields the requested formula.