This notebook develops a method for computing the normal ordering of expressions involving non-commuting variables, specifically, powers of a formal variable and its formal derivative , where . Using NonCommutativeMultiply, we preserve non-commutative structure and systematically apply distributivity rules. A key transformation rewrites expressions into a standard form , referred to as normal order. We implement a recursive tensor-based algorithm to efficiently compute the normal-ordered form of powers of such expressions. The correctness of the procedure is verified by comparing results against direct evaluation using a differential operator algebra. This approach enables symbolic manipulation and evaluation of complex operator expressions with guaranteed preservation of non-commutativity. Everything we’ve built for the pair carries over verbatim to the canonical pair of annihilation and creation operators since . With Mathematica V14.3 to be launched soon, I briefly showed the magic of new non-commutative functions for the particular case discussed here, and how easily they can be done only in a few lines of code; classic Mathematica magic!
x
D=d/dx
[D,x]=1
f(x,D)=
∑
n,m
c
nm
n
x
m
D
(x,D)
(a,)
†
a
[a,]=1
†
a
For formulas, we used MaTex . More info in how to use it can be found here: https://www.wolframcloud.com/obj/nikm/Published/TikZ.nb
Normal ordering of a function
Normal ordering of a function
Given a function with =d/dx we want to write it in the normal order
fx,
⋀
D
⋀
D
fx,=
⋀
D
N
∑
n=0
M
∑
m=0
c
n,m
n
x
m
⋀
D
We shall denote derivative by the formal D, and likewise x by the formal x, so they would never be assigned any value.
Write an expression in terms of x and its derivative:
In[]:=
exp=HoldFormD.x.+D.^2x.+1+x.^2+D.^3x.^2+D.+D.-D.-x.^3^2+x.+x.^2D.^3^2
c.
1
c.
2
c.
3
c.
5
c.
6
c.
7
c.
8
c.
9
Out[]=
c.
1
c.
2
2
D.
c.
3
2
(1+x.)
c.
5
c.
6
3
D.
2
x.
2
D.-
3
D.-x.
c.
7
c.
8
2
x.+
c.
9
2
x.
3
D.
Note that derivative of x and x do not commute, , so to incorporate this feature, one should use NonCommutativeMultiply instead of the normal multiplication Times
,x=1
⋀
D
In[]:=
timesToNC=exp//.{Power[x_,y_]:>NonCommutativeMultiply@@ConstantArray[x,y],Times[x_,y_]:>NonCommutativeMultiply[x,y]}//ReleaseHold
Out[]=
D.+D.+(-1)****D.+(-1)****x.****D.+(-1)****x.****D.+(-1)****x.**D.+(-1)****D.+(-1)****x.****D.+(-1)****x.****D.+(-1)****x.+x.+**x.**x.**D.**D.**D.**x.+**x.**x.**D.**D.**D.+**D.**x.+**(1+**x.)**(1+**x.)+**D.**D.**x.+**D.**D.**D.**x.**x.
c.
7
c.
8
c.
7
c.
8
c.
7
c.
8
c.
7
c.
8
c.
7
c.
8
c.
7
c.
8
c.
9
c.
9
c.
1
c.
3
c.
5
c.
5
c.
2
c.
6
Additionally, one must incorporate distributive features while preserving non-commutativity:
In[]:=
ClearAll[expandNCM]expandNCM[(h:NonCommutativeMultiply)[a___,b_Plus,c___]]:=Distribute[h[a,b,c],Plus,h,Plus,expandNCM[h[##]]&]expandNCM[(h:NonCommutativeMultiply)[a___,b_Times,c___]]:=Most[b]expandNCM[h[a,Last[b],c]]expandNCM[a_]:=ExpandAll[a]
Return a new expression that applies distributivity within the constraints of non-commutative multiplication:
In[]:=
newExp=timesToNC/.x_NonCommutativeMultiply:>expandNCM[x]/.NonCommutativeMultiply[x__]/;AnyTrue{x},FreeQ#,x.D.&:>Module{num,compl},num=Select{x},FreeQx.D.;compl=DeleteElements[{x},num];If[compl=={},Times@@num,(Times@@num)NonCommutativeMultiply@@compl]/.NonCommutativeMultiply[x__]:>NonCommutativeMultiply@@(Times@@@Split[{x}])/.NonCommutativeMultiply[x_]:>x
Out[]=
2
x.
2
D.
c.
1
2
D.
c.
2
c.
3
c.
3
c.
5
2
x.
c.
3
2
c.
5
3
D.
2
x.
c.
6
4
D.
3
c.
7
6
D.
6
c.
7
3
D.
2
c.
7
c.
8
3
D.
2
c.
7
c.
8
2
D.
2
c.
7
c.
8
2
D.
2
c.
7
c.
8
5
D.
5
c.
7
c.
8
5
D.
5
c.
7
c.
8
4
D.
5
c.
7
c.
8
2
D.
3
D.
5
c.
7
c.
8
3
D.
2
D.
5
c.
7
c.
8
4
D.
5
c.
7
c.
8
2
x.
2
D.
c.
7
2
c.
8
2
D.
2
x.
c.
7
2
c.
8
2
x.
c.
7
2
c.
8
c.
7
2
c.
8
c.
7
2
c.
8
2
x.
4
D.
4
c.
7
2
c.
8
4
D.
2
x.
4
c.
7
2
c.
8
4
D.
4
c.
7
2
c.
8
2
x.
3
D.
4
c.
7
2
c.
8
2
D.
2
x.
2
D.
4
c.
7
2
c.
8
3
D.
2
x.
4
c.
7
2
c.
8
3
D.
4
c.
7
2
c.
8
2
D.
2
D.
4
c.
7
2
c.
8
3
D.
4
c.
7
2
c.
8
3
D.
4
c.
7
2
c.
8
2
D.
2
D.
4
c.
7
2
c.
8
3
D.
4
c.
7
2
c.
8
2
D.
4
c.
7
2
c.
8
2
D.
4
c.
7
2
c.
8
2
D.
4
c.
7
2
c.
8
3
x.
3
c.
8
3
x.
3
c.
8
3
x.
3
D.
3
c.
7
3
c.
8
3
D.
3
x.
3
c.
7
3
c.
8
3
D.
2
x.
3
c.
7
3
c.
8
2
x.
3
D.
3
c.
7
3
c.
8
3
x.
2
D.
3
c.
7
3
c.
8
2
D.
3
x.
3
c.
7
3
c.
8
2
x.
2
D.
3
c.
7
3
c.
8
2
D.
2
x.
3
c.
7
3
c.
8
2
x.
2
D.
3
c.
7
3
c.
8
2
x.
2
D.
3
c.
7
3
c.
8
2
D.
2
x.
3
c.
7
3
c.
8
2
x.
2
D.
3
c.
7
3
c.
8
2
D.
2
x.
3
c.
7
3
c.
8
2
D.
2
x.
3
c.
7
3
c.
8
2
D.
3
c.
7
3
c.
8
2
D.
3
c.
7
3
c.
8
2
x.
3
c.
7
3
c.
8
2
x.
3
c.
7
3
c.
8
3
c.
7
3
c.
8
3
c.
7
3
c.
8
4
x.
2
D.
2
c.
7
4
c.
8
2
D.
4
x.
2
c.
7
4
c.
8
2
D.
3
x.
2
c.
7
4
c.
8
2
x.
2
D.
2
x.
2
c.
7
4
c.
8
3
x.
2
D.
2
c.
7
4
c.
8
4
x.
2
c.
7
4
c.
8
3
x.
2
c.
7
4
c.
8
2
x.
2
x.
2
c.
7
4
c.
8
3
x.
2
c.
7
4
c.
8
3
x.
2
c.
7
4
c.
8
2
x.
2
x.
2
c.
7
4
c.
8
3
x.
2
c.
7
4
c.
8
2
x.
2
c.
7
4
c.
8
2
x.
2
c.
7
4
c.
8
2
x.
2
c.
7
4
c.
8
5
x.
c.
7
5
c.
8
5
x.
c.
7
5
c.
8
4
x.
c.
7
5
c.
8
2
x.
3
x.
c.
7
5
c.
8
3
x.
2
x.
c.
7
5
c.
8
4
x.
c.
7
5
c.
8
6
x.
6
c.
8
3
x.
3
D.
c.
9
2
x.
3
D.
c.
9
2
x.
3
D.
2
x.
3
D.
2
c.
9
Given an expression as NonCommutativeMultiply[…], write it as ; or in other words, find {{, , ..., }, {, , ..., }}
p
∏
j=1
n
j
x
m
j
d
m
j
dx
n
1
n
2
n
p
m
1
m
2
m
p
In[]:=
ClearAll[nmGen];nmGen[pat_List]:=Module{pat1,pat2,coeff},pat1=IfMatchQpat[[1]],x.x.^_,pat,Prepend[pat,0];pat2=IfMatchQpat1[[-1]],D.D.^_,pat1,Append[pat1,0];coeff=Partition[pat2,2]/.x.^n_:>n,x.:>1,D.^n_:>n,D.:>1;Reverse/@Transpose[coeff]
Given an expression as , write it in the normal order (a,b) with (a,b)=a-+i,b-+i
where is the factorial power. Note that (0,0)=1 while all other elements (a,b) are zero.
p
∏
j=1
n
j
x
m
j
d
m
j
dx
∑
j
n
j
∑
a=0
∑
j
m
j
∑
b=0
ℋ
p
a
x
b
d
b
dx
ℋ
j+1
m
j+1
∑
i=0
ℋ
j
n
j+1
m
j+1
m j+1 |
i |
(i)
a-+1
n
j+1
(i)
α
ℋ
0
H
0
In[]:=
ClearAll[normalOrdering]normalOrdering[{n_,m_}]/;Length[m]==Length[n]:=Module{p=Length[m],=Total[n],ℳ=Total[m],ℋ},ℋ=Fold[Function[{ℋ,j},Table[Sum[If[0<=a-n[[j+1]]+i<=Total[n[[;;j]]]&&0<=b-(m[[j+1]]-i)<=Total[m[[;;j]]],ℋ[[a-(n[[j+1]]-i)+1,b-(m[[j+1]]-i)+1]]Binomial[m[[j+1]],i]FactorialPower[a-n[[j+1]]+i,i],0],{i,0,m[[j+1]]}],{a,0,Total[n[[;;j+1]]]},{b,0,Total[m[[;;j+1]]]}]],SparseArray[{1,1}->1,{+1,ℳ+1}],Range[0,p-1]];Sumℋ[[a+1,b+1]]*If[a==0||b==0,Times,NonCommutativeMultiply]x.^a,D.^b,{a,0,},{b,0,ℳ}//Expand
Given the expression with NonCommutativeMultiply, write it in the normal order:
In[]:=
expNO=newExp/.NonCommutativeMultiply[pat__]:>normalOrdering[nmGen[{pat}]]//Expand
Out[]=
2
x.
2
D.
c.
1
c.
1
c.
2
2
D.
c.
2
c.
3
c.
3
c.
5
2
x.
c.
3
2
c.
5
c.
6
2
D.
c.
6
2
x.
3
D.
c.
6
4
D.
3
c.
7
6
D.
6
c.
7
2
D.
2
c.
7
c.
8
3
D.
2
c.
7
c.
8
4
D.
5
c.
7
c.
8
5
D.
5
c.
7
c.
8
c.
7
2
c.
8
c.
7
2
c.
8
2
x.
2
D.
c.
7
2
c.
8
2
D.
4
c.
7
2
c.
8
3
D.
4
c.
7
2
c.
8
2
x.
4
D.
4
c.
7
2
c.
8
2
x.
3
c.
8
3
x.
3
c.
8
3
c.
7
3
c.
8
3
c.
7
3
c.
8
2
x.
2
D.
3
c.
7
3
c.
8
3
x.
3
D.
3
c.
7
3
c.
8
2
x.
2
c.
7
4
c.
8
3
x.
2
c.
7
4
c.
8
4
x.
2
D.
2
c.
7
4
c.
8
4
x.
c.
7
5
c.
8
5
x.
c.
7
5
c.
8
6
x.
6
c.
8
2
x.
2
D.
c.
9
3
x.
3
D.
c.
9
2
x.
4
D.
2
c.
9
3
x.
5
D.
2
c.
9
4
x.
6
D.
2
c.
9
Now let’s run some testing. We will use the following function which is given in the documentation page of NonCommutativeMultiply to represent composition in an algebra of differential operators.
In[]:=
DOperator[A_,expr_]/;FreeQ[A,D]:=AexprDOperator[L1_+L2_,expr_]:=DOperator[L1,expr]+DOperator[L2,expr]DOperator[A_L_,expr_]/;FreeQ[A,D]:=ADOperator[L,expr]DOperator[A:HoldPattern[D[__]&],expr_]:=A[expr]DOperator[L1__**L2_,expr_]:=Expand[DOperator[L1,DOperator[L2,expr]]]DOperator[L1_^n_,expr_]:=Nest[DOperator[L1,#1]&,expr,n]
First, evaluate the final result of the expression using NonCommutativeMultiply before applying distributivity:
In[]:=
DOperatornewExp/.D.:>#1&,f[x.]//Expand
∂
x.
Out[]=
2
x.
c.
1
c.
3
c.
3
c.
5
2
x.
c.
3
2
c.
5
c.
7
2
c.
8
2
x.
3
c.
8
3
c.
7
3
c.
8
2
x.
2
c.
7
4
c.
8
4
x.
c.
7
5
c.
8
6
x.
6
c.
8
′
f
c.
1
′
f
c.
2
′
f
c.
6
′
f
c.
7
2
c.
8
′
f
3
x.
3
c.
8
′
f
3
c.
7
3
c.
8
′
f
3
x.
2
c.
7
4
c.
8
′
f
5
x.
c.
7
5
c.
8
′
f
′′
f
c.
2
′′
f
c.
6
′′
f
2
c.
7
c.
8
′′
f
2
x.
c.
7
2
c.
8
′′
f
4
c.
7
2
c.
8
′′
f
2
x.
3
c.
7
3
c.
8
′′
f
4
x.
2
c.
7
4
c.
8
′′
f
2
x.
c.
9
′′
f
2
x.
c.
6
(3)
f
2
c.
7
c.
8
(3)
f
4
c.
7
2
c.
8
(3)
f
3
x.
3
c.
7
3
c.
8
(3)
f
3
x.
c.
9
(3)
f
3
c.
7
(4)
f
5
c.
7
c.
8
(4)
f
2
x.
4
c.
7
2
c.
8
(4)
f
2
x.
2
c.
9
(4)
f
5
c.
7
c.
8
(5)
f
3
x.
2
c.
9
(5)
f
6
c.
7
(6)
f
4
x.
2
c.
9
(6)
f
Now let’s get the outcome from the normal ordered expression:
In[]:=
DOperatorexpNO/.D.:>#1&,f[x.]
∂
x.
Out[]=
2
x.
c.
1
c.
3
c.
3
c.
5
2
x.
c.
3
2
c.
5
c.
7
2
c.
8
2
x.
3
c.
8
3
c.
7
3
c.
8
2
x.
2
c.
7
4
c.
8
4
x.
c.
7
5
c.
8
6
x.
6
c.
8
′
f
c.
1
′
f
c.
2
′
f
c.
6
′
f
c.
7
2
c.
8
′
f
3
x.
3
c.
8
′
f
3
c.
7
3
c.
8
′
f
3
x.
2
c.
7
4
c.
8
′
f
5
x.
c.
7
5
c.
8
′
f
′′
f
c.
2
′′
f
c.
6
′′
f
2
c.
7
c.
8
′′
f
2
x.
c.
7
2
c.
8
′′
f
4
c.
7
2
c.
8
′′
f
2
x.
3
c.
7
3
c.
8
′′
f
4
x.
2
c.
7
4
c.
8
′′
f
2
x.
c.
9
′′
f
2
x.
c.
6
(3)
f
2
c.
7
c.
8
(3)
f
4
c.
7
2
c.
8
(3)
f
3
x.
3
c.
7
3
c.
8
(3)
f
3
x.
c.
9
(3)
f
3
c.
7
(4)
f
5
c.
7
c.
8
(4)
f
2
x.
4
c.
7
2
c.
8
(4)
f
2
x.
2
c.
9
(4)
f
5
c.
7
c.
8
(5)
f
3
x.
2
c.
9
(5)
f
6
c.
7
(6)
f
4
x.
2
c.
9
(6)
f
Compare previous results:
In[]:=
%==%%
Out[]=
True
Normal ordering as a single function
Normal ordering as a single function
For the sake of efficiency, I will implement a slightly different code in this section.
In[]:=
ClearAll[expandNCM]expandNCM[(h:NonCommutativeMultiply)[a___,b_Plus,c___]]:=Distribute[h[a,b,c],Plus,h,Plus,expandNCM[h[##]]&]expandNCM[(h:NonCommutativeMultiply)[a___,b_Times,c___]]:=Most[b]expandNCM[h[a,Last[b],c]]expandNCM[a_]:=ExpandAll[a]
In[]:=
ClearAll[nmGen];nmGen[pat_List]:=Module{pat1,pat2,coeff},pat1=IfMatchQpat[[1]],x.x.^_,pat,Prepend[pat,0];pat2=IfMatchQpat1[[-1]],D.D.^_,pat1,Append[pat1,0];coeff=Partition[pat2,2]/.x.^n_:>n,x.:>1,D.^n_:>n,D.:>1;Transpose[coeff]
Create a normal-ordered expression:
Check the normal-ordered expression is the same as the original one:
With a function’s power, the above procedure for the normal ordering works, but won’t be efficient. In the next section, we will elaborate on an efficient way to handle those cases.
Normal ordering of a function’s power
Normal ordering of a function’s power
Generate parameters in above equation:
Generate an array with symbolic elements:
Note that the above expression is already normal-ordered
Generate the transition matrix:
Check its dimension:
Apply the transition matrix repeatedly:
Find the result by applying the function repeatedly (this may take some times if large values are assigned for parameters):
Compare the results:
Normal ordering as a single function
Normal ordering as a single function
Check the result is the same as before:
Magic of version 14.3
Magic of version 14.3
Mathematica 14.3 is launching soon, featuring several exciting updates, among them, new tools for non-commutative algebra. These enhancements allow you to replicate the above results in just a few lines of code; classic Wolfram Language magic. Note below functions are available only in V14.3.
Define variables:
Define scalars:
Define the non-commutative algebra:
Define algebraic relations:
Given algebra and the expression we used in the previous section, find the reduction of the polynomial in normal-ordered terms:
Given algebra and the expression we used in the previous section, find the reduction of the polynomial in normal-ordered terms for the power of that expression:
To compare the results, rewrite our result in the previous section in a similar format as non-commutative algebra:
Show the results are the same:
Annihilation and creation operators: same algebra
Annihilation and creation operators: same algebra
As noted in the abstract, the same algebra applies to annihilation and creation operators, with examples below, future work will refine this beyond basic replacement rules.
Generate an expression:
Find the normal-ordered form of the expression raised to power 2:
Generate an expression:
Find the normal-ordered form of the expression raised to power 2:
Generate an expression:
Find the normal-ordered form of the expression raised to power 3:
Generate an expression:
Find the normal-ordered form of the expression raised to power 2:
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Symbolic normal ordering of non-commutative differential operator expressions
by Mohammad Bahrami
Wolfram Community, STAFF PICKS, June 12, 2025
https://community.wolfram.com/groups/-/m/t/3478087
by Mohammad Bahrami
Wolfram Community, STAFF PICKS, June 12, 2025
https://community.wolfram.com/groups/-/m/t/3478087