Pull out scalars from NonCommutativeMultiply in commutator of SU2 spin algebra

Question

I am trying to work on spin chains in Mathematica. That is, N spin-1/2 on a one-dimensional lattice with periodic boundary condition (N+1=1). I need therefore the SU2 spin algebra $S_z, S_+, S_-$ with the commutation relations

$[S^z,S^±]=±S^±$

$[S^+,S^-]=2S^z$

My system is a product of Hilbert spaces labelled by the lattice site $j$, so that the commutation relations become:

$[S^z_i,S^±_j]=± \delta_{ij} S^±_i$

$[S^+_i,S^-_j]=2\delta_{ij}S^z_i$

I want to compute the following commutator (note the factor $j$ in the first sum):

$\left[\sum _j^3 j \left(\Delta \left(S^z _jS^z _{j+1}-\frac{1}{4}\right)+\frac{1}{2} \left(S^+_jS^-_{j+1}+S^-_jS^+_{j+1}\right)\right),\sum _k^3 \left(\Delta \left(S^z _kS^z _{k+1}-\frac{1}{4}\right)+\frac{1}{2} \left(S^+_kS^-_{k+1}+S^-_kS^+_{k+1}\right)\right)\right]$

I tried to use NonCommutativeMultiply to define a commutator but I got some problems: The factor $i$ (complex number) and the factor $\Delta$ don't get pulled out of the commutator, whereas other scalars do. Also, the commutator of two scalars doesn't give me zero and stays unevaluated. Any help on how to do this efficiently?

Answer 1

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

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

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

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]

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.

Answer 2

I will use NumericQ to test whether something is a scalar.

First, some NonCommutativeMultiply properties:

Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[H__, c_, T__] :=  c NonCommutativeMultiply[H, T] /; NumericQ[c]
NonCommutativeMultiply[H___, Times[c_, factors__], T___] :=  c NonCommutativeMultiply[
    H, Times[factors], T] /; NumericQ[c]
NonCommutativeMultiply[H___, 0, T___] := 0;
NonCommutativeMultiply[H___, 1, T___] := NonCommutativeMultiply[H, T]
NonCommutativeMultiply[H___, Plus[a_, addends__], T___] :=  NonCommutativeMultiply[H, 
    a, T] + NonCommutativeMultiply[H, Plus[addends], T]

Now, commutators.

Antisymmetry property:

Commutator[A_, A_] := 0;
Commutator /: Plus[Commutator[A_, B_], Commutator[B_, A_]] := 0;

Scalars commute with everything:

Commutator[A_, c_] := 0 /; NumericQ[c];
Commutator[c_, B_] := 0 /; NumericQ[c];

Linearity:

Commutator[Plus[a_, addends__], B_] := Commutator[a, B] + Commutator[Plus[addends], B]
Commutator[A_, Plus[b_, baddends__]] := Commutator[A, b] + Commutator[A, Plus[baddends]]
Commutator[Times[c_, factors__], B_] := c Commutator[Times[factors], B] /; NumericQ[c];
Commutator[A_, Times[c_, factors__]] := c Commutator[A, Times[factors]] /; NumericQ[c];

Leibniz law:

Commutator[NonCommutativeMultiply[A_, B_], C_] := A ** Commutator[B, C] + 
    Commutator[A, C] ** B
Commutator[A_, NonCommutativeMultiply[B_, C_]] := B ** Commutator[A, C] + 
    Commutator[A, B] ** C

Specifics for your problem:

\[CapitalDelta] /: NumericQ[\[CapitalDelta]] = True;
numSites=4;
Commutator[S[z][i_], S[plus][j_]] := KroneckerDelta[i, j] S[plus][i]
Commutator[S[z][i_], S[minus][j_]] := -KroneckerDelta[i, j] S[minus][i]
Commutator[S[plus][j_], S[z][i_]] := -KroneckerDelta[i, j] S[plus][i]
Commutator[S[minus][j_], S[z][i_]] := KroneckerDelta[i, j] S[minus][i]
Commutator[S[plus][i_], S[minus][j_]] := 2 KroneckerDelta[i, j] S[z][i]
Commutator[S[minus][i_], S[plus][j_]] := -2 KroneckerDelta[i, j] S[z][i]

Note that first assignment means ∆ will be treated as NumericQ and will match all the nice rules for scalars above. Also note that I suspect a typo in your third equation: there should be no ± inside the commutator! I have corrected this in the above definitions.

I think it is likely that you want all S operators at different lattice sites to commute (that is, locality). If so, you want

Commutator[S[labe1l_][i_], S[label2_][j_]] := 0 /; UnsameQ[i, j];

and you can possibly achieve some simplification by putting the operators into a canonical order. If the different sites commute, you can sort by site:

NonCommutativeMultiply[S[label1_][i_], S[label2_][j_]] := NonCommutativeMultiply[
    S[label2][j], S[label1][i]] /; i > j

and you can put the operators on a given site into a particular order by using information about the commutator. For example:

NonCommutativeMultiply[S[plus][i_], S[z][i_]] := NonCommutativeMultiply[
    S[z][i], S[plus][i]] - S[plus][i]
NonCommutativeMultiply[S[plus][i_], S[minus][i_]] := NonCommutativeMultiply[
    S[minus][i], S[plus][i]] + 2 S[z][i]
NonCommutativeMultiply[S[z][i_], S[minus][i_]] := NonCommutativeMultiply[
    S[minus][i], S[z][i]] - 2 S[z][i]

Optionally, you can implement periodic boundary conditions:

S[label_][i_] := S[label][Mod[i, numSites]] /; i > numSites - 1

Finally:

Commutator[
  Sum[j (\[CapitalDelta] (S[z][j] ** S[z][j + 1] - 1/4) + 
      1/2 (S[plus][j] ** S[minus][j + 1] + 
         S[minus][j] ** S[plus][j + 1])), {j, 0, numSites - 1}],
  Sum[(\[CapitalDelta] (S[z][k] ** S[z][k + 1] - 1/4) + 
      1/2 (S[plus][k] ** S[minus][k + 1] + 
         S[minus][k] ** S[plus][k + 1])), {k, 0, numSites - 1}]
  ] // ExpandAll

gives a monstrosity, but I think it's what you want. I will concede that it is notably more complicated than Jens' result, which makes me suspect that there might be a dramatic simplification that I'm missing. Also note that I wasn't sure if the sums start at 0 or 1. Generally, you might want a different convention for labeling the sites.