[proofplan] We write each angular momentum operator as $L_i=-i\hbar A_i$, where $A_i$ is the first-order rotational vector field $\sum_{a,b}\varepsilon_{iab}x_a\partial_{x_b}$. The commutator $[A_i,A_j]$ is computed directly by expanding both second-order differential expressions; the second-order terms cancel by symmetry, and the remaining first-order coefficient is simplified using the standard contraction identity for the Levi-Civita symbol. For the Hamiltonian, we compute separately the commutators of $L_i$ with the Laplacian and with the radial multiplication operator; the first vanishes because second derivatives commute and the antisymmetric Levi-Civita tensor kills the symmetric Hessian terms, while the second vanishes because the gradient of $V(|x|)$ is radial. [/proofplan]

[step:Introduce the rotational vector fields and record their domain invariance] For each $i \in \{1,2,3\}$, define the first-order differential operator $A_i:D\to C_c^\infty(\mathbb{R}^3 \setminus \{0\}; \mathbb{C})$ by \begin{align*} A_i\psi=\sum_{a,b=1}^3 \varepsilon_{iab} x_a \partial_{x_b}\psi . \end{align*} By the definition of $L_i$ in the theorem statement, $L_i=-i\hbar A_i$ on $D$. Since multiplication by the coordinate function $x_a:\mathbb{R}^3\to\mathbb{R}$ and differentiation $\partial_{x_b}$ preserve smoothness and do not enlarge support, each $A_i$ maps $D$ into $D$, and therefore each iterated expression $A_i(A_j\psi)$ and $L_i(L_j\psi)$ is defined for $\psi\in D$. [/step]

[step:Compute the commutator of the rotational vector fields]Fix $i,j\in\{1,2,3\}$ and $\psi\in D$. Let $\delta_{rs}$ denote the Kronecker delta on $\{1,2,3\}$, meaning $\delta_{rs}=1$ when $r=s$ and $\delta_{rs}=0$ when $r\ne s$. For each $a,b,c,d\in\{1,2,3\}$, the product rule gives \begin{align*} x_a\partial_{x_b}(x_c\partial_{x_d}\psi)=x_a\delta_{bc}\partial_{x_d}\psi+x_ax_c\partial_{x_b}\partial_{x_d}\psi . \end{align*} Using this identity in $A_i(A_j\psi)$ gives \begin{align*} A_i(A_j\psi)=\sum_{a,b,c,d=1}^3\varepsilon_{iab}\varepsilon_{jcd}x_a\delta_{bc}\partial_{x_d}\psi+\sum_{a,b,c,d=1}^3\varepsilon_{iab}\varepsilon_{jcd}x_ax_c\partial_{x_b}\partial_{x_d}\psi . \end{align*} Interchanging $i$ and $j$ gives the analogous expression for $A_j(A_i\psi)$. The second-order part of $[A_i,A_j]\psi$ is \begin{align*} \sum_{a,b,c,d=1}^3(\varepsilon_{iab}\varepsilon_{jcd}-\varepsilon_{jab}\varepsilon_{icd})x_ax_c\partial_{x_b}\partial_{x_d}\psi . \end{align*} This expression is zero because the coefficient is antisymmetric under the simultaneous exchange $(a,b,i)\leftrightarrow(c,d,j)$, while $x_ax_c\partial_{x_b}\partial_{x_d}\psi$ is symmetric under the same exchange, using equality of mixed partial derivatives for $\psi\in C^\infty$. Thus only the first-order part remains: \begin{align*} [A_i,A_j]\psi=\sum_{a,b,d=1}^3\varepsilon_{iab}\varepsilon_{jbd}x_a\partial_{x_d}\psi-\sum_{a,b,d=1}^3\varepsilon_{jab}\varepsilon_{ibd}x_a\partial_{x_d}\psi . \end{align*} Using the Levi-Civita contraction identity, obtained from the defining alternating rule for $\varepsilon_{iab}$, \begin{align*} \sum_{b=1}^3\varepsilon_{iab}\varepsilon_{jbd}=\delta_{id}\delta_{aj}-\delta_{ij}\delta_{ad} \end{align*} and the same identity with $i$ and $j$ exchanged, we obtain \begin{align*} [A_i,A_j]\psi=x_j\partial_{x_i}\psi-x_i\partial_{x_j}\psi . \end{align*} Finally, the identity \begin{align*} \sum_{k,a,b=1}^3\varepsilon_{ijk}\varepsilon_{kab}x_a\partial_{x_b}\psi=x_i\partial_{x_j}\psi-x_j\partial_{x_i}\psi \end{align*} shows that \begin{align*} [A_i,A_j]\psi=-\sum_{k=1}^3\varepsilon_{ijk}A_k\psi . \end{align*}[/step]

[guided]The purpose of introducing $A_i$ is to separate the algebra of rotations from the scalar factor $-i\hbar$. We compute $[A_i,A_j]$ first and restore the factor of $\hbar$ afterward. Fix $i,j\in\{1,2,3\}$ and $\psi\in D$. Let $\delta_{rs}$ denote the Kronecker delta on $\{1,2,3\}$, meaning $\delta_{rs}=1$ when $r=s$ and $\delta_{rs}=0$ when $r\ne s$. The operator $A_i$ differentiates in the direction of the vector field with components $\sum_a \varepsilon_{iab}x_a$. When $A_i$ is applied to $A_j\psi$, the derivative hits both the coordinate coefficient and the derivative of $\psi$. For each $a,b,c,d\in\{1,2,3\}$, the product rule gives \begin{align*} x_a\partial_{x_b}(x_c\partial_{x_d}\psi)=x_a\delta_{bc}\partial_{x_d}\psi+x_ax_c\partial_{x_b}\partial_{x_d}\psi . \end{align*} Therefore \begin{align*} A_i(A_j\psi)=\sum_{a,b,c,d=1}^3\varepsilon_{iab}\varepsilon_{jcd}x_a\delta_{bc}\partial_{x_d}\psi+\sum_{a,b,c,d=1}^3\varepsilon_{iab}\varepsilon_{jcd}x_ax_c\partial_{x_b}\partial_{x_d}\psi . \end{align*} The expression for $A_j(A_i\psi)$ is obtained by exchanging $i$ and $j$. We now subtract. The second-order terms vanish for a structural reason: the differential factor $x_ax_c\partial_{x_b}\partial_{x_d}\psi$ is symmetric under exchanging the pair $(a,b)$ with the pair $(c,d)$, because $x_ax_c=x_cx_a$ and $\partial_{x_b}\partial_{x_d}\psi=\partial_{x_d}\partial_{x_b}\psi$ for smooth $\psi$. The Levi-Civita coefficient in the commutator is antisymmetric under that same exchange. Hence the total second-order contribution cancels. The first-order contribution is \begin{align*} [A_i,A_j]\psi=\sum_{a,b,d=1}^3\varepsilon_{iab}\varepsilon_{jbd}x_a\partial_{x_d}\psi-\sum_{a,b,d=1}^3\varepsilon_{jab}\varepsilon_{ibd}x_a\partial_{x_d}\psi . \end{align*} We simplify the coefficients using the contraction identity for the Levi-Civita symbol, which follows from the defining alternating rule for $\varepsilon_{iab}$: \begin{align*} \sum_{b=1}^3\varepsilon_{iab}\varepsilon_{jbd}=\delta_{id}\delta_{aj}-\delta_{ij}\delta_{ad}. \end{align*} Applying this identity to the first sum and the corresponding exchanged identity to the second sum yields \begin{align*} [A_i,A_j]\psi=x_j\partial_{x_i}\psi-x_i\partial_{x_j}\psi . \end{align*} To identify this again as a rotational vector field, we use \begin{align*} \sum_{k,a,b=1}^3\varepsilon_{ijk}\varepsilon_{kab}x_a\partial_{x_b}\psi=x_i\partial_{x_j}\psi-x_j\partial_{x_i}\psi . \end{align*} Thus \begin{align*} [A_i,A_j]\psi=-\sum_{k=1}^3\varepsilon_{ijk}A_k\psi . \end{align*} This is exactly the Lie algebra computation for infinitesimal rotations in $\mathbb{R}^3$, written out at the level of differential operators.[/guided]

[step:Restore the quantum normalization to obtain the angular momentum commutator] Since $L_i=-i\hbar A_i$ and $L_j=-i\hbar A_j$, we have \begin{align*} [L_i,L_j]\psi=(-i\hbar)^2[A_i,A_j]\psi . \end{align*} Using the identity from the previous step, \begin{align*} [L_i,L_j]\psi=-\hbar^2\left(-\sum_{k=1}^3\varepsilon_{ijk}A_k\psi\right)=\hbar^2\sum_{k=1}^3\varepsilon_{ijk}A_k\psi . \end{align*} Since $L_k=-i\hbar A_k$, equivalently $A_k=(i/\hbar)L_k$, this becomes \begin{align*} [L_i,L_j]\psi=i\hbar\sum_{k=1}^3\varepsilon_{ijk}L_k\psi . \end{align*} [/step]

[step:Show that the Laplacian commutes with each angular momentum component] Fix $i\in\{1,2,3\}$ and $\psi\in D$. Since $\psi$ is smooth, all mixed partial derivatives appearing below commute. For each $a,b\in\{1,2,3\}$, applying the product rule twice gives \begin{align*} \Delta(x_a\partial_{x_b}\psi)=x_a\partial_{x_b}\Delta\psi+2\partial_{x_a}\partial_{x_b}\psi . \end{align*} Therefore \begin{align*} \Delta(A_i\psi)-A_i(\Delta\psi)=2\sum_{a,b=1}^3\varepsilon_{iab}\partial_{x_a}\partial_{x_b}\psi . \end{align*} The matrix with entries $\partial_{x_a}\partial_{x_b}\psi$ is symmetric in $a,b$, while $\varepsilon_{iab}$ is antisymmetric in $a,b$. Their contraction is zero, so \begin{align*} \Delta(A_i\psi)-A_i(\Delta\psi)=0 . \end{align*} Multiplying by $-i\hbar$ gives \begin{align*} \Delta(L_i\psi)-L_i(\Delta\psi)=0 . \end{align*} Hence the kinetic differential expression $-\frac{\hbar^2}{2m}\Delta$ commutes with $L_i$ on $D$. [/step]

[step:Show that radial multiplication commutes with each angular momentum component]Fix $i\in\{1,2,3\}$ and $\psi\in D$. Let $W:\mathbb{R}^3\setminus\{0\}\to \mathbb{R}$ be the radial potential on the punctured space, defined by \begin{align*} W(x)=V(|x|). \end{align*} Since $V\in C^1_{\mathrm{loc}}((0,\infty))$, the map $W$ is $C^1$ on $\mathbb{R}^3\setminus\{0\}$, and for each $b\in\{1,2,3\}$, \begin{align*} \partial_{x_b}W(x)=V'(|x|)\frac{x_b}{|x|}. \end{align*} Using the product rule, \begin{align*} A_i(W\psi)=\sum_{a,b=1}^3\varepsilon_{iab}x_a(\partial_{x_b}W)\psi+W A_i\psi . \end{align*} Thus \begin{align*} W A_i\psi-A_i(W\psi)=-\sum_{a,b=1}^3\varepsilon_{iab}x_a(\partial_{x_b}W)\psi . \end{align*} Substituting the derivative of $W$ gives \begin{align*} \sum_{a,b=1}^3\varepsilon_{iab}x_a(\partial_{x_b}W)\psi=\frac{V'(|x|)}{|x|}\left(\sum_{a,b=1}^3\varepsilon_{iab}x_ax_b\right)\psi . \end{align*} The tensor $x_ax_b$ is symmetric in $a,b$, while $\varepsilon_{iab}$ is antisymmetric in $a,b$, so the contraction is zero. Hence \begin{align*} W A_i\psi-A_i(W\psi)=0 . \end{align*} Since $L_i=-i\hbar A_i$, multiplication by $V(|x|)$ commutes with $L_i$ as a differential expression on $D$.[/step]

[guided]The potential term is multiplication by a function that depends only on the radius. Define $W:\mathbb{R}^3\setminus\{0\}\to\mathbb{R}$ by $W(x)=V(|x|)$. Because $V\in C^1_{\mathrm{loc}}((0,\infty))$ and the map $x\mapsto |x|$ is $C^1$ on $\mathbb{R}^3\setminus\{0\}$, the chain rule gives, for each $b\in\{1,2,3\}$, \begin{align*} \partial_{x_b}W(x)=V'(|x|)\frac{x_b}{|x|}. \end{align*} Now apply the product rule to $A_i(W\psi)$. The derivative hits both $W$ and $\psi$, so \begin{align*} A_i(W\psi)=\sum_{a,b=1}^3\varepsilon_{iab}x_a(\partial_{x_b}W)\psi+W A_i\psi . \end{align*} Hence \begin{align*} W A_i\psi-A_i(W\psi)=-\sum_{a,b=1}^3\varepsilon_{iab}x_a(\partial_{x_b}W)\psi . \end{align*} Substituting the derivative formula for $W$ gives \begin{align*} \sum_{a,b=1}^3\varepsilon_{iab}x_a(\partial_{x_b}W)\psi=\frac{V'(|x|)}{|x|}\left(\sum_{a,b=1}^3\varepsilon_{iab}x_ax_b\right)\psi . \end{align*} The expression $x_ax_b$ is symmetric in the indices $a,b$, while $\varepsilon_{iab}$ is antisymmetric in those same indices. Pairing the term indexed by $(a,b)$ with the term indexed by $(b,a)$ shows that their sum is zero. Therefore \begin{align*} W A_i\psi-A_i(W\psi)=0 . \end{align*} Multiplying this identity by the scalar factor $-i\hbar$ shows that multiplication by $V(|x|)$ commutes with $L_i=-i\hbar A_i$ as a differential expression on $D$.[/guided]

[step:Combine the kinetic and potential commutators] Let $i\in\{1,2,3\}$ and $\psi\in D$. The preceding two steps show that \begin{align*} -\frac{\hbar^2}{2m}\Delta(L_i\psi)-L_i\left(-\frac{\hbar^2}{2m}\Delta\psi\right)=0 \end{align*} and \begin{align*} V(|x|)L_i\psi-L_i(V(|x|)\psi)=0 \end{align*} as differential expressions on $\mathbb{R}^3\setminus\{0\}$. Adding these two identities gives \begin{align*} H(L_i\psi)-L_i(H\psi)=0 . \end{align*} Therefore $[H,L_i]\psi=0$ for every $\psi\in D$, which is the asserted commutation identity for the central Hamiltonian. [/step]