Return to Answer

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}$.

which in traditional form reads $2\text{Trace}(\sigma_i^{\dagger}\sigma_j) = \delta_{i,j}$. The ConjugateTranspose is needed to insure that the induced norm 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. This scalar product is really what answers the question: it pulls out the coefficients from the simplified operator expression in matrix form.

In this step, 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"]}
  ]

$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]]

{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"]}
  ]

$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[a.b - b.a]

In ladderReduce, I don't have to worry about commutation relations because the expression is already in a form where all operator products involve operators labeled by different site indices, which always commute. The replacement rule just substitutes a linear combination of $S_i^{\pm}$ for each $\sigma_i^{x/y}$, and this never introduces products of non-commuting opertators.

In this calculation, I make no assumptions about the spins on the chain having only nearest-neighbor interactions. The method is independent of this, but the sparseness of the matrices can be traced back to it. If one wanted to exploit the nearest-neighbor assumption, that would allow me to introduce another simplification: in the basis, only those products of spin matrices can occur that sit on adjacent sites.