The first two terms and the last two terms cancel independently. That the last two terms cancel follows from substituting \(r = \sqrt{x^2 + y^2 + z^2} \) and using the definition of Poisson brackets as a sum of partial derivatives.
I have oftentimes found the presentation of angular momentum operators skipping not-so-trivial steps (e.g.,
In classical mechanics the angular momentum of a particle with position \(\vec{r}\) and momentum \(\vec{p}\) is defined as \(\vec{L} = \vec{r} \times \vec{p} \). In cartesian coordinates we have
\begin{align} L_x &= yp_z - zp_y \\ L_y &= zp_x - xp_z \\ L_z &= xp_y - yp_x. \\ \end{align}A problem is said to possess spherical symmetry when the potential depends only on the raidal coordinate \(r\), $$V = V(r).$$ Such potentials are called central potentials. For central potentials we can write the Hamilitonian as $$H = H(x,y,z,p_x,p_y,p_z) = \frac{\vec{p}^2}{2m} + V(r)$$ with \(r = \sqrt{x^2 + y^2 + z^2}\). Angular momentum is a useful quantity in problems with spherical symmetry because it is a conserved quantity of such systems. In classical mechanics the condition that a quantity be conserved is that its Poisson bracket with the Hamiltonian vanish: $$\frac{dF}{dt}(q,p,t) = \frac{\partial F}{\partial t} + \frac{\partial F}{\partial q}\dot{q} + \frac{\partial F}{\partial p}\dot{p}= \frac{\partial F}{\partial t} + \Big \{ F,H \Big \}.$$ We can demonstrate directly that \(\vec{L}\) is conserved using Poisson brackets. Because a vector is conserved if and only if each of its components is conserved, we will show only that \(L_x\) is conserved. The proofs that \(L_y\) and \(L_z\) are conserved follow from cyclically permuting \(x,y,\) and \(z\). Using the definitions and the fundamental Poisson bracket we can easily check that $$\frac{dL_x}{dt} = \Big \{ L_x,H \Big \} = 0$$
The first two terms and the last two terms cancel independently. That the last two terms cancel follows from substituting \(r = \sqrt{x^2 + y^2 + z^2} \) and using the definition of Poisson brackets as a sum of partial derivatives.
To make use of the spherical symmetry inherent in a problem with a central potential we should express our equations in spherical coordinates. It will turn out that doing so leads to the convenient form for the Hamiltonian $$H = H(r,\theta,\phi, p_r, p_\theta, p_\phi) = \frac{p_r^2}{2m} +\frac{L^2}{2mr^2} + V(r)$$ And, noting that \(L\) is a constant, we see that the angular dependence drops out! This makes the Hamiltonian 1-dimensional instead of 3-dimensional — a great simplification. However, showing that this formula is true requires some work. What we are really doing is a canonical transformation from cartesian coordinates and momenta to spherical coordinates and momenta: $$x,y,z,p_x,p_y,p_z \rightarrow r, \theta, \phi, p_r,p_\theta, p_\phi $$ When we solve for the new Hamiltonian we know one leg of the transformation, namely, the coordinate transformation, but we do not know how the momenta transform \begin{align} x &= r\sin\theta \cos\phi & p_x =\; ?\\ y &= r\sin\theta \sin\phi & p_y =\; ? \\ z &= r\cos\theta & p_z = \;? \end{align} Any time we are doing this kind of transformation, where we know one leg but not the other, we need to find the generating function. In this case, where we know the old coordinates in terms of the new, we use a type 3 generating function. $$F_3 = F_3(\vec{p}, \vec{Q}) \qquad q_\alpha = - \frac{\partial F_3}{\partial p_\alpha} \quad P_\alpha = -\frac{\partial F_3}{\partial Q_\alpha}$$ (In practice the best way to decide the right generating function is either experience, or looking at all four and seeing which one fits.) Here we have \begin{align} -\frac{\partial F_3}{\partial p_x} = x &= r\sin\theta \cos\phi \\ -\frac{\partial F_3}{\partial p_y} = y &= r\sin\theta \sin\phi \\ -\frac{\partial F_3}{\partial p_z} = z &= r\cos\theta \end{align} Which we can easily integrate to determine $$F_3(p_x,p_y,p_z, r, \theta, \phi) = - p_x r\sin\theta \cos\phi - p_y r \sin\theta\sin\phi - p_zr\cos\theta$$ Now we can use the other leg to write the new momenta \begin{align} p_r &= p_x\sin\theta\cos\phi + p_y \sin\theta\sin\phi + p_z\cos\theta \\ p_\theta &= p_xr\cos\theta\cos\phi + p_y r\cos\theta\sin\phi - p_zr\sin\theta \\ p_\phi &= -p_xr\sin\theta\sin\phi + p_y r\sin\theta\cos\phi. \end{align} In this result you can notice that \(p_\phi = -p_xy + p_yx = L_z \), but this is getting a little ahead of ourselves. Inverting these is not too hard but does require a little work which I have outsourced to SymPy in code pasted at the bottom of this page. The end result is the spherical variable representation of the cartesian momenta \begin{align} p_x &= p_r\sin\theta\cos\phi + \frac{\cos\theta\cos\phi}{r}p_\theta - \frac{\sin\phi}{r\sin\theta}p_\phi \\ p_y &= p_r\sin\theta\sin\phi + \frac{\cos\theta\sin\phi}{r}p_\theta + \frac{\cos\phi}{r\sin\theta}p_\phi \\ p_z &= p_r\cos\theta - \frac{\sin\theta}{r}p_\theta \end{align}
This completes the second leg of the spherical polar transformation. We can now substitute these into the Hamiltonian. Doing so is fairly messy but the result simplifies nicely $$H = \frac{p_r^2}{2m} + \frac{p_\theta^2}{2mr^2} + \frac{p_\phi^2}{2mr^2\sin^2\theta} + V(r)$$ Likewise, writing the terms of angular momentum in spherical coordinates gives \begin{align} L_x &= \left(-\sin(\phi) p_\theta - \cot(\theta) \cos(\phi) p_\phi\right) \\ L_y &= \left( \cos(\phi) p_\theta - \cot(\theta) \sin(\phi) p_\phi\right) \\ L_z &= p_\phi \end{align} from which we can check $$L^2 = p_\theta^2 + \frac{1}{\sin^2\theta}p_\phi^2$$ which allows us to to reduce the Hamiltonian to 1 dimension, as above.
from sympy import *
init_printing()
r, theta, phi = symbols("r theta phi")
# The cartesion momentum to spherical momentum transformation is
#
# Pr = px * sin(th)*cos(phi) + py * sin(th)*sin(phi) + pz * cos(th)
# Pth = px * r*cos(th)*cos(phi) + py * r*cos(th)*sin(phi) - pz * r*sin(th)
# Pphi = -px * r*sin(th)*sin(phi) + py * r*sin(th)*cos(phi)
# In matrix notation we have
#
# spherical_momentum_vector = J * cartesian_momentum_vector
J = Matrix(( [ sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)],
[ r*cos(theta)*cos(phi), r*cos(theta)*sin(phi), -r*sin(theta)],
[-r*sin(theta)*sin(phi), r*sin(theta)*cos(phi), 0] ))
# The inverse transformation is
Jinverse = J.inv().applyfunc(simplify)
# ⎡ cos(φ)⋅cos(θ) -sin(φ) ⎤
# ⎢sin(θ)⋅cos(φ) ───────────── ────────⎥
# ⎢ r r⋅sin(θ)⎥
# ⎢ ⎥
# ⎢ sin(φ)⋅cos(θ) cos(φ) ⎥
# ⎢sin(φ)⋅sin(θ) ───────────── ────────⎥
# ⎢ r r⋅sin(θ)⎥
# ⎢ ⎥
# ⎢ -sin(θ) ⎥
# ⎢ cos(θ) ──────── 0 ⎥
# ⎣ r ⎦
pr, ptheta, pphi = symbols("p_r p_theta p_phi")
# Define (s)pherical coordinate transformation
x = r*sin(theta)*cos(phi)
y = r*sin(theta)*sin(phi)
z = r*cos(theta)
spherical_momentum_vector = Matrix(([pr, ptheta, pphi]))
px,py,pz = Jinverse * spherical_momentum_vector
# px =
# pᵩ⋅sin(φ) pₜₕₑₜₐ⋅cos(φ)⋅cos(θ)
# - ───────── + pᵣ⋅sin(θ)⋅cos(φ) + ────────────────────
# r⋅sin(θ) r
# py =
# pᵩ⋅cos(φ) pₜₕₑₜₐ⋅sin(φ)⋅cos(θ)
# ───────── + pᵣ⋅sin(φ)⋅sin(θ) + ────────────────────
# r⋅sin(θ) r
# pz =
# pₜₕₑₜₐ⋅sin(θ)
# pᵣ⋅cos(θ) - ─────────────
# r
Lx = (y*pz - z*py).simplify()
Ly = (z*px - x*pz).simplify()
Lz = (x*py - y*px).simplify()
# Lx =
# pᵩ⋅cos(φ)
# - ───────── - pₜₕₑₜₐ⋅sin(φ)
# tan(θ)
# Ly =
# pᵩ⋅sin(φ)
# - ───────── + pₜₕₑₜₐ⋅cos(φ)
# tan(θ)
# Lz = pᵩ
L2 = (Lx**2 + Ly**2 + Lz**2).simplify()
# 2
# pᵩ 2
# L2 = ─────── + pₜₕₑₜₐ
# 2
# sin (θ)
In quantum mechanics we make the substitution \( p_x \mapsto -i\hbar \frac{d}{dx} \). This gives us the quantum mechanical angular momentum operators \begin{align} \hat{L}_x &= -i\hbar \left(-\sin(\phi) \frac{\partial}{\partial \theta} - \cot(\theta) \cos(\phi) \frac{\partial}{\partial \phi} \right) \\ \hat{L}_y &= -i\hbar \left( \cos(\phi) \frac{\partial}{\partial \theta} - \cot(\theta) \sin(\phi) \frac{\partial}{\partial \phi} \right) \\ \hat{L}_z &= -i\hbar \frac{\partial}{\partial \phi} \end{align} as found for example in Sakuri § 3.6.
Despite having these forms available they are typically not very useful without knowing the derivative \(\frac{\partial}{\partial \theta} Y_l^m\). Oftentimes it is better to simply make the replacements \begin{align} L_x \rightarrow \frac{L_+ + L_-}{2} \\ L_y \rightarrow \frac{L_+ - L_-}{2i} \end{align} and leverage the algebraic properties of the raising and lowering operators.