Normal-ordering bosonic operators

F.P. Laussy
Russian Quantum Center
21 May 2022
Revised 31 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}}}]

If there is a single correlator, then one can use instead:

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

NormalOrderCorrelator[{1,1,0}]

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

NormalOrderCorrelator[{2,2,0}]

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,k},
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,k},
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,n,m},
(* 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]]]}]}
]

NormalOrderCorrelator[seq_]:=Module[{},
Select[FixedPoint[NormalOrderThroughOnce,{{{1,seq}},{}}][[2]],#[[1]]!=0&]
]

NormalOrder[lseq_]:=Module[{},
Select[FixedPoint[NormalOrderThroughOnce,{lseq,{}}][[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.

This finds that 2(SuperDagger[a])a (a^(†3))(a^2)(SuperDagger[a])a-(SuperDagger[a])a(SuperDagger[a])a(a^(†3))(a^2)-(a^(†3))(a^2)(SuperDagger[a])a(SuperDagger[a])a=-(a^(†3))(a^2):

NormalOrder[{
{2,{1,1,3,2,1,1}},
{-1,{1,1,1,1,3,2}},
{-1,{3,2,1,1,1,1}}
}]

This finds that a (SuperDagger[a])=1+(SuperDagger[a])a using the two main functions:

NormalOrderCorrelator[{1,1,0}]

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

NormalOrderCorrelator[{1,5,0}]

NormalOrderCorrelator[{1,1,1,1,0}]

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

The Stirling numbers:

NormalOrderCorrelator[{1,1,1,1,1,1,1,1,0}]

Checking analytical formulas:

k=.;l=.;

Flatten[Table[Sort[NormalOrderCorrelator[{1,1,k,l,1,1}]]==Sort[{{1,{k+2,l+2}},{k+l+1,{k+1,l+1}},{k l,{k,l}}}],{k,5},{l,3}]]//Union

NormalOrder[{
{2,{1,1,3,2,1,1}},
{-1,{1,1,1,1,3,2}},
{-1,{3,2,1,1,1,1}}
}]

k=2;l=3;

NormalOrder[{
{2,{1,1,k,l,1,1}},
{-1,{k,l,1,1,1,1}},
{-1,{1,1,1,1,k,l}}
}]

Flatten[Table[Sort[NormalOrderCorrelator[{1,1,1,1,k,l}]]==Sort[{{1,{k+2,l+2}},{1+2k,{k+1,l+1}},{k+k(k-1),{k,l}}}],{k,5},{l,3}]]//Union

Flatten[Table[Sort[NormalOrderCorrelator[{k,l,1,1,1,1}]]==Sort[{{1,{k+2,l+2}},{1+2l,{k+1,l+1}},{l+l(l-1),{k,l}}}],{k,5},{l,3}]]//Union

Table[{Sort[NormalOrderCorrelator[{1,1,k,l,1,1}]],Sort[{{1,{k+2,l+2}},{k+l+1,{k+1,l+1}},{k l,{k,l}}}]},{k,0,5},{l,3}]

Table[Sort[NormalOrderCorrelator[{1,1,k,l,1,1}]]==Sort[{{1,{k+2,l+2}},{k+l+1,{k+1,l+1}},{k l,{k,l}}}],{k,0,5},{l,3}]

Sort[NormalOrderCorrelator[{1,1,1,1,1,1}]]

Sort[{{1,{k+2,l+2}},{k+l+1,{k+1,l+1}}}/.{k->1,l->1}]

Sort[%]

Sort[{{1,{5,4}},{6,{3,2}},{6,{4,3}}}]

Sort[({{1,{k+2,l+2}},{l+k+1,{k+1,l+1}},{k l,{k,l}}})]

{{1,{5,4},{6,{4,3}},{6,{3,2}}}}

Sort[{{1,{5,4},{6,{4,3}},{6,{3,2}}}}]

Sort[{{1,{5,4},{6,{4,3}},{6,{3,2}}}}]

Sort[{{1,{5,4}},{6,{4,3}},{6,{3,2}}}]

Sort[%]

Union[{{1,{5,4},{6,{4,3}},{6,{3,2}}}}]

%==%%

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

Table[Sort[NormalOrder[{1,{1,k,l,1,0}}]]==DropLeadingZeros/@Sort[{{1,{k+1,l+1}},{k+l+1,{k,l}},{k l,{k-1,l-1}}}/.{{0,0}->{0}}],{k,1,6},{l,1,6}]

Table[Sort[NormalOrder[{1,{1,k+1,l}}]]==DropLeadingZeros/@Sort[{{1,{k+1,l+1}},{k+1,{k,l}}}/.{{0,0}->{0}}],{k,1,6},{l,1,6}]

Table[Sort[NormalOrder[{1,{k,l+1,1,0}}]]==DropLeadingZeros/@Sort[{{1,{k+1,l+1}},{l+1,{k,l}}}/.{{0,0}->{0}}],{k,1,6},{l,1,6}]

Table[Sort[NormalOrder[{1,{1,1,k,l}}]]==DropLeadingZeros/@Sort[{{1,{k+1,l+1}},{k,{k,l}}}/.{{0,0}->{0}}],{k,1,6},{l,1,6}]

Table[Sort[NormalOrder[{1,{k,l,1,1}}]]==DropLeadingZeros/@Sort[{{1,{k+1,l+1}},{l,{k,l}}}/.{{0,0}->{0}}],{k,1,6},{l,1,6}]