A program for finding the summation of an analytical expression

Question

Consider the Baker-Hausdorff formula for two operators $A$ and $B$: $$e^BA e^{-B} =A+[B,A]+\frac{1}{2!}[B,[B,A]]+\frac{1}{3!}[B,[B,[B,A]]]+....,$$ where $[A,B]=AB-BA$. In the case of my problem, $[B,a_1]=\alpha\, a_1+\beta \,a_2+\gamma,$ and $[B,a_2]=\alpha\, a_2+\beta \,a_1,$ where the parameters are either imaginary $(\alpha, \beta)$ or complex $(\gamma)$. $a_1$ and $a_2$ are operators. So, I need to find the following: $$ a_1\,+(\alpha a_1+\beta a_2+\gamma)+\frac{1}{2!}\Big\{\alpha\Big(\alpha a_1+\beta a_2+\gamma\Big)+\beta\Big(\alpha\, a_2+\beta \,a_1\Big)\Big\}+\frac{1}{3!}\Big\{\alpha\Big(\alpha (\alpha\, a_1+\beta \,a_2+\gamma)+\beta (\alpha\, a_2+\beta \,a_1)+\gamma\Big)+\beta\Big(\alpha\, (\alpha\, a_2+\beta \,a_1)+\beta \,(\alpha\, a_1+\beta \,a_2+\gamma)\Big)\Big\}+....$$

In other words, in each new term, one should replace $a_1$ by $(\alpha\, a_2+\beta \,a_1+\gamma)$ and $a_2$ with $(\alpha\, a_2+\beta \,a_1)$, and repeat this procedure $N$ times, and finally calculate the limit when $N\rightarrow \infty$. Since I'm not familiar with Mathematica programming, I tried to find the "pattern" of successive terms and thereby find the limit of n'th term as $N\rightarrow\infty$. I tried

f0 = a1;

f1 = Collect[α a1 + β a2 + γ /. {a1 -> α a1 + β a2 + γ, a2 -> α a2 + β a1}, {a1, a2}] // Simplify

f2 = Collect[% /. {a1 -> α a1 + β a2 + γ, a2 -> α a2 + β a1}, {a1, a2}] // Simplify

f3 = Collect[% /. {a1 -> α a1 + β a2 + γ, a2 -> α a2 + β a1}, {a1, a2}] // Simplify

f4 = Collect[% /. {a1 -> α a1 + β a2 + γ, a2 -> α a2 + β a1}, {a1, a2}] // Simplify

Collect[f0 + f1 + 1/2 f2 + 1/6 f3 + 1/4! f4, {a1, a2}]

But, unfortunately, the pattern is not so clear to me. So I think it would be better to find it systematically. The answer is probably either exponential or trigonometric.

In fact, I'm trying to find equation 6 of the following reference: arxiv.org/pdf/1609.00075.pdf

Answer 1

Here is a systematic way of generating the nested commutators:

comm[b, a_Times] := comm[b, Expand@a]
comm[b, a_Plus] := comm[b, #] & /@ a
comm[b, x_ a[y_]] := x comm[b, a[y]]
comm[b, x_ /; FreeQ[x, a]] := 0
comm[b, a[1]] := α a[1] + β a[2] + γ
comm[b, a[2]] := α a[2] + β a[1]

Then, for instance, one can do

terms = Table[Nest[comm[b, #] &, a[1], n], {n, 0, 3}]
(* {a[1], γ + α a[1] + β a[2],
     β (β a[1] + α a[2]) + α (γ + α a[1] + β a[2]),
     2 α β (β a[1] + α a[2]) + α^2 (γ + α a[1] + β a[2]) + β^2 (γ + α a[1] + β a[2])} *)

Finally, we make a transformation of the operators to symmetric and anti-symmetric combinations using

First@Solve[{a[p] == (a[1] + a[2])/Sqrt[2], a[m] == (a[1] - a[2])/Sqrt[2]}, {a[1], a[2]}]
(* {a[1] -> a[m]/Sqrt[2] + a[p]/Sqrt[2], a[2] -> 1/2 (-Sqrt[2] a[m] + Sqrt[2] a[p])} *)

Applying this transformation yields the following:

terms = Table[Nest[comm[b, #] &, a[1], n], {n, 0, 3}];
terms /. First@Solve[{a[p] == (a[1] + a[2])/Sqrt[2], a[m] == (a[1] - a[2])/Sqrt[2]}, {a[1], a[2]}];
Collect[%, _a, Factor]
(* {a[m]/Sqrt[2] + a[p]/Sqrt[2],
    γ + ((α - β) a[m])/Sqrt[2] + ((α + β) a[p])/Sqrt[2],
    α γ + ((α - β)^2 a[m])/Sqrt[2] + ((α + β)^2 a[p])/Sqrt[2],
    (α^2 + β^2) γ + ((α - β)^3 a[m])/Sqrt[2] + ((α + β)^3 a[p])/Sqrt[2]} *)

From there, I think it is apparent how some of the exponentials in the linked paper arise.