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.