**Edit: a sparse array version for larger chains is given at the end.**
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}$ `KroneckerProduct`s, 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.
**Update: reducing memory consumption**
For larger `n`, the exponential growth of the basis size makes the computation very memory hungry and slow. But this can be alleviated by a combination of two additional tricks:
**First**, I use `SparseArray` in all the matrices. This means I had to rewrite all the definitions above. The explanations don't change, so I'll just list the modified definitions below.
**Second**, I avoid the use of `LinearSolve` completely by using the orthogonality relation for the Pauli matrices,
Table[2 Tr[ConjugateTranspose[#[[i]]].#[[j]]], {i, 1, 4}, {j, 1, 4}] &[
1/2 PauliMatrix[{0, 1, 2, 3}]] // MatrixForm
> $$\left(
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)$$
which in traditional form reads $2\text{Trace}(\sigma_i^{\dagger}\sigma_j) = \delta_{i,j}$. The `ConjugateTranspose` is needed to insure that the result is always non-negative in the vector space of the complex numbers, as required for a unitary scalar product. This means that the coefficient of any particular product of Pauli matrices in the complete expression can be obtained by a projection using this scalar product. As a result, the computation first determines the expansion in the Pauli matrices, not the ladder operators. Then, I add a function `ladderReduce` that substitutes the ladder operators $S_i^{\pm}$ for $\sigma_i^{x/y}$.
{sx, sy, sz} = Map[SparseArray, 1/2 PauliMatrix[{1, 2, 3}]];
sP = sx + I sy;
sM = sx - I sy;
sZ[n_, j_] :=
KroneckerProduct @@
Insert[SparseArray@Table[IdentityMatrix[2], {n - 1}], sz,
Mod[j, n, 1]]
sPlus[n_, j_] :=
KroneckerProduct @@
Insert[SparseArray@Table[IdentityMatrix[2], {n - 1}], sP,
Mod[j, n, 1]]
sMinus[n_, j_] :=
KroneckerProduct @@
Insert[SparseArray@Table[IdentityMatrix[2], {n - 1}], sM,
Mod[j, n, 1]]
id[n_] := KroneckerProduct @@ SparseArray@Table[IdentityMatrix[2], {n}]
singleBasis = {SparseArray[1/2 IdentityMatrix[2]], sx, sy, sz};
singleBasisSymbols = {1/2 &, Subsuperscript["σ", #, "x"] &,
Subsuperscript["σ", #, "y"] &,
Subsuperscript["σ", #, "z"] &};
ClearAll[basis];
basis[n_] :=
Function[{mat},
2 Map[Tr[ConjugateTranspose[
Apply[KroneckerProduct, singleBasis[[#]]]].mat] &,
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][matrix] basisNames[n]]]
]
ladderReduce[pauli_] :=
Simplify[pauli /. {Subsuperscript["σ", i_,
"x"] :> (Subsuperscript["S", i, "+"] +
Subsuperscript["S", i, "-"]),
Subsuperscript["σ", i_, "y"] :>
1/I (Subsuperscript["S", i, "+"] - Subsuperscript["S", i, "-"]),
Subsuperscript["σ", i_, "z"] :> Subsuperscript["S", i, "z"]}
]
This set of functions can now handle the example mentioned in the comment, where `n =7`:
$Assumptions = Δ > 0;
a = With[{n = 7},
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 = 7},
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}]];
pauliExpression=pauliReduce[comm[a,b]]
![sigmas][1]
The calculation takes a few seconds, but is much more efficient than before (with `n=7`, the first version of the code ran so long that I aborted it before it completed). Finally, convert this to the basis of ladder operators:
ladderReduce[pauliExpression]
![ladder][2]
[1]: https://i.sstatic.net/k41qa.png
[2]: https://i.sstatic.net/hCZYf.png