Normal-ordering bosonic operators

F.P. Laussy
Russian Quantum Center
21 May 2022

Introduction

The implementation relies on a “weighted sequence” (wseq), which is of the type

{w,{k,l,m,n}}

where w is the weight and k, l, m, n, etc. (the list can be of any size) are the exponents of a products of annihilation and creation operators, always starting with an annihilation one on the right (normal order), so, e.g.,

{1,{1,1}}

is (SuperDagger[a])a while

{2,{1,1,0}}

is 2a (SuperDagger[a]). Of course one can make more complex sequences, e.g.,

{-1,{1,2,3,1}}

is -(SuperDagger[a])(a^2) (a^(†3))a. Various such weighted sequences in a list are understood as a sum of them, e.g.

{{-1,{1,2,3,1}},{2,{1,1,0}}}

is -(SuperDagger[a])(a^2) (a^(†3))a + 2 a (SuperDagger[a]). Etc. What this notebook does is, given an expression in any order, it returns it normally-ordered version, e.g.:

NormalOrder[{1,{1,1,0}}]

This is the basic a (SuperDagger[a])=1+(SuperDagger[a])a

NormalOrder[{1,{2,2,0}}]

This is (a^2)(a^(†2))=2+4(SuperDagger[a])a +(a^(†2))(a^2) , etc.

Code

There’s an error generated when turning operators into constant, but it’s okay just to ignore it:

Off[Part::partw]

Drop leading zeros from a wseq: this effectively trim out the 0th power of operators.

DropLeadingZeros[wseq_]:=Module[{cz},cz=1;
While[wseq[[2,cz]]==0,cz++];
{wseq[[1]],Drop[wseq[[2]],cz-1]/. {}->{0}}]

Normal order wseqs which first element is not normally ordered, i.e., starting with (a^(†n))

FirstNotOrdered[myseq_]:=Module[{wseq,mseq,nseq,minseq,newseq},
wseq=myseq[[1]];
nseq=myseq[[2,1]];
mseq=myseq[[2,2]];
minseq=Min[mseq,nseq];
newseq=Table[{wseq,myseq[[2]]},{k,0,minseq}];
Do[
newseq[[k+1]][[1]]=k!wseq Binomial[mseq,k]Binomial[nseq,k];
newseq[[k+1]][[2,1]]=0;
newseq[[k+1]][[2,2]]=mseq-k;
newseq[[k+1]][[2,3]]+=nseq-k;
,
{k,0,minseq}];
newseq
]

FirstOrdered[myseq_]:=Module[{wseq,mseq,nseq,minseq,newseq},
wseq=myseq[[1]];
nseq=myseq[[2,2]];
mseq=myseq[[2,3]];
minseq=Min[mseq,nseq];
newseq=Table[{wseq,myseq[[2]]},{k,0,minseq}];
Do[
newseq[[k+1]][[1]]=k!wseq Binomial[mseq,k]Binomial[nseq,k];
newseq[[k+1]][[2,2]]=0;
newseq[[k+1]][[2,3]]=mseq-k;
newseq[[k+1]][[2,4]]+=nseq-k;
,
{k,0,minseq}];
newseq
]

Contract[myseq_]:=Module[{},
{myseq[[1]],Flatten[{Total[myseq[[2,{1,3}]]],myseq[[2]][[Range[4,Length[myseq[[2]]]]]]}]}
]

NormalOrderThroughOnce[seqinput_]:=Module[{remainseq,remainseqother,doneseq,leadseq,newseq,num},
(* Entries that can be processed without normal-ordering *)
(* This should apply only on 1st call if input is not optimum *)
remainseq=seqinput[[1]];
doneseq=seqinput[[2]];
doneseq=Join[doneseq,Select[remainseq,Length[#[[2]]]<=2&]];
remainseq=Select[remainseq,Length[#[[2]]]>2&];
remainseq=DropLeadingZeros/@remainseq;
(* Working first on seq whose first-element is not-normally ordered, i.e., starting with (a^(†n)) *)
leadseq=Select[remainseq,OddQ[Length[#[[2]]]]&];
remainseq=Select[remainseq,EvenQ[Length[#[[2]]]]&];
newseq=Flatten[FirstNotOrdered/@leadseq,1];
newseq=DropLeadingZeros/@newseq;
doneseq=Join[doneseq,Select[newseq,Length[#[[2]]]<=2&]];
remainseqother=Select[Join[remainseq,newseq],Length[#[[2]]]>2 && EvenQ[Length[#[[2]]]]&];
remainseq=Select[Join[remainseq,newseq],Length[#[[2]]]>2 && OddQ[Length[#[[2]]]]&];
(* Working then on seq whose first-element is a dagger operator*)
newseq=Flatten[FirstOrdered/@remainseqother,1];
(* merge the terms that are settled *)
num[{n_,m_}]=0;
num[{n_}]=0;
(num[#[[2]]]+=#[[1]])&/@doneseq;
(* That's the result after one iteration of the normal-ordering *)
{Join[remainseq,Contract/@newseq],Table[{num[i],i},{i,Union[doneseq[[All,2]]]}]}
]

NormalOrder[lseq_]:=Module[{fullseq,doneseq},
fullseq={lseq};doneseq={};
Select[FixedPoint[NormalOrderThroughOnce,{fullseq,{}}][[2]],#[[1]]!=0&]
]

Examples

Remember that encoding is such that the rightmost index is the exponent for a, then next is (SuperDagger[a]), then a, (SuperDagger[a]), a,(SuperDagger[a]), etc.

NormalOrder[{1,{1,1,0}}]

NormalOrder[{1,{1,5,0}}]

NormalOrder[{1,{1,1,1,1,0}}]

The Stirling numbers:

NormalOrder[{1,{1,1,1,1,1,1,1,1,0}}]

Checking analytical formulaes:

NormalOrder[{1,{7,5,0}}]

m=7;n=5;

Table[{k!Binomial[m,k]Binomial[n,k],{m-k,n-k}},{k,0,Min[n,m]}]

NormalOrder[{1,{2,2,2,2,2,2,0}}]

NormalOrder[{1,{3,3,3,3,3,3,0}}]