Normal-ordering bosonic operators

F.P. Laussy
Russian Quantum Center
http://laussy.org/wiki/BosonNormalOrder

v°0.11: 21 May 2022
v°0.2: Revised 31 May 2022
v°1.0: Bug fixed and wikilaussy page created: 13 June 2023

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 product 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]). One can make more complex sequences in this way, 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.:

BosonNormalOrder[{{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

BosonNormalOrderCorrelator[{1,1,0}]

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

BosonNormalOrderCorrelator[{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]];
remainseq=DropLeadingZeros/@remainseq;
doneseq=Join[doneseq,Select[remainseq,Length[#[[2]]]<=2&]];
remainseq=Select[remainseq,Length[#[[2]]]>2&];
(* 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]]]}]}
]

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

BosonNormalOrder[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):

BosonNormalOrder[{
{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:

BosonNormalOrderCorrelator[{1,1,0}]

BosonNormalOrder[{{1,{1,1,0}}}]

This finds that a (a^(†5))=5 (a^(†4))+(a^(†5))a

BosonNormalOrderCorrelator[{1,5,0}]

The Stirling numbers:

BosonNormalOrderCorrelator[{1,1,1,1,1,1,1,1,0}]

Checking analytical formulas

We check up to some max numbers; importantly, this should work when the exponents are zero (nmin=0), i.e., the operators are  not actually present:

nmin=0;nmax=10;

When we check the right-hand-side (explicit result) one must be careful to make sure that, e.g., {0,0} is contracted as {0} (this is done by DropLeadingZeros/@) and that results with negative weights are removed (this is done by Select[]):

The greyish areas are the initial expression (top line) and its explicit formula (bottom line):

This checks that a (a^(†k))(a^l) (SuperDagger[a])=(a^(† k+1))(a^(l+1))+(k+l+1) (a^(† k))(a^l)+k l (a^(† k-1))(a^(l-1))

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,k,l,1,0}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+1,l+1}},{k+l+1,{k,l}},{k l,{k-1,l-1}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (SuperDagger[a])a (a^(†k))(a^l)(SuperDagger[a])a=(a^(† k+2))(a^(l+2))+(k+l+1) (a^(† k+1))(a^(l+1))+k l (a^(† k))(a^l)
(same as previous multiplied left and right by normally-ordered results)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,1,k,l,1,1}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+2,l+2}},{k+l+1,{k+1,l+1}},{k l,{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (SuperDagger[a]) a (SuperDagger[a]) a (a^(† k)) (a^l)=(a^(†[k+2])) (a^(l+2))+(1+2 k) (a^(†[k+1])) (a^(l+1))+(k^2) (a^(† k)) (a^l)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,1,1,1,k,l}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+2,l+2}},{1+2k,{k+1,l+1}},{((k )^2),{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that a (SuperDagger[a]) a (a^(† k))=(a^(†[k+1])) (a^2)+(1+2 k) (a^(† k)) a+(k^2) (a^(† k-1))

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,1,1,k,0}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+1,2}},{1+2k,{k,1}},{((k )^2),{k-1,0}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (a^(† k)) (a^l) ((((SuperDagger[a]) a))^2)=(a^(†[k+2])) (a^(l+2))+(2 l+1) (a^(†[k+1])) (a^(l+1))+(l^2) (a^(† k)) (a^l)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{k,l,1,1,1,1}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+2,l+2}},{1+2l,{k+1,l+1}},{(l^2),{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (a^l) (SuperDagger[a]) a (SuperDagger[a])=(a^(†2)) (a^(l+1))+(2 l+1) (SuperDagger[a]) (a^l)+(l^2) (a^(l-1))

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{l,1,1,1,0}]
]==Sort[Select[DropLeadingZeros/@
{{1,{2,l+1}},{1+2l,{1,l}},{(l^2),{l-1}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (SuperDagger[a]) a (a^(† k)) (a^l) (SuperDagger[a]) a=(a^(†[k+2])) (a^(l+2))+(k+l+1) (a^(†[k+1])) (a^(l+1))+k l (a^(† k)) (a^l)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,1,k,l,1,1}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+2,l+2}},{1+k+l,{k+1,l+1}},{k l,{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that a (a^(† k)) (a^l) (SuperDagger[a])==(a^(†[k+1])) (a^(l+1))+(k+l+1) (a^(† k)) (a^l)+k l (a^(†[k-1])) (a^(l-1))

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,k,l,1,0}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+1,l+1}},{1+k+l,{k,l}},{k l,{k-1,l-1}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that a (a^(†[k+1])) (a^l)=(a^(†[k+1])) (a^(l+1))+(1+k) (a^(† k)) (a^l)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,k+1,l}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+1,l+1}},{1+k,{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (a^(† k)) (a^(l+1)) (SuperDagger[a])=(a^(†[k+1])) (a^(l+1))+(1+l) (a^(† k)) (a^l)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{k,l+1,1,0}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+1,l+1}},{1+l,{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (SuperDagger[a]) a (a^(† k)) (a^l)=(a^(†[k+1])) (a^(l+1))+k (a^(† k)) (a^l)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{1,1,k,l}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+1,l+1}},{k,{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

This checks that (a^(† k)) (a^l) (SuperDagger[a]) a=(a^(†[k+1])) (a^(l+1))+l (a^(† k)) (a^l)

Union@Flatten[Table[Sort[
BosonNormalOrderCorrelator[{k,l,1,1}]
]==Sort[Select[DropLeadingZeros/@
{{1,{k+1,l+1}},{l,{k,l}}}
,#[[1]]!=0&]],{k,nmin,nmax},{l,nmin,nmax}]]

ToExpression["a^{\\dagger k}a^{l}a^\\dagger a=a^{\\dagger(k+1)}a^{l+1}+la^{\\dagger k}a^l",TeXForm,HoldForm]