Square of angular momentum operator in spherical coordinates I want to compute the square of the angular momentum operator in spherical coordinates. I already know how the cartesian components look like: 
\begin{align}
L_x &= -i\hbar \left(-\sin\phi\,\partial_{\theta} - \cos\phi\,\cot\theta \,\partial_{\phi}\right)\\
L_y &= -i\hbar \left(\cos\phi\,\partial_{\theta} - \sin\phi\,\cot\theta \,\partial_{\phi}\right)\\
L_z &= -i\hbar\,\partial_{\phi}
\end{align}
The square should then be given by: $L^2 = L_x^2 + L_y^2 + L_z^2$
With that, what I would get is this: 
$$L^2 = -\hbar^2 \left(\partial_{\theta}^2 + \cot^2\theta\,\partial_{\phi}^2 +\partial_{\phi}^2 \right) = -\hbar^2 \left(\partial_{\theta}^2 + (\cot^2\theta + 1)\partial_{\phi}^2\right) = -\hbar^2\left(\partial_{\theta}^2 + \frac{1}{\sin^2\theta}\partial_{\phi}^2\right)$$
From the 3 terms in the bracket, the first term comes from $\cos^2\phi + \sin^2\phi = 1$, the second term from the same rule, and the 3rd term is the square of the $L_z$ component. The mixed terms you get, when you square $L_x$ and $L_y$ should cancel: 
$$2\,\sin\phi\,\partial_{\theta}\,\cos\phi\,\cot\theta\,\partial_{\phi} - 2\,\cos\phi\,\partial_{\theta}\,\sin\phi\,\cot\theta\,\partial_{\phi} = 0 $$
However, the solution I find everywhere is: 
$$L^2 = \left(\partial_{\theta}^2 + \cot\theta\,\partial_{\theta} + \frac{1}{\sin^2\theta}\partial_{\phi}^2\right)$$
I can't find where the error is.
 A: The error is that
\begin{align}
&L_x^2f(\theta,\phi)\\
&=\left( -i\hbar \left(-\sin\phi\,\partial_{\theta} - \cos\phi\,\cot\theta \,\partial_{\phi}\right)\right)
\left( -i\hbar \left(-\sin\phi\,\partial_{\theta} - \cos\phi\,\cot\theta \,\partial_{\phi}\right)\right)f(\theta,\phi)\\
&\ne-\hbar^2\left(\sin^2\phi\partial^2_{\theta^2}+2\cos\phi\sin\phi\cos\theta 
\partial^{2}_{\phi,\theta}
+\cos^2\phi\cot^2\theta\partial^2_{\phi^2}\right)
f(\theta,\phi) 
\end{align}
since the derivatives $\partial_\theta$ and $\partial_\phi$ do not commute with the functions $\sin\phi$ and $\cos\phi\cot\theta$ in the expressions of the operators.  In other words, there are additional cross-terms in $L_x^2$ beyond $2\cos\phi\sin\phi\cos\theta 
\partial^{2}_{\phi,\theta}\partial_{\phi\theta}$; these extra cross-terms are linear in $\partial_\theta$ and $\partial_\phi$ and come - for instance - from the product rule applied to
\begin{align}
-\sin\phi\partial_\theta
\left(-\cos\phi\cot\theta\partial_\phi f(\theta,\phi)\right)&=
-\sin\phi\cos\phi\csc^2\theta\partial_\phi f(\theta,\phi)\\
&\qquad +\sin\phi\cot\theta\cos\phi\partial^2_{\theta\phi}f(\theta,\phi)\, .
\end{align}
The same argument applies to $L_y^2$.  $L_z^2$ is easy since the coefficients in front of the derivative are constant.
