A continuación, calcularemos las relaciones de conmutación entre los distintos operadores del momento angular, para determinar cuales de ellas se pueden conocer simultáneamente. \begin{equation} [\hat{l}_x , \hat{l}_y]=\hat{l}_x \hat{l}_y - \hat{l}_y \hat{l}_x \end{equation} Para que no resulten expresiones demasiado largas calculemos los dos términos por separado \begin{equation} \hat{l}_x\hat{l}_y f=\left[-i\hbar\left(y\frac{\partial}{\partial z}-z\frac{\partial}{\partial y}\right)\right]\left[-i\hbar\left(z\frac{\partial f}{\partial x}-x\frac{\partial f}{\partial z}\right)\right]= \end{equation} \begin{equation} =-\hbar^2\left(y\frac{\partial f}{\partial x}+yz\frac{\partial^2 f}{\partial z \partial x}-yz\frac{\partial^2 f}{\partial z^2}-z^2\frac{\partial^2 f}{\partial y \partial x}+zx\frac{\partial^2 f}{\partial y \partial z}\right) \end{equation}
Ahora calculamos el segundo termino del conmutador \begin{equation} \hat{l}_y\hat{l}_x f=\left[-i\hbar\left(z\frac{\partial}{\partial x}-x\frac{\partial}{\partial z}\right)\right]\left[-i\hbar\left(y\frac{\partial f}{\partial z}-z\frac{\partial f}{\partial y}\right)\right]= \end{equation} \begin{equation} =-\hbar^2 \left( zy\frac{\partial^2 f}{\partial x \partial z}-z^2\frac{\partial^2 f}{\partial x \partial y}-zy\frac{\partial^2 f}{\partial z^2}+x\frac{\partial f}{\partial y}+xz\frac{\partial^2 f}{\partial z \partial y}\right) \end{equation} Sustituyendo en el conmutador: \begin{equation} [\hat{l}_x , \hat{l}_y]=\hat{l}_x \hat{l}_y - \hat{l}_y \hat{l}_x =-\hbar^2\left(y\frac{\partial}{\partial x}-x\frac{\partial}{\partial y}\right)=-\frac{\hbar}{i}\left(-i\hbar\right)\left(x\frac{\partial}{\partial y}-y\frac{\partial}{\partial x}\right)=i\hbar\hat{l}_z \end{equation} De forma análoga se calculan los comutadores entre las restantes componentes del momento angular. Resumiendo: \begin{eqnarray} \left[\hat{l}_x , \hat{l}_y\right] &=& i\hbar\hat{l}_z\\ \left[\hat{l}_y , \hat{l}_z\right] &=& i\hbar\hat{l}_x\\ \left[\hat{l}_z , \hat{l}_x\right] &=& i\hbar\hat{l}_y \end{eqnarray} Vamos a evaluar los conmutadores de $\hat{l}^2$
\begin{equation} \left[\hat{l}^2, \hat{l}_x \right]=\left[\hat{l}^2_x + \hat{l}^2_y + \hat{l}^2_z ,\hat{l}_x\right]=\underbrace{\left[\hat{l}^2_x , \hat{l}_x\right]}_{0}+\left[\hat{l}^2_y , \hat{l}_x\right]+\left[\hat{l}^2_z , \hat{l}_x\right]=\left[\hat{l}_y\hat{l}_y , \hat{l}_x\right]+\left[\hat{l}_z\hat{l}_z , \hat{l}_x\right]= \end{equation}
Aplicando la propiedad $\left[\hat{A}\hat{B},\hat{C}\right]=\left[\hat{A},\hat{C}\right]\hat{B} + \hat{A}\left[\hat{B},\hat{C}\right]$
\begin{equation} =\left[\hat{l}_y,\hat{l}_x \right]\hat{l}_y + \hat{l}_y\left[\hat{l}_y,\hat{l}_x \right]\hat{l}_y + \left[\hat{l}_z,\hat{l}_x \right]\hat{l}_z + \hat{l}_z\left[\hat{l}_z,\hat{l}_x \right]\hat{l}_z =\\-i\hbar\hat{l}_z\hat{l}_y + \hat{l}_y (-i\hbar\hat{l}_z) +i\hbar\hat{l}_y \hat{l}_z + \hat{l}_z(i\hbar\hat{l}_y)=0 \end{equation}