This computational essay develops a systematic implementation of the Magnus expansion for solving time-dependent operator differential equations with non-commuting Hamiltonians at different times. The notebook introduces efficient algorithms to generate closed-form expressions for the nested-commutator structure of the Magnus series up to arbitrary order, with flexible options for commutator forms, symbolic time integration, and Gröbner basis–based non-commutative polynomial reduction. Applications are presented for physically relevant Hamiltonians, including (i) bosonic creation/annihilation operators, (ii) Pauli spin operators under arbitrary time-dependent drives, and (iii) driven two-level systems. The calculations explicitly demonstrate simplifications in higher-order Magnus terms and provide symbolic and numerical verification against exact time-ordered evolution. This framework provides a versatile computational tool for analyzing quantum dynamics, Floquet theory, quantum control, and effective Hamiltonian engineering, enabling exploration of driven quantum systems where non-commutativity plays a crucial role.
We want to solve
Y(t)
t
=A(t)Y(t)
with
Y(0)=I
, with generally non-commutating
A
operator,
[A(
t
1
),A(
t
2
)]≠0
. The exact solution can be written as
Y(t)=
Ω(t)
e
. The Magnus expansion expresses
Ω(t)
as an infinite series of time-ordered integrals of nested commutators:
Ω(t)=
∞
∑
n=1
Ω
n
(t)
, where
Ω
n
(t)
is homogeneous of “order
n
” in
A
. It contains
n
time integrals and
n−1
commutators. For the case of quantum unitary evolution, i.e.
U(t)
t
=-H(t)U(t)
, one needs to replace
A(t)
with
-H(t)
.
An efficient closed form (right-nested, independent commutators) for
Ω
n
is given as follows:
Ω
n
(t)=
∑
σ∈

n-1
d(σ)
(-1)
n
n-1
d(σ)

t
∫
0

t
1
t
1
∫
0

t
2
…
t
n-1
∫
0

t
n
A
t
σ(1)
,A
t
σ(2)
,…,A
t
σ(n-1)
,A
t
σ(n)
…
​where

n-1
is permutations of
{1,2,…,n-1}
, and
d(σ)
is the number of descents of σ (indices
i
with
σ(i)>σ(i+1)
), and
c
n
=
d(σ)
(-1)
n
n-1
d(σ)

.
Define two functions, one for calculating
d(σ)
and the other one to calculate the coefficient
c
n
=
d(σ)
(-1)
n
n-1
d(σ)

:
In[]:=
ClearAll[descentCount,coeff];​​descentCount[σ_List]:=Count[Differences[σ],_?Negative];​​coeff[n_,σ_List]:=With{d=descentCount[σ]},
(-1)^d
n Binomial[n-1,d]

Ω
n
(t)=
t
∫
0

t
1
t
1
∫
0

t
2
…
t
n-1
∫
0

t
n
Ω
n,com
(t)​​
Ω
n,com
(t)=
1
n
∑
σ∈

n-1
d(σ)
(-1)

n-1
d(σ)

A
t
σ(1)
,A
t
σ(2)
,…,A
t
σ(n-1)
,A
t
σ(n)
…
​​MagnusTerms the order-
n
of Magnus term. The first input argument is n and the 2nd one is the corresponding operator. Given n, MagnusTerms automatically assign time variables. If a third argument is given, then it can be Dot, Composition, TensorProduct or NonCommutativeMultiply.​MagnusTerms takes three options: - “CommutatorForm”: if True, it hold the form of commutators; if False, it will do the operation-“TimeIntegration”: if True, it returns
Ω
n
(t)
in above formula; if False, it returns only
Ω
n,com
(t)
​- “Reduce”: if not give, it does nothing, but if a specific function is given (as NonCommutativePolynomialReduce)
In[]:=
ClearAll[MagnusTerms]​​Options[MagnusTerms]={"CommutatorForm"->True,"TimeIntegration"->False,"Reduce"->None};​​MagnusTerms[n_,op_,OptionsPattern[]]:=Module{σList=Permutations[Range[n-1]],ΩCommutator},​​(*
Ω
n,com
(t)=
1
n
∑
σ∈

n-1
d(σ)
(-1)

n-1
d(σ)

A
t
σ(1)
,A
t
σ(2)
,…,A
t
σ(n-1)
,A
t
σ(n)
…*)​​ΩCommutator=Sum[coeff[n,j]Fold[If[OptionValue["CommutatorForm"],HoldForm[Commutator[#2,#1]]&,Commutator[#2,#1]&],Reverse@Append[Table[op[
t
i
],{i,j}],op[
t
n
]]],{j,σList}];​​IfOptionValue["TimeIntegration"],​​(*
Ω
n
(t)=
t
∫
0

t
1
t
1
∫
0

t
2
…
t
n-1
∫
0

t
n
Ω
n,com
(t)*)​​Integrate[If[MatchQ[OptionValue["Reduce"],None],ΩCommutator,OptionValue["Reduce"][ΩCommutator]],Sequence@@Flatten/@Transpose[{Table[
t
i
,{i,n}],Prepend[Table[{0,
t
i
},{i,n-1}],{0,t}]}]],If[MatchQ[OptionValue["Reduce"],None],ΩCommutator,OptionValue["Reduce"][ΩCommutator]]​​​​MagnusTerms[n_,op_,alg_,OptionsPattern[]]:=Module{σList=Permutations[Range[n-1]],ΩCommutator},​​(*
Ω
n,com
(t)=
1
n
∑
σ∈

n-1
d(σ)
(-1)

n-1
d(σ)

A
t
σ(1)
,A
t
σ(2)
,…,A
t
σ(n-1)
,A
t
σ(n)
…*)​​ΩCommutator=Sum[coeff[n,j]Fold[If[OptionValue["CommutatorForm"],HoldForm[Commutator[#2,#1,alg]]&,Commutator[#2,#1,alg]&],Reverse@Append[Table[op[
t
i
],{i,j}],op[
t
n
]]],{j,σList}];​​(*
Ω
n
(t)=
t
∫
0

t
1
t
1
∫
0

t
2
…
t
n-1
∫
0

t
n
Ω
n,com
(t)*)​​If[OptionValue["TimeIntegration"],Integrate[If[MatchQ[OptionValue["Reduce"],None],ΩCommutator,OptionValue["Reduce"][ΩCommutator]],Sequence@@Flatten/@Transpose[{Table[
t
i
,{i,n}],Prepend[Table[{0,
t
i
},{i,n-1}],{0,t}]}]],If[MatchQ[OptionValue["Reduce"],None],ΩCommutator,OptionValue["Reduce"][ΩCommutator]]]​​
Generate
Ω
1,com
(t)
which will be
A(t)
In[]:=
MagnusTerms1,A.
Out[]=
A.[
t
1
]
Generate
Ω
1
(t)
In[]:=
MagnusTerms1,A.,"TimeIntegration"->True
Out[]=
t
∫
0
A.[
t
1
]
t
1
Generate
Ω
2,com
(t)
In[]:=
MagnusTerms2,A.//TraditionalForm
Out[]//TraditionalForm=
1
2
[A.(
t
1
),A.(
t
2
)]
Check it is the same as
1
2
[A(
t
1
),A(
t
2
)]
In[]:=
MagnusTerms2,A.,"CommutatorForm"->False-
1
2
CommutatorA.[
t
1
],A.[
t
2
]//NonCommutativeExpand
Out[]=
0
Do the commutation operation (default):
In[]:=
MagnusTerms2,A.,"CommutatorForm"->False
Out[]=
1
2
(A.[
t
1
]**A.[
t
2
]-A.[
t
2
]**A.[
t
1
])
Do the commutation operation by Dot:
In[]:=
MagnusTerms2,A.,Dot,"CommutatorForm"->False
Out[]=
1
2
(A.[
t
1
].A.[
t
2
]-A.[
t
2
].A.[
t
1
])
Generate
Ω
2
(t)
In[]:=
MagnusTerms2,A.,"TimeIntegration"->True//TraditionalForm
Out[]//TraditionalForm=
1
2
t
∫
0
t
1
∫
0
[A.(
t
1
),A.(
t
2
)]
t
2

t
1
Generate
Ω
3,com
(t)
In[]:=
MagnusTerms3,A.//TraditionalForm
Out[]//TraditionalForm=
1
3
[A.(
t
1
),[A.(
t
2
),A.(
t
3
)]]-
1
6
[A.(
t
2
),[A.(
t
1
),A.(
t
3
)]]
Check it is the same as
1
3
[A(
t
1
),[A(
t
2
),A(
t
3
)]]-
1
6
[A(
t
2
),[A(
t
1
),A(
t
3
)]]
In[]:=
MagnusTerms3,A.,"CommutatorForm"->False-
1
3
CommutatorA.[
t
1
],CommutatorA.[
t
2
],A.[
t
3
]-
1
6
CommutatorA.[
t
2
],CommutatorA.[
t
1
],A.[
t
3
]//NonCommutativeExpand
Out[]=
0
Generate
Ω
3
(t)
In[]:=
MagnusTerms3,A.,"TimeIntegration"->True//TraditionalForm
Out[]//TraditionalForm=
t
∫
0
t
1
∫
0
t
2
∫
0
1
3
[A.(
t
1
),[A.(
t
2
),A.(
t
3
)]]-
1
6
[A.(
t
2
),[A.(
t
1
),A.(
t
3
)]]
t
3

t
2

t
1
Generate
Ω
4,com
(t)
In[]:=
MagnusTerms4,A.//TraditionalForm
Out[]//TraditionalForm=
1
4
[A.(
t
1
),[A.(
t
2
),[A.(
t
3
),A.(
t
4
)]]]-
1
12
[A.(
t
1
),[A.(
t
3
),[A.(
t
2
),A.(
t
4
)]]]-
1
12
[A.(
t
2
),[A.(
t
1
),[A.(
t
3
),A.(
t
4
)]]]-
1
12
[A.(
t
2
),[A.(
t
3
),[A.(
t
1
),A.(
t
4
)]]]-
1
12
[A.(
t
3
),[A.(
t
1
),[A.(
t
2
),A.(
t
4
)]]]+
1
12
[A.(
t
3
),[A.(
t
2
),[A.(
t
1
),A.(
t
4
)]]]
Check it is the same as which will be
1
4
[A(
t
1
),[A(
t
2
),[A(
t
3
),A(
t
4
)]]]-
1
12
[A(
t
1
),[A(
t
3
),[A(
t
2
),A(
t
4
)]]]-
1
12
[A(
t
2
),[A(
t
1
),[A(
t
3
),A(
t
4
)]]]-
1
12
[A(
t
2
),[A(
t
3
),[A(
t
1
),A(
t
4
)]]]-
1
12
[A(
t
3
),[A(
t
1
),[A(
t
2
),A(
t
4
)]]]+
1
12
[A(
t
3
),[A(
t
2
),[A(
t
1
),A(
t
4
)]]]

Some applications

For the following applications, I will consider cases that the operators are not assigned explicitly but only their algebra (non-commutative one) is define.
The Hamiltonian above is expressed in the interaction picture; see Sec. 4.1.1 of the preprint (arXiv:0810.5488) for details.
Define the variables within the non-commutative algebra together with their corresponding commutation relations:
For scalar variables in above definitions, we use the same format as in the Magnus series defined in the previous section.
A function that simplifies and reduces the non-commutative polynomial of the annihilation and creation operators:
We use Nest with ReleaseHold since the Magnus terms defined in the previous section are expressed in the hold commutator form.
As noted earlier, we chose a Hamiltonian that does not commute at different times. Show the Hamiltonian at different times do not commute:
It is interesting that the third them and so on are all zero. Therefore, only first and 2nd terms contribute to the exponential:
Let’s check the commutator relationship for Pauli operators (we will use this to define our non-commutative algebra later)
In the codes above, we assigned explicit values to the Pauli matrices. Later, we will use only their algebraic relations, treating the operators symbolically without explicit assignment.
Define the variables within the non-commutative algebra together with their corresponding commutation relations:
A function that simplifies and reduces the non-commutative polynomial of Pauli matrices:
As noted earlier, we chose a Hamiltonian that does not commute at different times. Show the Hamiltonian at different times do not commute:
In above code, we have applied NonCommutativePolynomialReduce two times to take care of all reductions (note the difference in the order of variables).
Define corresponding algebra and Groebner Basis:
Define a function to simplify corresponding polynomials:
For each term, find the coefficient for Pauli matrices (keep in mind that we will have only linear terms):
Now let’s repeat the same process, but this time assigning explicit values to Pauli matrics and doing the commutation/Dot calculations
As expected, the coefficients in front of Pauli matrices are the same the ones we obtained from non - commutative algebraic calculations
Some numerical exploration
Calculate the unitary operator:

References

https://arxiv.org/pdf/0810.5488
https://arxiv.org/pdf/1710.10851

CITE THIS NOTEBOOK

Efficient computation of Magnus series for non-commuting operators in quantum dynamics​
by Mohammad Bahrami​
Wolfram Community, STAFF PICKS, August 21, 2025
​https://community.wolfram.com/groups/-/m/t/3534231