I also know that L and S commute, but I am unsure why. I've heard that it is simply because they act on difference variables, but I don't understand exactly what this means. Is there a way to show this explicitly?
Suppose we have two Hilbert spaces $H_1$ and $H_2$, an operator $A_1$ acting on $H_1$, and an operator $A_2$ acting on $H_2$. Let $H = H_1 \otimes H_2$. Then we can define $A_1$ and $A_2$ on $H$ by defining
\begin{align}
A_1(|a_1\rangle \otimes |a_2\rangle) &= A_1|a_1\rangle \otimes |a_2\rangle \\
A_2(|a_1\rangle \otimes |a_2\rangle) &= |a_1\rangle \otimes A_2|a_2\rangle
\end{align}
where $|a_1\rangle \in H_1, |a_2\rangle \in H_2$; and extending linearly to all of $H$. Then
\begin{align}
A_1 \circ A_2(|a_1\rangle \otimes |a_2\rangle) &= A_1(|a_1\rangle \otimes A_2|a_2\rangle) \\
&= A_1|a_1\rangle \otimes A_2|a_2\rangle \\
&= A_2(A_1|a_1\rangle \otimes |a_2\rangle) \\
&= A_2 \circ A_1(|a_1\rangle \otimes |a_2\rangle)
\end{align}
so the commutator vanishes on all pure tensors, and hence on all of $H$.
This is precisely the situation we have with the operators $L$ and $S$. In general, the wave function of a particle lives in a tensor product space. The spatial part of the wave function lives in one space, that of square-integrable functions on $\mathbb{R}^3$. The spin part, on the other hand, lives in a spinor space, i.e., some representation of $SU(2)$. $L$ acts on the spatial part, whereas $S$ acts on the spin part.
What are the remaining commutation relations between $J$, $J^2$, $L$, $L^2$, $S$, and $S^2$?
You should be able to work these out on your own, using the commutation and anti-commutation relations you already know, and properties of commutators and anti-commutators. For example,
$$[J_i, L_j] = [L_i + S_i, L_j] = [L_i, L_j] + [S_i, L_j] = i\hbar\epsilon_{ijk} L_k$$
Likewise:
\begin{align}
[J_i^2, L_j] &= J_i[J_i, L_j] + [J_i, L_j]J_i \\
&= J_i(i\hbar \epsilon_{ijk}L_k) + (i\hbar\epsilon_{ijk}L_k)J_i \\
&= i\hbar\epsilon_{ijk}\{J_i, L_k\} \\
&= i\hbar\epsilon_{ijk}\{L_i + S_i, L_k\} \\
&= i\hbar\epsilon_{ijk}(\{L_i, L_k\} + \{S_i, L_k\}) \\
&= i\hbar\epsilon_{ijk}(2\delta_{ik} I + 2S_i L_k) \\
&= 2i\hbar\epsilon_{ijk} S_i L_k
\end{align}
where we have used the fact that $L_i$ and $S_j$ commute, the linearity of $[,]$ and $\{,\}$, and the identity
$$[AB, C] = A[B, C] + [A, C]B$$
Best Answer
Since $L_i = \epsilon_{ijk} x_jp_k$ (operators) one has
$$ [L_i,L_j] = \epsilon_{iab}\epsilon_{jcd}[x_ap_b,x_cp_d] = \epsilon_{iab}\epsilon_{jcd}(x_a[p_b,x_c]p_d + x_c[x_a,p_d]p_b) $$
This first step relies on the following property of the commutator: [AB,CD] = A[B,CD] + [A,CD]B + C[AB,D] + [AB,C]D, and then performing the expansion again. The only terms that 'survive' are those involving the canonical conjugate variables. terms like $[x_a,x_b] =0$. So,
$$ [L_i,L_j] = \epsilon_{iab}\epsilon_{jcd}(x_a\underbrace{[p_b,x_c]}_{-i\hbar \delta_{b,c}}p_d + x_c\underbrace{[x_a,p_d]}_{i\hbar \delta_{ad}}p_b) = i\hbar \epsilon_{iab}\epsilon_{jcd}(-x_ap_d \delta_{bc} + x_cp_b\delta_{ad}) $$
Because of the definition of levi-civita tensor, you can absorb a minus sign by just permuting any two neighboring indices. Furthermore, after carying out the deltas, I like to rename $x_cp_b$ to $x_ap_d$ in the second term. This leads to
$$ [L_i,L_j] =i\hbar(\epsilon_{iab}\epsilon_{bjd} + \epsilon_{dib}\epsilon_{bja})x_ap_d$$
Keep in mind that any index apart from i and j are summed over $$[L_i,L_j] = i\hbar(\delta_{ij}\delta_{ad} - \delta_{id}\delta_{aj} + \delta_{dj}\delta_{ia} - \delta_{da}\delta_{ij})x_ap_d = i\hbar(x_ip_j - x_jp_i) = i\hbar \epsilon_{ijk}L_k $$
I suggest you work out the missing parts to understand how this levi-civita business works.