For a specific quantum mechanical problem I need to multiply out operators in order to calculate a trace by hand. For example I need a Hamiltonian squared with $H^2$. The Hamiltonian contains of a few single terms and flip-flop terms embedded in sums like $\frac{1}{2} \sum_{m=1}^N S_0^+ S_m^- + S_0^- S_m^+$.
My first idea was simply to write down the sums and ask Mathematica to multiply them out by expanding them. Unfortunately, this yields a result I could have imagined myself with $\sum_m (...) \times \sum_n (...)$. What I need is a sum of those product terms for example of the form $\sum_{m=1,n=1}^N S_0^+ S_m^-S_0^+ S_n^- + \sum_{m=1,n=1}^N (...)$ as I want to evaluate each operator separately by hand.
Changing the sum parenthesis doesn't seem to have any effect. Expand refuses to work.
Is there a possibility to let Mathematica multiply out all those terms individually?
Remark: As specific operators don't commute, e.g. $[S^+,S^-]$ for the same index I need Mathematica to give out each term. Simplifying by changing operator positions could lead to false results, so something like $(a+b)*(a+b)=a^2+ab+ba+b^2$ is needed here.
Edit: Here is a minimum example of what doesn't work as expected:
Expand[(-h Subscript[S^z, 0] -
Sum[Subscript[J,
m] (Subscript[S, 0]^z Subscript[S, m]^z +
1/2 SuperPlus[Subscript[S, 0]] SuperMinus[Subscript[S, m]] +
1/2 SuperMinus[Subscript[S, 0]] SuperPlus[Subscript[S,
m]]), {m, 1, NN}]) (-h Subscript[S^z, 0] -
Sum[Subscript[J,
n] (Subscript[S, 0]^z Subscript[S, n]^z +
1/2 SuperPlus[Subscript[S, 0]] SuperMinus[Subscript[S, n]] +
1/2 SuperMinus[Subscript[S, 0]] SuperPlus[Subscript[S,
n]]), {n, 1, NN}])]
NN
known, or does it need to remain asNN
? $\endgroup$