[proofplan]
The proof expands the commutator of $J_i$ with the quadratic operator $J^2$ using the product rule for commutators. The angular momentum commutation relations convert each first-order commutator into a linear combination of the components $J_k$. The resulting double sum contains a factor symmetric in the indices $j,k$ multiplied by the Levi-Civita symbol, which is antisymmetric in those same indices; pairing the terms indexed by $(j,k)$ and $(k,j)$ gives cancellation. The argument is purely algebraic on the common invariant domain $D$.
[/proofplan]
[step:Expand the commutator with the quadratic operator $J^2$]
Fix $i \in \{1,2,3\}$ and $\psi \in D$. Since each $J_a$ maps $D$ into $D$, every product of the operators $J_1,J_2,J_3$ appearing below is defined on $\psi$. Recall that $J^2:D\to D$ is the [linear map](/page/Linear%20Map) defined by
\begin{align*}
J^2\phi=\sum_{j=1}^{3}J_j^2\phi
\end{align*}
for every $\phi\in D$.
For linear maps $A,B,C : D \to D$, the commutator product identity on $D$ is
\begin{align*}
[AB,C]\psi = A[B,C]\psi + [A,C]B\psi .
\end{align*}
Indeed,
\begin{align*}
A[B,C]\psi + [A,C]B\psi = A(BC\psi-CB\psi)+(ACB\psi-CAB\psi).
\end{align*}
Cancelling the middle two terms gives
\begin{align*}
A[B,C]\psi + [A,C]B\psi = ABC\psi-CAB\psi = [AB,C]\psi .
\end{align*}
Applying this identity with $A=B=J_j$ and $C=J_i$, and using linearity of the commutator in its first argument, gives
\begin{align*}
[J^2,J_i]\psi = \sum_{j=1}^{3} [J_j^2,J_i]\psi .
\end{align*}
For each $j$,
\begin{align*}
[J_j^2,J_i]\psi = J_j[J_j,J_i]\psi + [J_j,J_i]J_j\psi .
\end{align*}
Therefore
\begin{align*}
[J^2,J_i]\psi = \sum_{j=1}^{3}\left(J_j[J_j,J_i]\psi + [J_j,J_i]J_j\psi\right).
\end{align*}
[guided]
Fix one component index $i \in \{1,2,3\}$ and one vector $\psi \in D$. The point of assuming $J_a : D \to D$ for each $a$ is that all iterated expressions such as $J_jJ_k\psi$, $J_jJ_kJ_i\psi$, and $J_iJ_j^2\psi$ remain meaningful as vectors in $D \subset H$. The quadratic angular momentum operator $J^2:D\to D$ is defined by
\begin{align*}
J^2\phi=\sum_{j=1}^{3}J_j^2\phi
\end{align*}
for every $\phi\in D$.
We first prove the algebraic identity that lets us differentiate a product inside a commutator. Let $A,B,C : D \to D$ be linear maps. By the definition of the commutator,
\begin{align*}
[AB,C]\psi = AB(C\psi)-C(AB\psi).
\end{align*}
We insert and subtract the same middle term $ACB\psi$:
\begin{align*}
AB(C\psi)-C(AB\psi) = ABC\psi-ACB\psi+ACB\psi-CAB\psi .
\end{align*}
Grouping the first two terms and the last two terms gives
\begin{align*}
ABC\psi-ACB\psi+ACB\psi-CAB\psi = A(BC\psi-CB\psi)+(ACB\psi-CAB\psi).
\end{align*}
Since $ACB\psi-CAB\psi=(AC-CA)B\psi=[A,C]B\psi$, this becomes the commutator product identity:
\begin{align*}
[AB,C]\psi = A[B,C]\psi + [A,C]B\psi .
\end{align*}
Now apply this identity to the product $J_j^2 = J_jJ_j$ and the operator $J_i$. For each $j \in \{1,2,3\}$, take $A=J_j$, $B=J_j$, and $C=J_i$. This gives
\begin{align*}
[J_j^2,J_i]\psi = J_j[J_j,J_i]\psi + [J_j,J_i]J_j\psi .
\end{align*}
Since $J^2$ is defined by
\begin{align*}
J^2\psi = \sum_{j=1}^{3} J_j^2\psi ,
\end{align*}
linearity of the commutator in its first entry gives
\begin{align*}
[J^2,J_i]\psi = \sum_{j=1}^{3}\left(J_j[J_j,J_i]\psi + [J_j,J_i]J_j\psi\right).
\end{align*}
This is the desired expansion: the commutator with the quadratic operator has been reduced to commutators between individual angular momentum components.
[/guided]
[/step]
[step:Substitute the angular momentum commutation relations]
Let $\hbar\in\mathbb{R}$ denote the reduced Planck constant, and let $\varepsilon_{abc}$ denote the Levi-Civita symbol on $\{1,2,3\}$, so $\varepsilon_{abc}=1$ for even permutations of $(1,2,3)$, $\varepsilon_{abc}=-1$ for odd permutations, and $\varepsilon_{abc}=0$ when two indices coincide. By the defining angular momentum commutation relations assumed in the theorem statement, for each $j \in \{1,2,3\}$ the relation with indices $j,i$ gives
\begin{align*}
[J_j,J_i]\psi = i\hbar \sum_{k=1}^{3} \varepsilon_{jik}J_k\psi .
\end{align*}
Because $J_j$ is linear and maps $D$ into $D$,
\begin{align*}
J_j[J_j,J_i]\psi = i\hbar \sum_{k=1}^{3} \varepsilon_{jik}J_jJ_k\psi .
\end{align*}
Also, since $J_j\psi \in D$, applying the same commutation relation to the vector $J_j\psi$ gives
\begin{align*}
[J_j,J_i]J_j\psi = i\hbar \sum_{k=1}^{3} \varepsilon_{jik}J_kJ_j\psi .
\end{align*}
Substituting both identities into the expansion from the previous step yields
\begin{align*}
[J^2,J_i]\psi = i\hbar \sum_{j=1}^{3}\sum_{k=1}^{3}\varepsilon_{jik}(J_jJ_k+J_kJ_j)\psi .
\end{align*}
[guided]
We now replace the abstract commutators in the expansion by the angular momentum relations. Here $\hbar\in\mathbb{R}$ is the reduced Planck constant, and $\varepsilon_{abc}$ is the Levi-Civita symbol on $\{1,2,3\}$: it equals $1$ on even permutations of $(1,2,3)$, equals $-1$ on odd permutations, and equals $0$ when two indices coincide. The theorem statement assumes that these angular momentum commutation relations hold on the common domain $D$. Thus, for each $j \in \{1,2,3\}$, the relation with ordered indices $j,i$ states that, for the fixed vector $\psi \in D$,
\begin{align*}
[J_j,J_i]\psi = i\hbar \sum_{k=1}^{3} \varepsilon_{jik}J_k\psi .
\end{align*}
The domain condition is used twice here. First, because each $J_k\psi$ lies in $D$, the finite sum on the right is a vector in $D$. Second, because $J_j:D\to D$ is linear, applying $J_j$ to that finite sum gives
\begin{align*}
J_j[J_j,J_i]\psi = i\hbar \sum_{k=1}^{3}\varepsilon_{jik}J_jJ_k\psi .
\end{align*}
For the other term, the input to the commutator is not $\psi$ but $J_j\psi$. This is legitimate because $J_j\psi\in D$. Applying the same angular momentum relation to this vector gives
\begin{align*}
[J_j,J_i]J_j\psi = i\hbar \sum_{k=1}^{3}\varepsilon_{jik}J_kJ_j\psi .
\end{align*}
Substituting these two identities into the expansion of $[J^2,J_i]\psi$ gives
\begin{align*}
[J^2,J_i]\psi = i\hbar \sum_{j=1}^{3}\sum_{k=1}^{3}\varepsilon_{jik}(J_jJ_k+J_kJ_j)\psi .
\end{align*}
The proof has now reduced the desired commutation statement to a finite algebraic cancellation indexed by the two component labels $j$ and $k$.
[/guided]
[/step]
[step:Cancel the symmetric operator factor against the antisymmetric Levi-Civita symbol]
For fixed $j,k \in \{1,2,3\}$, define the operator $S_{jk}:D \to D$ by
\begin{align*}
S_{jk}(\phi) = (J_jJ_k+J_kJ_j)\phi
\end{align*}
for every $\phi \in D$. Then $S_{jk}=S_{kj}$ by definition. The Levi-Civita symbol is antisymmetric in the first and third displayed indices here:
\begin{align*}
\varepsilon_{kij} = -\varepsilon_{jik}.
\end{align*}
Thus, for every $j,k \in \{1,2,3\}$,
\begin{align*}
\varepsilon_{jik}S_{jk}\psi+\varepsilon_{kij}S_{kj}\psi = \varepsilon_{jik}S_{jk}\psi-\varepsilon_{jik}S_{jk}\psi = 0 .
\end{align*}
Terms with $j=k$ also vanish because $\varepsilon_{jij}=0$. Hence the whole double sum is zero:
\begin{align*}
\sum_{j=1}^{3}\sum_{k=1}^{3}\varepsilon_{jik}(J_jJ_k+J_kJ_j)\psi = 0 .
\end{align*}
Therefore
\begin{align*}
[J^2,J_i]\psi = 0 .
\end{align*}
Since $i \in \{1,2,3\}$ and $\psi \in D$ were arbitrary, $J^2$ commutes with each angular momentum component on $D$.
[guided]
The remaining expression is
\begin{align*}
i\hbar \sum_{j=1}^{3}\sum_{k=1}^{3}\varepsilon_{jik}(J_jJ_k+J_kJ_j)\psi .
\end{align*}
The scalar factor $i\hbar$ is irrelevant for vanishing, so we focus on the double sum. For fixed $j,k \in \{1,2,3\}$, define $S_{jk}:D\to D$ by
\begin{align*}
S_{jk}(\phi) = (J_jJ_k+J_kJ_j)\phi
\end{align*}
for every $\phi\in D$. This operator is symmetric in $j$ and $k$: exchanging the two labels gives
\begin{align*}
S_{kj}(\phi) = (J_kJ_j+J_jJ_k)\phi = S_{jk}(\phi)
\end{align*}
for every $\phi\in D$.
The Levi-Civita symbol has the opposite behaviour: it is antisymmetric when two indices are exchanged. With the middle index $i$ fixed, swapping $j$ and $k$ gives
\begin{align*}
\varepsilon_{kij} = -\varepsilon_{jik}.
\end{align*}
Therefore the two terms in the double sum indexed by $(j,k)$ and $(k,j)$ cancel:
\begin{align*}
\varepsilon_{jik}S_{jk}\psi+\varepsilon_{kij}S_{kj}\psi = \varepsilon_{jik}S_{jk}\psi-\varepsilon_{jik}S_{jk}\psi = 0 .
\end{align*}
When $j=k$, the term vanishes separately because the Levi-Civita symbol is zero whenever two indices coincide:
\begin{align*}
\varepsilon_{jij}=0.
\end{align*}
Thus every term in the finite double sum is either paired with an equal and opposite term or is zero on its own. Hence
\begin{align*}
\sum_{j=1}^{3}\sum_{k=1}^{3}\varepsilon_{jik}(J_jJ_k+J_kJ_j)\psi = 0 .
\end{align*}
Substituting this into the expression for $[J^2,J_i]\psi$ gives
\begin{align*}
[J^2,J_i]\psi = 0 .
\end{align*}
Because the component index $i\in\{1,2,3\}$ and vector $\psi\in D$ were arbitrary, this proves that $J^2$ commutes with each angular momentum component on the common invariant domain $D$.
[/guided]
[/step]
