I am computing $\nabla \nabla$ of a function $\mu_{\gamma}=-\frac{e^{-\gamma}}{4 \pi r}$, and get the following (since the function only depends on $r$).
$$ \begin{split} \nabla \mu_{\gamma} &= \frac{\textbf{r}}{|\textbf{r}|} \frac{\text{d}\mu_{\gamma}}{\text{d}|\textbf{r}|}.\\ \\ \nabla (\nabla \mu_{\gamma}) &= \frac{\nabla (\textbf{r})}{|\textbf{r}|}\frac{\partial}{\partial r}\mu_{\gamma} + \hat{r} \nabla \bigg( \frac{\partial}{\partial r} \mu_{\gamma} \bigg)\\ &=\frac{\textbf{r}}{|\textbf{r}|} \frac{\text{d} \mu_{\gamma}}{\text{d} |\textbf{r}|} + \frac{\textbf{r}}{|\textbf{r}|}\frac{\textbf{r}}{|\textbf{r}|} \frac{\partial }{\partial r} \bigg( \frac{\partial}{\partial r} \mu_{\gamma} \bigg) \\ &= \frac{\delta_{ij}}{|\textbf{r}|}\frac{\text{d} \mu_{\gamma}}{\text{d}|\textbf{r}|} + \frac{\textbf{r} \otimes \textbf{r}}{|\textbf{r}|^2} \frac{\text{d}^2 \mu_{\gamma}}{\text{d} |\textbf{r}|^2}. \end{split}$$ However, once I compute the first and second derivatives of the function $\mu_{\gamma}$, what I get does not match the result with Mathematica, in the sense that when I subtract the result of my computation from the result which Mathematica gives, the result is not $0$. Is there something that I am doing wrong?
Edit: I realised that I had a minus sign in the wrong place and for that reason the derivatives were coming out slightly wrong and the expressions were not matching the Mathematica result.
$\newcommand{\d}{\operatorname{d}\!}$ It is important to realise that the gradient is an operator that acts only on scalar functions, and produces a tensor of degree 1. After that, you need to use the covariant derivative. That's probably where you went wrong, since the covariant derivative of the basis is not always zero: that is Christoffel symbols are for. If $\mu$ is a scalar function, $\d{x}^i$ are the basis 1-forms, and $\Gamma^{i}{}_{jk}$ are the Christoffel symbols, the equations $$\nabla \mu = \sum_i\frac{\partial \mu}{\partial x^i}\d x^{i} \tag{1}$$ $$\nabla (\d{x}^{i}) = -\sum_{j,k}\Gamma^{i}{}_{jk}\d{x}^j\otimes\d{x}^{k} \tag{2}$$ are valid in any coordinate system. In the case polar coordinates for the plane, the Christoffel symbols are $\Gamma^{2}{}_{21}=\Gamma^{2}{}_{12}=\frac{1}{r}$, $\Gamma^{1}{}_{22}=-r$ and every other entry is zero. Hence, for the case that $\mu$ depends only on $r$, by virtue of the Leibniz rule and equation $(1)$, we have $$\nabla\mu = \frac{\partial \mu}{\partial r}\d{r}$$ $$\nabla\nabla\mu = \nabla\left(\frac{\partial \mu}{\partial r}\right)\d{r} + \frac{\partial \mu}{\partial r}\nabla(\d{r})$$ now from $(2)$, we have $\nabla(\d{r}) \equiv \nabla(\d{x}^{1}) = -\sum_{i,j}\Gamma^{1}{}_{ij}\d{x}^i\otimes\d{x}^{j}=-\Gamma^{1}{}_{22}\d{x}^2\otimes\d{x}^2\equiv r\d{\theta}\otimes\d{\theta}$; hence we get $$\nabla\nabla\mu = \frac{\partial^{2} \mu}{\partial r^{2}}\d{r}\otimes\d{r} + r\frac{\partial \mu}{\partial r}\d{\theta}\otimes\d{\theta}$$ As you see, since the Hessian is a bilinear form, it is an element of the second tensor power of the cotangent space at each point. You could turn the Hessian into a contravariant object using the musical isomorphism $\sharp$ provided by the metric: send $\d{x}^{i}\mapsto g^{ij}v_{j}$ (where $v_i$ are the holonomic basis vectors). In polar coordinates this gives $\d{r}\equiv\d{x}^{1}\mapsto g^{11}v_{1} \equiv v_{r}$ and $\d{\theta}\equiv\d{x}^{2}\mapsto g^{22}v_{2} \equiv \frac{1}{r^{2}}v_{\theta}$. Hence you get $$(\nabla\nabla\mu)^{\sharp} = \frac{\partial^{2} \mu}{\partial r^{2}}v_{r}\otimes v_{r} + \frac{1}{r^{3}}\frac{\partial \mu}{\partial r}v_{\theta}\otimes v_{\theta}$$
Moreover, people often like to normalize their basis. We have $|v_{r}|=1$ and $|v_\theta|=r$, so the normal basis will be $e_r = v_r$ and $e_\theta = \frac{1}{r}v_\theta$. Hence you get
$$(\nabla\nabla\mu)^{\sharp} = \frac{\partial^{2} \mu}{\partial r^{2}}e_{r}\otimes e_{r} + \frac{1}{r}\frac{\partial \mu}{\partial r}e_{\theta}\otimes e_{\theta}$$