Derivation of Squared Angular Momentum in Spherical Coordinates While reading my textbook, I found the following:

I tried to prove the above equation by doing the following.
Knowing that : $$(\vec{A}\times\vec{B}).(\vec{C}\times\vec{D})=(\vec{A}.\vec{C})(\vec{B}.\vec{D})-(\vec{A}.\vec{D})(\vec{B}.\vec{C})$$
And making the paropiate substitutions I'm left with this:
$$\nabla^2-(\hat{r}.\vec{\nabla})(\vec{\nabla}.\hat{r})$$
I am working in spherical coordinates so we have:
$$\hat{r}.\vec{\nabla}=\frac{\partial}{\partial r}$$
$$\vec{\nabla}.\hat{r}=?$$
where in place of the question mark I found $\frac{2}{r}$ but I know it must be wrong. I used the following partial derivatives of unit vector properties in spherical coordiantes:

Also Knowing the gradient for spherical coordinates:

I know I must be wrong because the quatruple product stands for vectors, but it might not be the case for operators. I tried taking in consideration the order of the products, as a last resource i tried to take into account the commutator of the unit vector and the gradient but failed to find an expression because i'm still green in this area. I would greatly apreciate if someone could clear my doubts, I also know I could've taken the cross products directly but is there a way of doing this by taking this path? It's always enlightning to be able to do computations in different ways.  
 A: You must remember that $\textbf{r}$ is an operator and to compute $\nabla\cdot\hat r$ you must act it on a function of coordinates. Here is how I derived it.\begin{equation}
\textbf{L}^2=(\textbf{r}\times\textbf{p})\cdot(\textbf{r}\times\textbf{p})
\end{equation} Using the formula $\textbf{A}\cdot(\textbf{B}\times\textbf{C})=\textbf{C}\cdot(\textbf{A}\times\textbf{B})$ twice, we get, \begin{equation}
\textbf{L}^2=\textbf{r}\cdot(\textbf{p}\times(\textbf{r}\times\textbf{p}))
\end{equation} Using the formula for vector triple product we get,
\begin{equation}
\textbf{L}^2=\textbf{r}\cdot(p^2\textbf{r}-\textbf{p}(\textbf{p}\cdot\textbf{r}))
\end{equation} Using $[\textbf{r},p^2]=2i\hbar\textbf{p}$ and $\textbf{r}\cdot\textbf{p}-\textbf{p}\cdot\textbf{r}=3i\hbar 1$, (to prove these commutation relations just use $[x_i,p_j]=i\hbar\delta_{ij}$) we get
\begin{equation}
\textbf{L}^2=r^2p^2+i\hbar\textbf{r}\cdot\textbf{p}-(\textbf{r}\cdot\textbf{p})(\textbf{r}\cdot\textbf{p})
\end{equation} Now use the fact that $\hat{r}\cdot\nabla f(\textbf{r})=\frac{\partial f}{\partial r}$ to obtain\begin{equation}
\textbf{L}^2=r^2p^2+\hbar^2\frac{\partial }{\partial r}(r^2\frac{\partial}{\partial r})=-\hbar^2r^2[\nabla^2-\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\frac{\partial}{\partial r})]
\end{equation}
A: I worked on your line of reasoning and was able to obtain the answer. I use $\hat{r}$ to mean an operator (not unit vector), which to compute, you have to act on a function (from right to left as operators always do), as SRS have said.
After applying the quadruple  product for vectors, you obtained:
$$
\hat{L}^2 f = -\hbar^2 r^2 \nabla^2 f + \hbar^2\left(r\frac{\partial}{\partial r}\right)\left(\nabla\cdot(\hat{r}f)\right)
$$
Using the identity
$$
\frac{1}{r} \frac{\partial^2}{\partial r^2}(rf) \equiv \frac{1}{r^2} \frac{\partial}{\partial r} \left(r^2 \frac{\partial f}{\partial r} \right)
$$
I was able to tame the second term on the RHS and obtain the desired form.
Disclaimer: I am unsure if my resolution of the divergence is legitimate though.
