Sorting function for non commuting bosons

Question

I am trying to write a sorting function which will sort expressions involving products of bosonic objects which do not commute. For example, I can have objects like $a,\ a^\dagger,\ b,\ b^\dagger$ where $[a,a^\dagger] = 1$, $[b,b^\dagger] = 1$ and all other commutators vanish. So the sorting function should for example sort a term like $a a^\dagger b a b^\dagger$ to $a + a^\dagger a^2 + ab^\dagger b+a^\dagger a^2 b^\dagger b$. So I want to sort them first lexicographically and then such that dagger objects will be to the left.

I have been using a function from the notebook provided by Simon Tyler at arXiv. So I have created a bosonic object using

Clear[Boson, BosonC, BosonA]
Boson /: MakeBoxes[
 Boson[cr : (True | False), sym_], fmt_] := 
  With[{sbox = If[StringQ[sym], sym, ToBoxes[sym]]}, 
   With[{abox = If[cr, SuperscriptBox[#, FormBox["\[Dagger]", Bold]], #] &@sbox},
    InterpretationBox[abox, Boson[cr, sym]]
   ]
  ]
BosonA[sym_:"a"] := Boson[False, sym]
BosonC[sym_:"a"] := Boson[True, sym]

Next following the code, I want to turn the atomic expressions into lists and use Mathematica's ordering after that. So I have functions like

NCOrder[Boson[cr : True | False, sym_]] := 
NCOrder[Boson[cr, sym]] = {{sym, Boole[cr]}}

NCSort[expr_] := Module[{temp, NCM},
 temp = expr /. NonCommutativeMultiply -> NCM;   
 temp = temp /. (NCM[a : _Boson ..] :> (NCM[a][[Ordering[#]]] &[NCOrder /@ {a}]));
 temp = temp /. NCM[singleArg_] :> singleArg;   
 temp = temp /. NCM -> NonCommutativeMultiply
]

This NCSort function right now is just ordering it lexicographically. I want to modify it so that it would sort and order the expressions as I explained above.

Ordering of creation/annihilation operators isn't a sort in a traditional sense as you are not just rearranging terms, but modifying them in the process. Most sorting assumes a stable length, and this process does not have that. — rcollyer, Dec 13 '12 at 19:41
See section "Some noncommutative algebraic manipulation" in this nb. It has a couple of examples that use commutators and do canonicalization. — Daniel Lichtblau, Dec 13 '12 at 20:51

Leonid Shifrin · Accepted Answer · 2012-12-14T12:36:03.663

9

I suggest the following function, which looks rather ugly and is not really efficient, but seems to do the job:

expand[expr_] :=
 Module[{times},
   expr /.
     NonCommutativeMultiply[b___] :> times[b] //.
     {
       times[left___, cnum_ /; FreeQ[cnum, _Boson], right___] :> 
             cnum*times[left, right],
       times[left___, fst : Boson[_, s_], middle__, sec : Boson[_, s_],right___] /;
          FreeQ[{middle}, Boson[_, s]] &&
              OrderedQ[{Last@Union[Cases[{middle}, Boson[_, sym_] :> sym]], s}] :>
                  times[left, middle, fst, sec, right],
       times[left___, Boson[False, s_], Boson[True, s_], right___] :>
          times[left, right] + times[left, Boson[True, s], Boson[False, s], right],
       times[b_Boson] :> b
     } //.
    times[left___, p : Longest[PatternSequence[(b : Boson[type_, s_]) ..]],right___] :>
      times[left, b^Length[{p}], right] /. 
    times[arg__] :> NonCommutativeMultiply[arg] /. times[] -> OverHat[1]
];

You can use it as

expand[NonCommutativeMultiply[BosonA["a"], BosonC["a"], BosonA["b"], BosonA["a"], BosonC["b"]]]

(*
    a+a^\[Dagger]**a^2+a**b^\[Dagger]**b+a^\[Dagger]**a^2**b^\[Dagger]**b
*)

edited Dec 14 '12 at 12:36

answered Dec 13 '12 at 21:25

Leonid Shifrin

114,335
15
329
420

@Jens Works for me expand[NonCommutativeMultiply[BosonA["a"],BosonC["a"]]] gives NonCommutativeMultiply[]+a^\[Dagger]**a – Leonid Shifrin Dec 13 '12 at 21:52
@Jens Right, I forgot to mention - should be If instead of f. I corrected it in the post. – Leonid Shifrin Dec 13 '12 at 22:05
Thanks - I just saw that it had already been edited but I must have had the old definition. Miscommunication. – Jens Dec 13 '12 at 22:11
@Jens No, I edited already after your comment on this, so you could not have used the edited version when you tried. – Leonid Shifrin Dec 13 '12 at 22:14
Good, so I was only half blind. Anyway, I already gave you +1 and I've tested some more, also combining with Plus, without finding anything wrong. – Jens Dec 13 '12 at 22:22
@Jens Thanks for testing and an upvote, good to know that it works. It basically does what we usually do by hand. But, if you ask me, it is esier to switch to path integrals, no operators - no problems :) (I am joking - operator formalism certainly has its merits) – Leonid Shifrin Dec 13 '12 at 22:29
@Leonid I just saw this and I am still trying to understand the code but I tried it and it works great for the example that I mentioned however I tried another example like a**a**a^\dagger**a^\dagger which should expand to 2 + 4 a^\dagger**a + a^\dagger**a^\dagger*a*a but it doesn't expand anything. – Bilentor Dec 13 '12 at 22:29
@Karan Your example does expand exactly as expected for me. – Jens Dec 13 '12 at 22:40
@Karan Works for me too: expand[NonCommutativeMultiply[BosonA["a"], BosonA["a"], BosonC["a"], BosonC["a"]]] results in 2 NonCommutativeMultiply[] + 4 Boson[True, "a"] ** Boson[False, "a"] + (Boson[True, "a"])^2 ** Boson[False, "a"]^2. Make sure you use the corrected definitions for the boson operator formatting: I edited your post, since you had f instead of If in one place. – Leonid Shifrin Dec 13 '12 at 22:44
1

Leonid, here's a suggestion for improvement: how about replacing the last rule in your expand module by times[arg__] :> NonCommutativeMultiply[arg] /. times[] -> OverHat[1] so that the identity operator is displayed in a more readable form? – Jens Dec 13 '12 at 22:50
It works now. I had modified the properties of NonCommutativeMultiply in my original notebook to make it distributive and other things which I believe was conflicting somehow with expand. – Bilentor Dec 13 '12 at 22:54
@Jens Thanks, edited. – Leonid Shifrin Dec 13 '12 at 22:58
So I ran into a trouble. The expressions I need to simplify contain c-numbers in front of the oscillators. The function expand will only simplify if the input arguments are purely Boson objects. Is there to modify the function so that it will factor out the c-numbers and then expand? – Bilentor Dec 14 '12 at 03:52
@Karan See the edited version. Try to understand the logic here, this part with c-numbers was not hard at all. – Leonid Shifrin Dec 14 '12 at 12:36
@Leonid I tried to use this times[left___, Times[n_?NumericQ, b : Boson[type_, s_]], right___] :> n*times[left, b, right]. Can you tell me what is wrong with this line? – Bilentor Dec 14 '12 at 13:44
@Leonid Also the edited version evaluates a**2**b to 2 a**b but not a**(2*b). – Bilentor Dec 14 '12 at 13:54
@Karan I don't have much time now, but Times sometimes evaluates, so you often need HoldPattern. I also did not include Times in the rule, that's why your input does not evaluate. – Leonid Shifrin Dec 14 '12 at 14:05
@Leonid I appended this line in the function times[left___, cnum_ /; (! FreeQ[cnum, Times[n_?NumericQ, _Boson]]), right___] :> Extract[FactorList[cnum], {1, 1}]*times[left,Extract[FactorList[cnum], {2, 1}], right] and it does the job however there are some other issues. When I have fractions, it doesn't work. Also the function isn't distributive over addition. I tried a similar command but that is not very general. – Bilentor Dec 14 '12 at 20:29
@Karan Good that you managed. In general, such problems are expected. If you want to make it rather general-purpose, you have to add more rules. – Leonid Shifrin Dec 14 '12 at 20:32

Sorting function for non commuting bosons

1 Answers1

Linked