Since you're specifically asking about $SU(2)$, we have a convenient representation in the form of the $2\times 2$ Pauli matrices. This makes it seem very beneficial to map your problem onto a matrix multiplication. Here is how you can do it.
First I define the three single spin matrices, then embed them into the product space of the $N$ spin system using the functions sZ, sPlus and sMinus which take the dimension N as the first argument and the site index j as the second argument.
There is an ambiguity here in that we can come up with different ways of defining the larger matrices satisfying the given commutation relations. But I will use some additional physics to assume that we are working in the $2^N$ dimensional space of spins of the sites. This means I will use KroneckerProduct to construct matrices that act like the identity on all but one spin, and there will be an identity of dimension $2^N$ that multiplies the factors $1/4$ in your expression, too. Adding this, I get the following:
{sx, sy, sz} = 1/2 PauliMatrix[{1, 2, 3}];
sP = sx + I sy;
sM = sx - I sy;
sZ[n_, j_] :=
KroneckerProduct @@
Insert[Table[IdentityMatrix[2], {n - 1}], sz, Mod[j, n, 1]]
sPlus[n_, j_] :=
KroneckerProduct @@
Insert[Table[IdentityMatrix[2], {n - 1}], sP, Mod[j, n, 1]]
sMinus[n_, j_] :=
KroneckerProduct @@
Insert[Table[IdentityMatrix[2], {n - 1}], sM, Mod[j, n, 1]]
id[n_] := KroneckerProduct @@ Table[IdentityMatrix[2], {n}]
Now verify for an example (site 1 in a 2-spin system) that the commutation relations hold, including the boundary condition that maps j=3 to j=1:
sZ[2, 3].sPlus[2, 1] - sPlus[2, 1].sZ[2, 1] == sPlus[2, 1]
(* ==> True *)
sZ[2, 1].sMinus[2, 1] - sMinus[2, 1].sZ[2, 1] == -sMinus[2, 1]
(* ==> True *)
sPlus[2, 1].sMinus[2, 1] - sMinus[2, 1].sPlus[2, 1] ==
2 sZ[2, 1]
(* ==> True *)
Finally, define the commutator and the two operators of interest:
comm[a_, b_] := a.b - b.a
a = With[{n = 3},
Sum[j (Δ (sZ[n, j].sZ[n, j + 1] -
1/4 id[n]) +
1/2 (sPlus[n, j].sMinus[n, j + 1] +
sMinus[n, j].sPlus[n, j + 1])), {j, 1, n}]];
b = With[{n = 3},
Sum[Δ (sZ[n, j].sZ[n, j + 1] -
1/4 id[n]) +
1/2 (sPlus[n, j].sMinus[n, j + 1] +
sMinus[n, j].sPlus[n, j + 1]), {j, 1, n}]];
comm[a, b]
(*
==> {{0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1/2 - Δ/2, 0,
1/4 - Δ/4, 0, 0, 0}, {0, -(1/2) + Δ/2,
0, 0, -(1/4) + Δ/4, 0, 0, 0}, {0, 0, 0, 0, 0,
1/4 - Δ/4, -(1/4) + Δ/4,
0}, {0, -(1/4) + Δ/4, 1/4 - Δ/4, 0, 0,
0, 0, 0}, {0, 0, 0, -(1/4) + Δ/4, 0,
0, -(1/2) + Δ/2, 0}, {0, 0, 0,
1/4 - Δ/4, 0, 1/2 - Δ/2, 0, 0}, {0, 0,
0, 0, 0, 0, 0, 0}}
*)
This is the result for the commutator.
Edit
Now the goal of the question is to get an algebraic result, not a big matrix. However, this can also be achieved within the above matrix method, by reversing the translation from operator language to matrices. The idea is that a basis of the single-spin $2\times 2$ matrix space is formed by the four elements of singleBasis defined below. From this, we can construct a basis of the $2^N$ dimensional matrix space acting on the product of $n$ spin spaces, again using KroneckerProduct. This works because there are $4$ elements to choose from for each of the $N$ spins, giving $4^N = 2^{2N}$ KroneckerProducts, just as many as there are entries in a $2^N$ dimensional square matrix.
The central ingredient is the function basis which creates a LinearSolveFunction based on the KroneckerProduct basis matrices just described. I use memoization to remember the calculation once it has been called for a particular dimension $N$, because presumably one wants to simplify more than one expression with the same $N$. The linear solver is invoked here to find the expansion coefficients of a given matrix like the output from the commutator above, in the basis of the product matrix space.
The next step is then done in pauliReduce, which takes a specific more or less complicated large matrix coming from a calculation like the one above, and then uses the basis expansion coefficients to rewrite that matrix as a compact symbolic sum over the spin operators.
The output uses the function basisNames which creates a list of symbols in the same order as the list of KroneckerProduct basis matrices. The formatted set of names for the spin operators is defined in singleBasisSymbols. The output uses regular sums and commutative products because after the simplifications each term in the sum will have only one operator for each site, and the operators for different sites commute anyway.
singleBasis = {IdentityMatrix[2], sP, sM, sz};
singleBasisSymbols = {1 &, \!\(
\*SubsuperscriptBox[\("\<S\>"\), \(#\), \("\<+\>"\)] &\), \!\(
\*SubsuperscriptBox[\("\<S\>"\), \(#\), \("\<-\>"\)] &\), \!\(
\*SubsuperscriptBox[\("\<S\>"\), \(#\), \("\<z\>"\)] &\)};
ClearAll[basis];
basis[n_] :=
basis[n] =
LinearSolve@
Transpose[
Map[Flatten,
KroneckerProduct @@ singleBasis[[#]] & /@
Tuples[Table[Range[1, 4], {n}]]]]
names[indexList_List] :=
Times @@ MapIndexed[singleBasisSymbols[[#]][#2[[1]]] &, indexList]
basisNames[n_] := names /@ Tuples[Table[Range[1, 4], {n}]]
Clear[pauliReduce];
pauliReduce[matrix_?MatrixQ] := Module[{n = Log[2, Length[matrix]]},
Simplify[Total[basis[n][Flatten[matrix]] basisNames[n]]]
]
Next, do some tests, among them some commutation relations. Then apply it to the commutator expression in the question:
pauliReduce[sZ[3, 1].sMinus[3, 1]]
$-\frac{\text{S}_1^-}{2}$
pauliReduce[comm[sPlus[3, 1], sMinus[3, 1]]]
$2 \text{S}_1^{\text{z}}$
pauliReduce[comm[sPlus[3, 2], sMinus[3, 1]]]
$0$
pauliReduce[comm[a, b]]
$$-\frac{1}{2} (\Delta -1) \left(\text{S}_1^- \left(\text{S}_2^+ \text{S}_3^{\text{z}}-\text{S}_3^+ \text{S}_2^{\text{z}}\right)+\text{S}_1^+ \left(\text{S}_3^- \text{S}_2^{\text{z}}-\text{S}_2^- \text{S}_3^{\text{z}}\right)+2 \text{S}_1^{\text{z}} \left(\text{S}_3^- \text{S}_2^+-\text{S}_2^- \text{S}_3^+\right)\right)$$
This is the most reduced form in which the commutator expression can be written, and I didn't need to define any commutator algebra. The only price to pay is that the matrix space whose basis I have to keep track of becomes exponentially large with increasing number of sites, $N$. The benefit is that all the algebra is done as matrix manipulations and hence requires no additional thought per se.