Expanding Exponentials in Noncommutative Algebras: 
Zassenhaus, Baker–Campbell–Hausdorff and Magnus Formulas

Mads Bahrami, Wolfram Research
quantum@wolfram.com

#WolframTechConf

ABSTRACT

This talk presents an in-depth exploration of the Zassenhaus, Baker–Campbell–Hausdorff and Magnus formulas, which decompose exponentials of operator sums in noncommutative settings. We detail their derivations, applications and recent algorithmic formulations for high-degree expansions. Special attention is given to symbolic implementations in Wolfram Language, showcasing efficient computation using custom-built functions and algebraic structures. Applications to operator theory, Lie algebras and symbolic matrix analysis will be discussed.

Baker-Campbell-Hausdorff (BCH) series:

(e^(Subscript[𝒜, 1]))(e^(Subscript[𝒜, 2]))((…e)^(Subscript[𝒜, m]))=(e^((∑_(j=1)^m)(Subscript[𝒜, j])+(∑_(n=2)^∞)(Subscript[𝒵, n])))

Zassenhaus series:

(e^((∑_(j=1)^m)(Subscript[𝒜, m])))=(∏_(j=1)^m)(e^(Subscript[𝒜, j]))(∏_(k=2))(e^(Subscript[𝒲, k]))

Magnus series:

((dY(t))/(dt))=A(t)Y(t) ⇒Y(t)=(e^((∑_(n=1)^∞)(Subscript[Ω, n])(t)))

Initialization

ClearAll[BakerCampbellHausdorffTerms,ZassenhausTerms,MagnusTerms]
BakerCampbellHausdorffTerms=ResourceFunction["https://www.wolframcloud.com/obj/mohammadb/DeployedResources/Function/BakerCampbellHausdorffTerms/"];
ZassenhausTerms=ResourceFunction["ZassenhausTerms"];
MagnusTerms=ResourceFunction["MagnusTerms"];

Adjoint function

ClearAll[adjointAction];
adjointAction[a_,b_,n_:1]:=adjointAction[a,b,n,NonCommutativeMultiply]
adjointAction[a_,b_,0,alg:(_NonCommutativeAlgebra|Dot|{Dot,_}|Composition|TensorProduct|NonCommutativeMultiply)]:=b
adjointAction[a_,b_,n_:1,alg:(_NonCommutativeAlgebra|Dot|{Dot,_}|Composition|TensorProduct|NonCommutativeMultiply)]:=
	Nest[If[MatchQ[alg,NonCommutativeMultiply],Commutator[a,#],Commutator[a,#,alg]]&,b,n]/;n>=1

Define a function to reduce/simplify a NCA polynomial of a and (SuperDagger[a]):

ClearAll[NCPolynomialReductionBosonic]
NCPolynomialReductionBosonic[exp_,scalars_:{}]:=NonCommutativePolynomialReduce[exp,{Commutator[a,(SuperDagger[a])]-1},{a,(SuperDagger[a])},NonCommutativeAlgebra[<|"ScalarVariables"->scalars|>]][[-1]]

Given the algebra and variables, define a function to reduce a polynomial in (Subscript[J, j]) (Cartesian basis):

ClearAll[NCPolynomialReductionSU2]
NCPolynomialReductionSU2[exp_,vars_,scalars_:{},scale_:1]/;Length[vars]==3:=NonCommutativePolynomialReduce[exp,DeleteCases[0]@Flatten@(Outer[Commutator,vars,vars]-I scale TensorContract[TensorProduct[LeviCivitaTensor[3],vars],{{3,4}}]),{vars},NonCommutativeAlgebra[<|"ScalarVariables"->scalars|>]][[-1]]

With (Subscript[J, ±])=(Subscript[J, 1])±i (Subscript[J, 2]), one finds [(Subscript[J, 3]),(Subscript[J, ±])]=±(Subscript[J, ±]) and [(SubPlus[J]),(SubMinus[J])]=2(Subscript[J, 3])
So it means {(SubPlus[J]),(SubMinus[J]),(Subscript[J, 3])} also have SU(2) algebra. They are usually referred to as Cartan–Weyl basis.

ClearAll[NCPolynomialReductionSU2Weyl]
NCPolynomialReductionSU2Weyl[exp_,vars_,scalars_:{}]/;Length[vars]==3:=NonCommutativePolynomialReduce[exp,{Commutator@@vars[[;;2]]-2vars[[-1]],Commutator[vars[[-1]],vars[[1]]]-vars[[1]],Commutator[vars[[-1]],vars[[2]]]+vars[[2]]},{vars},NonCommutativeAlgebra[<|"ScalarVariables"->scalars|>]][[-1]]

Adjoint action

The adjoint action of A on another operator B is define as the commutation [A,B]. A direct interesting result is the exponential adjoint action (e^A)B (e^(-A))=(e^(Subscript[ad, A]))B with (Subscript[ad, A])[B]=[A,B]

adjointAction[A,B]

The j-th iterated adjoint of A is defined as the n-fold nested commutator of A with B, denoted by (Subscript[ad, A]^((j)))[B]= [A, [A,…[A,B]…]]:

adjointAction[A,B,2]==Commutator[A,Commutator[A,B]]

adjointAction[A,B,3]==Commutator[A,Commutator[A,Commutator[A,B]]]

We shall study some examples from quantum optics. To do so, let us first define a function to reduce noncommutative polynomials of annihilation and creation operators.

Examples

Consider (e^(λ (SuperDagger[a])a))a (e^(-λ (SuperDagger[a])a)). Compute 10 terms in the exponential adjoint:

NCPolynomialReductionBosonic[Sum[NCPolynomialReductionBosonic[adjointAction[λ (SuperDagger[a])**a,a,j],{λ}]/j!,{j,0,10}],{λ}]

Verify the summation of 10 terms is the same as Taylor expansion of (e^(-λ))

Normal[Series[(E^(-λ)),{λ,0,10}]]==CoefficientRules[%,{a,(SuperDagger[a])}][[1,-1]]

So we managed to verify (e^(λ (SuperDagger[a])a))a (e^(-λ (SuperDagger[a])a))=(e^(-λ))a

Consider (e^((1/2)((SuperStar[ζ])(a^2)-ζ ((SuperDagger[a])^2))))a (e^(-(1/2)((SuperStar[ζ])(a^2)-ζ ((SuperDagger[a])^2)))) with ζ=r (e^(iϕ)). Compute 10 terms in the exponential adjoint:

NCPolynomialReductionBosonic[Sum[NCPolynomialReductionBosonic[adjointAction[((𝓇 (E^(-I ϕ)) a**a-𝓇 (E^(I ϕ)) (SuperDagger[a])**(SuperDagger[a]))/2),a,j]/j!,{𝓇,ϕ}],{j,0,10}],{𝓇,ϕ}]

Verify (e^((1/2)((SuperStar[ζ])(a^2)-ζ ((SuperDagger[a])^2))))a (e^(-(1/2)((SuperStar[ζ])(a^2)-ζ ((SuperDagger[a])^2))))=u a+v (SuperDagger[a]) with u=cosh(r) and v=-(e^(i ϕ))sinh(r):

Normal[Series[{Cosh[𝓇],(E^(I ϕ))Sinh[𝓇]},{𝓇,0,10}]]==CoefficientRules[%,{a,(SuperDagger[a])}][[All,-1]]

Baker-Campbell-Hausdorff (BCH) series

The BCH series expresses the product of two or more exponentials of non-commuting operators as a single exponential:  (e^(Subscript[𝒜, 1]))(e^(Subscript[𝒜, 2]))((…e)^(Subscript[𝒜, m]))=(e^((∑_(j=1)^m)(Subscript[𝒜, j])+(∑_(n=2)^∞)(Subscript[𝒵, n])))

For two operators, (e^(X+Y)), one gets the Dynkin’s formula: (Subscript[𝒵, m])=(∑_(n=1)^m)((((-1))^(n-1))/n)(∑_((Subscript[p, 1])+(Subscript[q, 1])>0,…,(Subscript[p, n])+(Subscript[q, n])>0,(∑_(j=1)^n)((Subscript[p, j])+(Subscript[q, j]))=m-1))(((Subscript[ad, X]^(Subscript[p, 1]))(Subscript[ad, Y]^(Subscript[q, 1]))…(Subscript[ad, X]^(Subscript[p, n]))(Subscript[ad, Y]^(Subscript[q, n])))/((Subscript[p, 1])!(Subscript[q, 1])!… (Subscript[p, n])!(Subscript[q, n])!))(X+Y).

Show a few terms of (e^A)(e^B)=(e^(A+B+(∑_(n=2))(Subscript[𝒵, n])))

ClearAll[A,B];
TableForm[Table[{n,TraditionalForm@BakerCampbellHausdorffTerms[{ A,B},n,"CommutatorForm"->True]},{n,2,4}],TableHeadings->{None, {"n","(Subscript[𝒵, n])"}},TableAlignments->Left]

Most quantum operators belong to small Lie algebras (Heisenberg–Weyl, SU(2), SU(1,1)), where the commutator chains truncate or close, allowing the BCH formula to yield clean, closed-form exponential identities.

Normal ordering and Weyl relations

We will explore some examples involving photonic operators, where the single-mode photon field satisfies the Bosonic commutation relation [(OverHat[a]),(SuperDagger[(OverHat[a])])]=1.

Annihilation and creation operator:

Consider the expression: (e^(μ (SuperDagger[(OverHat[a])])))(e^(ν (OverHat[a])))=(e^(μ (SuperDagger[(OverHat[a])])+ν (OverHat[a])+(∑_n)(Subscript[𝒵, n]))) 
Compute the order-2 term in BCH series:

BakerCampbellHausdorffTerms[{μ a,ν (SuperDagger[a])},2]

Reduce/simplify it:

NCPolynomialReductionBosonic[%,{μ,ν}]

Note the result is a c-number only.

Compute the order-3 term of BCH series:

BakerCampbellHausdorffTerms[{μ a,ν (SuperDagger[a])},3]

Reduce it:

NCPolynomialReductionBosonic[%,{μ,ν}]

Likewise, the other terms disappear as well. Thus one finds: (e^(μ (SuperDagger[(OverHat[a])])))(e^(ν (OverHat[a])))=(e^(μ (SuperDagger[(OverHat[a])])+ν (OverHat[a])+μ ν/2))

Canonical position and momentum

Consider canonical position and momentum operators: 
(OverHat[q])=(1/(Sqrt[2]))((OverHat[a])+(SuperDagger[(OverHat[a])])) and (OverHat[p])=(1/(i(Sqrt[2])))((OverHat[a])-(SuperDagger[(OverHat[a])]))

Define canonical position and momentum operators:

q=(a+(SuperDagger[a]))/(Sqrt[2]);p=(a-(SuperDagger[a]))/(I(Sqrt[2]));

Verify [(OverHat[q]),(OverHat[p])]=i

NCPolynomialReductionBosonic@Commutator[q,p]

Now let’s try to compute (e^(i α(OverHat[q]) ))(e^(i β(OverHat[p]) )).

Compute the second term in the BCH series:

NCPolynomialReductionBosonic[BakerCampbellHausdorffTerms[{I α q,I β p},2],{α,β}]

Show that the higher-order terms are zero:

Table[NCPolynomialReductionBosonic[BakerCampbellHausdorffTerms[{I α q,I β p},n],{α,β}],{n,3,6}]

Thus it means (e^(i α(OverHat[q]) ))(e^(i β(OverHat[p]) ))=(e^(-i α β/2))(e^(i (α(OverHat[q])+ β(OverHat[p)]) ))

Displacement operator

Displacement operator is defined as 𝒟(α)=(e^(α (SuperDagger[a])-(SuperStar[α])a)). Show that 𝒟(α)𝒟(β)=(e^((α (SuperStar[β])-(SuperStar[α])β)/2))𝒟(α+β)

NCPolynomialReductionBosonic[BakerCampbellHausdorffTerms[{α (SuperDagger[a])-Conjugate[α]a,β (SuperDagger[a])-Conjugate[β]a},2],{α,β}]

Note that the result is a c-number.

Show that the higher-order terms are zero:

Table[NCPolynomialReductionBosonic[BakerCampbellHausdorffTerms[{α (SuperDagger[a])-Conjugate[α]a,β (SuperDagger[a])-Conjugate[β]a},n],{α,β}],{n,3,6}]

Compute the second through fifth terms in the BCH series of 𝒟((Subscript[α, 1]))𝒟((Subscript[α, 2]))𝒟((Subscript[α, 3])) and verify 𝒟((Subscript[α, 1]))𝒟((Subscript[α, 2]))𝒟((Subscript[α, 3]))=(e^((∑_j)((Subscript[α, j])(SuperDagger[a])-(SuperStar[(Subscript[α, j])])a)+(1/2)(∑_(i<j))((Subscript[α, i])(SuperStar[(Subscript[α, j])])-(Subscript[α, j])(SuperStar[(Subscript[α, i])]))))

Table[NCPolynomialReductionBosonic[BakerCampbellHausdorffTerms[Table[(Subscript[α, j]) (SuperDagger[a])-Conjugate[(Subscript[α, j])]a,{j,3}],n],Table[(Subscript[α, j]) ,{j,3}]],{n,2,5}]

SU(2) spin rotations

Set up the variables and their SU(2) algebra [(Subscript[J, i]),(Subscript[J, j])]=(Subscript[iϵ, ijk])(Subscript[J, k]) with ϵ the 3-dimensional Levi-Civita totally antisymmetric tensor. Note that in our relations we do not include quadratic identities such as ((Subscript[J, j])^2)=𝕀, because these are representation-dependent (e.g., for Pauli matrices) rather than algebraic identities.

Log[(e^(-i δt (Subscript[θ, 1])(Subscript[J, 1])))(e^(-i δt (Subscript[θ, 2])(Subscript[J, 2])))]=-i δt (Subscript[H, eff])
Our goal to find effective Hamiltonian describing the above rotations.

We will use rotation angles (Subscript[θ, j]) and also time variable δt. We should set them as scalar variable in the algebra.

Compute the second term of the BCH series:

vars=Table[(Subscript[J, j]),{j,3}];
scalars=Append[Table[(Subscript[θ, j]),{j,3}],δt];

Compute (Subscript[𝒵, 2]) for (e^(-i δt (Subscript[θ, 1])(Subscript[J, 1])))(e^(-i δt (Subscript[θ, 2])(Subscript[J, 2]))):

NCPolynomialReductionSU2[#,vars,scalars]&@BakerCampbellHausdorffTerms[{-I δt (Subscript[θ, 1]) (Subscript[J, 1]),-I δt (Subscript[θ, 2])(Subscript[J, 2])},2]

Compute (Subscript[H, eff]) up to the second order of δt:

Sum[NCPolynomialReductionSU2[BakerCampbellHausdorffTerms[{-I δt (Subscript[θ, 1]) (Subscript[J, 1]),-I δt (Subscript[θ, 2])(Subscript[J, 2])},n],vars,scalars],{n,3}]/(-I δt)//Expand

Numerical example

Note that in the previous examples, all operators were symbolic. One can use the same thing with numerical operators, too.

(e^x)(e^y)=(e^A) with x=(({{α, 0}, {0, -α}})), y=(({{0, β}, {0, 0}}))
and the result is A=x+((2α)/(1-(e^(-2α))))y=(({{α, β(((2α)/(1-(e^(-2α)))))}, {0, -α}}))

Let’s calculate 10 terms for A:

x=α PauliMatrix[3];y={{0,β},{0,0}};
A=Sum[BakerCampbellHausdorffTerms[{x,y},n,Dot],{n,11}];
A//MatrixForm

Focus on (Subscript[A, 12]) element. One can show that the Taylor series of the function below follows the same result as above.

Normal[β Series[((2α)/(1-Exp[-2α])),{α,0,10}]]==A[[1,2]]

Zassenhaus series

The Zassenhaus formula is defined as (e^((∑_(j=1)^m)(Subscript[𝒜, m])))=(∏_(j=1)^m)(e^(Subscript[𝒜, j]))(∏_(k=2))(e^(Subscript[𝒲, k])) with (Subscript[𝒲, n]) the order-n term that involves all possible Lie monomials of degree n in terms of nested commutators of the operators. Zassenhaus ensures every factor (Subscript[𝒲, n]) you generate remains inside the same Lie algebra (closed under commutators). In other words, if (Subscript[𝒜, m])∈𝕘 then all (Subscript[𝒲, n])∈𝕘. Thus the whole infinite product has exponents that stay inside the same Lie algebra.

Show Zassenhaus terms for (e^(A+B))=(e^A)(e^B)(∏_(k=2))(e^(Subscript[𝒲, k])):

ClearAll[A,B];
TableForm[Table[{n,TraditionalForm@ZassenhausTerms[{ A,B},n,"CommutatorForm"->True]},{n,2,4}],TableHeadings->{None, {"n","(Subscript[𝒲, n])"}},TableAlignments->Left]

Displacement operator:

Consider the displacement operator 𝒟(α)=(e^(α (SuperDagger[a])-(SuperStar[α])a)).

Compute the reduced form of the order-2 term in Zassenhaus series:

NCPolynomialReductionBosonic[ZassenhausTerms[{α (SuperDagger[a]),-Conjugate[α]a},2],{α}]

Note the result is a c-number, which implies higher order terms are zero:

Table[NCPolynomialReductionBosonic[ZassenhausTerms[{α (SuperDagger[a]),-Conjugate[α]a},n],{α}],{n,3,6}]

Therefore, one gets: 𝒟(α)=(e^(-(1/2)|α((|)^2)))(e^(α (SuperDagger[a])))(e^(-(SuperStar[α])a))

Now let’s switch the first and second operator and get a new result as 𝒟(α)=(e^((1/2)|α((|)^2)))(e^(-(SuperStar[α])a))(e^(α (SuperDagger[a]))):

Compute the second term in Zassenhaus series:

NCPolynomialReductionBosonic[ZassenhausTerms[{-Conjugate[α]a,α (SuperDagger[a])},2],{α}]

Compute a few terms (3rd and more) and show they are zero:

Table[NCPolynomialReductionBosonic[ZassenhausTerms[{-Conjugate[α]a,α (SuperDagger[a])},n],{α}],{n,3,6}]

SU(2) Trotter error

H=(1/2)((Subscript[θ, 1])(Subscript[J, 1])+(Subscript[θ, 2])(Subscript[J, 2])+(Subscript[θ, 3])(Subscript[J, 3]))
(e^(-i δt H))=(e^(-i δt ((Subscript[θ, 1])/2)(Subscript[J, 1])))(e^(-i δt ((Subscript[θ, 2])/2)(Subscript[J, 2])))(e^(-i δt ((Subscript[θ, 1])/2)(Subscript[J, 3])))(e^( ((δt)^2) (Subscript[𝒜, 2])))(e^( ((δt)^3) (Subscript[𝒜, 3])))…

Compute the second term (order of ((δt)^2))

NCPolynomialReductionSU2[ZassenhausTerms[{-I δt (Subscript[θ, 1])/2 (Subscript[J, 1]),-I δt (Subscript[θ, 2])/2 (Subscript[J, 2]),-I δt (Subscript[θ, 3])/2 (Subscript[J, 3])},2],vars,scalars]

Compute the third term (order of ((δt)^3))

NCPolynomialReductionSU2[ZassenhausTerms[{-I δt (Subscript[θ, 1])/2 (Subscript[J, 1]),-I δt (Subscript[θ, 2])/2 (Subscript[J, 2]),-I δt (Subscript[θ, 3])/2 (Subscript[J, 3])},3],vars,scalars]

SU(2) Gauss/Wei–Norman decomposition

With (Subscript[J, ±])=(Subscript[J, 1])±i (Subscript[J, 2]), one finds [(Subscript[J, 3]),(Subscript[J, ±])]=±(Subscript[J, ±]) and [(SubPlus[J]),(SubMinus[J])]=2(Subscript[J, 3])
So it means {(SubPlus[J]),(SubMinus[J]),(Subscript[J, 3])} also have SU(2) algebra. They are usually referred to as Cartan–Weyl basis.

U(θ)=(e^(-i θ (Subscript[J, 1])))=(e^(-((i θ)/2) ((SubPlus[J])+(SubMinus[J]))))=(e^(-((i θ)/2) (SubPlus[J])))(e^(-((i θ)/2) (SubMinus[J])))(∏_(k=2))(e^(Subscript[𝒲, k]))

The series above is an infinite right-ordered product. But if 𝕘 is finite-dimensional (which is the case here, dim(SU(2))=3), one can re-express that product as a finite product of exponentials of a chosen basis of 𝕘 with scalar (c-number) coefficients. This is the content of the Wei–Norman theorem on Lie groups (Gauss/Wei–Norman factorization):

U(θ)=(e^(α(θ) (SubPlus[J])))(e^(β(θ) (Subscript[J, 3])))(e^(γ(θ) (SubMinus[J])))

Additionally for (U^(-1)), we get:

(U^(-1))=(e^(-γ (SubMinus[J])))(e^(-β (Subscript[J, 3])))(e^(-α (SubPlus[J])))

Putting them together, one finds:

((dU)/(dθ))(U^(-1))=-i (Subscript[J, 1])=(d/(dθ))((e^(α(θ) (SubPlus[J])))(e^(β(θ) (Subscript[J, 3])))(e^(γ(θ)) (SubMinus[J]))))(e^(-γ(θ) (SubMinus[J])))(e^(-β(θ) (Subscript[J, 3])))(e^(-α(θ) (SubPlus[J])))
=α'[θ](SubPlus[J])+β'[θ]((e^(α(θ) (SubPlus[J])))(Subscript[J, 3])(e^(-α(θ) (SubPlus[J]))))+γ'[θ]((e^(α(θ) (SubPlus[J])))(e^(β(θ) (Subscript[J, 3])))(SubMinus[J])(e^(-β(θ) (Subscript[J, 3])))(e^(-α(θ) (SubPlus[J]))))
=α'[θ](SubPlus[J])+β'[θ](e^(Subscript[Ad, α(θ) (SubPlus[J])]))(Subscript[J, 3])+γ'[θ](e^(Subscript[Ad, β(θ) (Subscript[J, 3])]))(e^(Subscript[Ad, α(θ) (SubPlus[J])]))(SubMinus[J])

Define noncommutative variables:

vars2={(SubPlus[J]),(SubMinus[J]),(Subscript[J, 3])};

Calculate (e^(Subscript[Ad, α(θ) (SubPlus[J])]))(Subscript[J, 3]) up to (10^th) term and simplify it:

NCPolynomialReductionSU2Weyl[Sum[NCPolynomialReductionSU2Weyl[(1/(n!))adjointAction[α[θ] (SubPlus[J]),(Subscript[J, 3]),n],vars2,{α[θ]}],{n,0,10}],vars2,{α[θ]}]

As one can see, the 3rd and higher terms are zero, therefore one finds: (e^(Subscript[Ad, α(θ) (SubPlus[J])]))(Subscript[J, 3])=(Subscript[J, 3])-γ[θ](SubMinus[J]).

Calculate (e^(α[θ](SubPlus[J]) ))(e^(β[θ] (Subscript[J, 3])))(SubMinus[J])(e^(-β[θ] (Subscript[J, 3])))(e^(-α[θ](SubPlus[J]) ))=(e^(Subscript[Ad, β(θ) (Subscript[J, 3])]))(e^(Subscript[Ad, α(θ) (SubPlus[J])]))(SubMinus[J]) up to (10^th) term and simplify it:

Sum[NCPolynomialReductionSU2Weyl[(1/(m!))adjointAction[α[θ](SubPlus[J]),Sum[NCPolynomialReductionSU2Weyl[(1/(n!))adjointAction[β[θ] (Subscript[J, 3]),(SubMinus[J]),n],vars2,{β[θ]}],{n,0,10}],m],vars2,{α[θ],β[θ]}],{m,0,10}]

Given the series expansion, one can verify that above equation is the same as (SubMinus[J]) (e^(-β[θ]))+2 (Subscript[J, 3]) α[θ] (e^(-β[θ]))-(SubPlus[J]) ((α[θ])^2) (e^(-β[θ])).

Putting above results together and simplify it, one finds:

finalResult=NCPolynomialReductionSU2Weyl[γ'[θ]((SubMinus[J]) (E^(-β[θ]))+2 (Subscript[J, 3]) α[θ] (E^(-β[θ]))-(SubPlus[J]) ((α[θ])^2) (E^(-β[θ])))+β'[θ]((Subscript[J, 3])-(SubPlus[J]) α[θ])+(SubPlus[J]) α'[θ],vars2,{α[θ],β[θ],γ[θ],γ'[θ],β'[θ],α'[θ]}]

((dU)/(dθ))(U^(-1))=-i (Subscript[J, 1]), expand it in terms of noncommutative variables:

eq1=Lookup[Table[UnitVector[3,j],{j,3}]][CoefficientRules[-I(((SubPlus[J])+ (SubMinus[J]))/2),vars2]]/._Missing:>0

So the same for the final result we obtained above:

eq2=Lookup[Table[UnitVector[3,j],{j,3}]][CoefficientRules[finalResult,vars2]]/._Missing:>0

Equate coefficients:

des=Thread[ eq1== eq2]

Solve them:

sol=DSolveValue[des,{α[θ],β[θ],γ[θ]},θ]//FullSimplify

Set an initial value for (Subscript[c, 1]):

Rest[sol/.{θ->0,(C[1])->0}]=={0,0}//Reduce

Plug the result back into the solution:

sol/.{(C[1])->0,(C[2])->2 Log[2],(C[3])->-1}//FullSimplify

Putting all results together, one finds U(θ)=(e^(-i θ (Subscript[J, 1])))=(e^(-i Tan[(θ/2)](SubPlus[J])))(e^(-Log[((Cos[(θ/2)])^2)](Subscript[J, 3])))(e^(-i Tan[(θ/2)](SubMinus[J])))

Magnus series

Magnus series is used for the solution of matrix linear initial-value problem ((dY(t))/(dt))=A(t)Y(t) with Y(0)=I, where in general [A((Subscript[t, 1])),A((Subscript[t, 2]))]≠0. The exact solution can be written as Y(t)=(e^(Ω(t))).

The Magnus expansion expresses Ω(t) as an infinite series of time-ordered integrals of nested commutators: Ω(t)=(∑_(n=1)^∞)(Subscript[Ω, n])(t), where (Subscript[Ω, 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. ((dU(t))/(dt))=-i H(t)U(t), one needs to replace A(t) with -i H(t).

An efficient closed form (right-nested, independent commutators) for (Subscript[Ω, n]) is given as follows: (Subscript[Ω, n])(t)=(∫_0^t)d(Subscript[t, 1])(∫_0^(Subscript[t, 1]))d(Subscript[t, 2])…(∫_0^(Subscript[t, n-1]))d(Subscript[t, n])(Subscript[Ω, n,com])(t) 
with (Subscript[Ω, n,com])(t)=(∑_(σ∈(Subscript[𝒢, n-1])))(Subscript[c, n]) [A((Subscript[t, σ(1)])),[A((Subscript[t, σ(2)])),…,[A((Subscript[t, σ(n-1)])),A((Subscript[t, σ(n)]))]…]] where (Subscript[𝒢, n-1]) is permutations of {1,2,…,n-1}, and (Subscript[c, n])=((((-1))^(d(σ)))/(n (({{n-1}, {d(σ)}})))) with d(σ) is the number of descents of σ (indices i with σ(i)>σ(i+1))

Show three terms of Magnus series (without time integration):

ClearAll[A];
TableForm[Table[{n,TraditionalForm@MagnusTerms[ A,n,"CommutatorForm"->True]},{n,2,4}],TableHeadings->{None, {"n","(Subscript[Ω, n])"}},TableAlignments->Left]

Examples

H(t)=f(t)((e^(i t))a+(e^(-i t))(SuperDagger[a]))

Define a time-dependent Hamiltonian as H(t)=f(t)((e^(i t))a+(e^(-i t))(SuperDagger[a])) :

ClearAll[h1,f]
h1[t_]:= f[t](a Exp[I t]+(SuperDagger[a]) Exp[-I t])

The Hamiltonian above is expressed in the interaction picture; see Sec. 4.1.1 of the preprint (arXiv:0810.5488) for details. To avoid any conflict, we shall use formal symbol of t for the time.

Find commutation term [H((Subscript[t, 1])),H((Subscript[t, 2]))] in the second term of Magnus expansion (Subscript[Ω, 2])(t)=(1/2) (∫_0^t)(∫_0^(Subscript[t, 1]))[H((Subscript[t, 1])),H((Subscript[t, 2]))]d(Subscript[t, 2])d(Subscript[t, 1]):

With[{n=2},MagnusTerms[h1,n,NCPolynomialReductionBosonic[#,Flatten[{#,f[#]}&/@Table[(Subscript[\[FormalT], i]),{i,n}]]]&]]//FullSimplify

Find commutation term, [H((Subscript[t, 1])),[H((Subscript[t, 2])),H((Subscript[t, 3]))]], in the third term of Magnus expansion (Subscript[Ω, 3])(t)= (∫_0^t)(∫_0^(Subscript[t, 1]))(∫_0^(Subscript[t, 2]))((1/3) [H((Subscript[t, 1])),[H((Subscript[t, 2])),H((Subscript[t, 3]))]]-(1/6) [H((Subscript[t, 2])),[H((Subscript[t, 1])),H((Subscript[t, 3]))]])d(Subscript[t, 3])d(Subscript[t, 2])d(Subscript[t, 1]):

With[{n=3},MagnusTerms[h1,n,NCPolynomialReductionBosonic[#,Flatten[{#,f[#]}&/@Table[(Subscript[\[FormalT], i]),{i,n}]]]&]]

Ω=(∑_n)(((-i))^n)(Subscript[Ω, n])=-i (Subscript[Ω, 1])-(Subscript[Ω, 2]):

Sum[(((-I))^n)MagnusTerms[h1,n,NCPolynomialReductionBosonic[#,{\[FormalCapitalT],(Subscript[\[FormalT], 1]),f[(Subscript[\[FormalT], 1])],(Subscript[\[FormalT], 2]),f[(Subscript[\[FormalT], 2])]}]&,"TimeIntegration"->True],{n,2}]//FullSimplify

It is interesting that the third them and so on are all zero. Therefore, only first and 2nd terms contribute to the exponential:
U(t)=Exp[-i((∫_0^t)d(Subscript[t, 1])f[(Subscript[t, 1])](e^(i (Subscript[t, 1]))))a-i((∫_0^t)d(Subscript[t, 1])f[(Subscript[t, 1])](e^(-i (Subscript[t, 1]))))(SuperDagger[a])-i (∫_0^t)(∫_0^(Subscript[t, 1]))d(Subscript[t, 2])d(Subscript[t, 1]) f[(Subscript[t, 1])] f[(Subscript[t, 2])] Sin[(Subscript[t, 1])-(Subscript[t, 2])]]] which can be rewritten as
U(t)=Exp[ -i β]Exp[ (α a-(SuperStar[α])(SuperDagger[a]))]
with α=-i(∫_0^t)d(Subscript[t, 1])f[(Subscript[t, 1])](e^(i (Subscript[t, 1]))) and β=(∫_0^t)(∫_0^(Subscript[t, 1]))d(Subscript[t, 2])d(Subscript[t, 1]) f[(Subscript[t, 1])] f[(Subscript[t, 2])] Sin[(Subscript[t, 1])-(Subscript[t, 2])]:

Table[MagnusTerms[h1,n,NCPolynomialReductionBosonic[#,Flatten[{#,f[#]}&/@Table[(Subscript[\[FormalT], i]),{i,n}]]]&],{n,3,8}]

H(t)=(Δ/2) (Subscript[σ, 3])+A Cos[ω t] (Subscript[σ, 1])

Define the Hamiltonian as H(t)=(Δ/2)(Subscript[σ, 3])+A cos[ω t](Subscript[σ, 1]):

ClearAll[h2,A,Δ,ω,σ]
h2[t_]:=Δ/2 (Subscript[σ, 3])+ A Cos[ω t] (Subscript[σ, 1])

Calculate {(Subscript[Ω, 1])(t),(Subscript[Ω, 2])(t),(Subscript[Ω, 3])(t),(Subscript[Ω, 4])(t)}:

Ω1To4=Table[MagnusTerms[h2,n,NCPolynomialReductionSU2[#,{(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])},Join[{Δ,A,ω},Table[(Subscript[\[FormalT], i]),{i,n}]],2]&,"TimeIntegration"->True],{n,4}];

For each term, find the coefficient for Pauli matrices (keep in mind that we will have only linear terms):

Ω1To4CoefficientRules=CoefficientRules[Total@Ω1To4,{(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])}];

Now let's repeat the same process, but this time assigning explicit values to Pauli matrics and doing the commutation/dot calculations:

ClearAll[h3,Δ,ω,A,σ]
h3[t_]:=Δ/2 PauliMatrix[3]+ A Cos[ω t] PauliMatrix[1]

Ω1To4Manual=(Subscript[a, 1])(Subscript[σ, 1])+(Subscript[a, 2])(Subscript[σ, 2])+(Subscript[a, 3])(Subscript[σ, 3])

Ω1To4Manual=Sum[MagnusTerms[h3,n,Dot,"TimeIntegration"->True],{n,4}];

Note (Subscript[a, i])=Tr[(Subscript[σ, i]).Ω1To4Manual]/2. One can calculate it and compare with the previous result:

Table[Tr[Ω1To4Manual.PauliMatrix[j]]/2,{j,3}]-Ω1To4CoefficientRules[[All,-1]]//FullSimplify

As expected, the coefficients in front of Pauli matrices are the same the ones we obtained from non - commutative algebraic calculations.

Some numerical exploration. To find approximation to unitary operator determined by H(t), calculate {-i (Subscript[Ω, 1])(t),(((-i))^2)(Subscript[Ω, 2])(t),(((-i))^3)(Subscript[Ω, 3])(t),(((-i))^4)(Subscript[Ω, 4])(t)}:

omega=(Table[(-I)^n,{n,4}]Ω1To4)/.{A->.2,Δ->.1,ω->1,(Subscript[σ, 1])->PauliMatrix[1],(Subscript[σ, 2])->PauliMatrix[2],(Subscript[σ, 3])->PauliMatrix[3]};

Calculate the unitary operator:

ham=(Δ/2 PauliMatrix[3]+ A Cos[ω t] PauliMatrix[1]/.{A->.2,Δ->.1,ω->1});
u=NDSolveValue[{\[FormalU]'[t]==-I ham.\[FormalU][t],\[FormalU][0]==IdentityMatrix[2]},\[FormalU],{t,0,tf=4π}]

Calculate {(Subscript[Ω, 1])(t),(∑_(n=1)^2)(Subscript[Ω, n])(t),(∑_(n=1)^3)(Subscript[Ω, n])(t),(∑_(n=1)^4)(Subscript[Ω, n])(t)}:

Ωn=Accumulate[omega];

Calculate (e^((∑_n)(Subscript[Ω, n])(t))) and compare it with the solution 𝒯((e^(-i(∫_0^t)H(s)ds))):

Legended[Show[Table[LogPlot[Norm[MatrixExp[Ωn[[j]]]-u[\[FormalCapitalT]]],{\[FormalCapitalT],10^-3,tf},PlotStyle->ColorData[97][j]],{j,4}],PlotRange->All,GridLines->Automatic,PlotLabel->"Norm difference of unitary operator from Magnus approximation",AspectRatio->1/2,Frame->True,FrameLabel->{"time", "||(E^((∑_n)(Subscript[Ω, n]) (t)))-𝒯(\!\(\*SuperscriptBox[\(E\), \(\(-I\)\\ \(\*SubsuperscriptBox[\(∫\), \(0\), \(t\)]H \((s)\) \[DifferentialD]s\)\)]\))||"}],Placed[LineLegend[Table[ColorData[97][j],{j,4}],Table["Magnus terms up to n="<>ToString[j],{j,4}]],{Right,Bottom}]]