# A unifying framework for perturbative exponential factorizations

Ana Arnal Email<EMAIL_ADDRESS>Fernando Casas Email<EMAIL_ADDRESS>Cristina Chiralt Email<EMAIL_ADDRESS>José Angel Oteo Email<EMAIL_ADDRESS>
###### Abstract

We propose a framework where Fer and Wilcox expansions for the solution of
differential equations are derived from two particular choices for the initial
transformation that seeds the product expansion. In this scheme intermediate
expansions can also be envisaged. Recurrence formulas are developed. A new
lower bound for the convergence of the Wilcox expansion is provided as well as
some applications of the results. In particular, two examples are worked out
up to high order of approximation to illustrate the behavior of the Wilcox
expansion.

Keywords: Sequences of linear transformations; Wilcox expansion; Fer
expansion; Zassenhaus formula; Bellman problem.

∗, ${\dagger}$,${\ddagger}$ Departament de Matemàtiques and IMAC, Universitat
Jaume I, 12071 Castellón, Spain,
$\S$ Departament de Física Teòrica, Universitat de València and Institute for
Integrative Systems Biology (I2SysBio), 46100-Burjassot, Valencia, Spain

##  1 Introduction

Linear differential equations of the form

$\dot{x}\equiv\frac{dx}{dt}=A(t)x,\qquad x(0)=x_{0}\in\mathbb{C}^{d},$ (1)

with $A(t)$ a $d\times d$ matrix whose entries are integrable functions of
$t$, are ever-present in many branches of science, the fundamental evolution
equation of Quantum Mechanics, the Schrödinger equation, being a particular
case. In consequence, solving Eq. (1) is of the greatest importance. In spite
of their apparent simplicity, however, they are seldom solvable in terms of
elementary functions, and so different procedures have been proposed along the
years to render approximate solutions. These are specially useful in the
analytical treatment of perturbative problems, such as those arising in the
time evolution of quantum systems [1], control theory or problems where time-
ordered products are involved [2]. Among them, exponential perturbative
expansions have received a great deal of attention, due to some remarkable
properties they possess. In particular, if Eq. (1) is defined in a Lie group,
the approximations they furnish also evolve in the same Lie group. As a
consequence, important qualitative properties of the exact solution are also
preserved by the approximations. Thus, if Eq. (1) represents the time-
dependent Schrödinger equation, then the approximate evolution operator is
still unitary, and as a consequence the total sum of (approximate) transition
probabilities is the unity, no matter where the expansion is truncated. There
are in fact many physical problems (in non-linear mechanics, optical
spectroscopy, magnetic resonance, etc.) involving periodic fast-oscillating
external fields that are also modeled by Eq. (1), with $A(t)$ periodic. In
that case, especially tailored expansions incorporating the well-known Floquet
theorem [3], such as the average Hamiltonian theory [4] and the Floquet–Magnus
expansion [5, 6] have also been proposed.

When dealing with the general problem (1), one of the most widely used
exponential approximations corresponds to the Magnus expansion [7],

$x(t)=\mathrm{e}^{\Omega(t)}x_{0},$ (2)

where $\Omega$ is an infinite series,

$\Omega(t)=\sum_{k=1}^{\infty}\Omega_{k}(t)$ (3)

whose terms are linear combinations of time-ordered integrals of nested
commutators of $A$ evaluated at different times (see [8] for a review,
including applications to several physical and mathematical problems). What is
more interesting for our purposes here is that this expansion can be related
with a coordinate transformation $x\longmapsto X$ rendering the original
system (1) into the trivial equation

$\frac{dX}{dt}=0,$ (4)

with the static solution $X(t)=X(0)=x_{0}$, and that the transformation is
given precisely by $x(t)=\exp(\Omega(t))X(t)$ [9].

In contrast to the Magnus expansion, the Floquet-Magnus expansion gets the
solution with two exponential transformations when $A(t)$ is periodic, whereas
other exponential perturbative expansions are based on infinite product
factorizations of $x(t)$,

$x(t)=\mathrm{e}^{\Omega_{1}(t)}\,\mathrm{e}^{\Omega_{2}(t)}\cdots\mathrm{e}^{\Omega_{n}(t)}\cdots\,x(0),$
(5)

such as those proposed by Fer and Wilcox. In fact, as pointed out in [10],
both expansions have a curious history worth to be sketched. It was Fer who
proposed the expansion that bears his name in [11], although he never applied
it to solve any specific problem. Bellman reviewed this paper in the
Mathematical Reviews111MR0104009., and even proposed the expansion as an
exercise in [12]. Nevertheless, Wilcox identified it in [1] with an
alternative factorization, Eq. (5), which was indeed a different and new type
of expansion. From them on, their historical trajectories move apart. Thus,
Fer expansion was rediscovered by Iserles [13] as a tool for the numerical
integration of linear differential equations and later on used in quantum
mechanics [14] and solid-state nuclear magnetic resonance [15], but also as a
Lie-group integrator [16, 17, 18]. On the other hand, Wilcox expansion has
been rediscovered several times in the literature, in particular in [19] in
the context of nonlinear control systems, and in [20] as a general tool for
approximating the time evolution operator in quantum mechanics.

The first goal in this work is to recast both infinite product expansions
within a unifying framework. This is done by considering, instead of just
_one_ exponential transformation as in the case of the Magnus expansion, a
_sequence_ of such transformations,
$\mathrm{e}^{\Omega_{1}(t)},\ldots,\mathrm{e}^{\Omega_{k}(t)},\ldots$, chosen
to satisfy certain requirements. To be more specific, suppose one replaces
$A(t)$ in Eq. (1) by $\lambda A(t)$, where $\lambda>0$ is a parameter. Then,
if the transformations are chosen so that each $\Omega_{k}(t)$ is proportional
to $\lambda^{k}$, we recover the Wilcox expansion, whereas we end up with the
Fer expansion when each $\Omega_{k}(t)$ is an infinite series in $\lambda$
whose first term is proportional to $\lambda^{2^{k-1}}$.

We also show that further alternative descriptions yield new factorizations.
This additional degree of freedom can be indeed used to deal better with the
features of the matrix $A(t)$ as in Floquet-Magnus, when $A(t)$ is periodic.

One might then consider this sequence of linear transformations as a
generalization of the concept of picture in Quantum Mechanics when Eq. (1)
refers to the Schrödinger equation.

Our second goal consists in obtaining, on the basis of this framework, new
results concerning Wilcox expansion. Thus, we develop a recursive procedure to
obtain every order of approximation in terms of nested commutators, as well as
a convergence radius bound. We also establish a formal connection of the
Wilcox expansion with the Zassenhaus formula [7, 21].

Eventually, the important problem of expanding the exponential
$\exp(A+\varepsilon B)$ for $\varepsilon>0$ and $A$, $B$, two generic non-
commuting operators will be addressed, and two applications of the results
given.

##  2 A sequence of transformations: the general case

Given the initial value problem (1), let us consider a linear change of
variables $x\longmapsto X_{1}$ of the form

$x(t)=\mathrm{e}^{\Omega_{1}(t)}X_{1}(t),\qquad\Omega_{1}(0)=0$ (6)

transforming the original system into

$\frac{dX_{1}}{dt}=B_{1}(t)X_{1}.$ (7)

For the time being, the generator $\Omega_{1}(t)$ of the transformation is not
specified. Then, $B_{1}(t)$ can be expressed in terms of $A(t)$ and
$\Omega_{1}(t)$ as follows. First, by inserting Eq. (6) into Eq. (1), and
taking Eq. (7) into account, one can obtain

$\frac{d}{dt}\exp(\Omega_{1})=A(t)\exp(\Omega_{1})-\exp(\Omega_{1})B_{1}(t),$
(8)

whence

$B_{1}(t)=\mathrm{e}^{-\Omega_{1}}A(t)\mathrm{e}^{\Omega_{1}}-\mathrm{e}^{-\Omega_{1}}\frac{d}{dt}\mathrm{e}^{\Omega_{1}}.$
(9)

The derivative of the matrix exponential can be written as [8]

$\frac{d}{dt}\exp(\Omega_{1}(t))=d\exp_{\Omega_{1}(t)}(\dot{\Omega}_{1}(t))\,\exp(\Omega_{1}),$
(10)

where the symbol $d\exp_{\Omega}(C)$ stands for the (everywhere convergent)
power series

$d\exp_{\Omega}(C)=\sum_{k=0}^{\infty}\frac{1}{(k+1)!}\mathrm{ad}_{\Omega}^{k}(C)\equiv\frac{\exp(\mathrm{ad}_{\Omega})-I}{\mathrm{ad}_{\Omega}}(C).$
(11)

Here $\mathrm{ad}_{\Omega}^{0}C=C$,
$\mathrm{ad}_{\Omega}^{k}C=[\Omega,\mathrm{ad}_{\Omega}^{k-1}C]$, and
$[\Omega,C]$ denotes the usual commutator. Therefore

$B_{1}(t)=\mathrm{e}^{-\Omega_{1}}\left(A(t)-d\exp_{\Omega_{1}(t)}(\dot{\Omega}_{1}(t))\right)\mathrm{e}^{\Omega_{1}}\equiv\mathrm{e}^{-\mathrm{ad}_{\Omega_{1}}}(B_{0}-G_{1})$
(12)

where

$B_{0}(t)\equiv A(t),\qquad G_{1}(t)\equiv
d\exp_{\Omega_{1}(t)}(\dot{\Omega}_{1}(t))$ (13)

and

$\mathrm{e}^{-\mathrm{ad}_{\Omega_{1}}}F=\sum_{k\geq
0}\frac{(-1)^{k}}{k!}\mathrm{ad}_{\Omega_{1}}^{k}F=\mathrm{e}^{-\Omega_{1}}F\,\mathrm{e}^{\Omega_{1}}.$
(14)

Of course, nothing prevents us from repeating the whole procedure above and
introduce a second transformation to Eq. (7) of the form

$X_{1}(t)=\mathrm{e}^{\Omega_{2}(t)}X_{2}(t),\qquad\Omega_{2}(0)=0,$ (15)

so that the new variables verify

$\frac{dX_{2}}{dt}=B_{2}(t)X_{2}.$ (16)

In general, for the $n$-th such linear transformation

$X_{n-1}(t)=\mathrm{e}^{\Omega_{n}(t)}X_{n}(t),\qquad\Omega_{n}(0)=0$ (17)

with

$\frac{dX_{n}}{dt}=B_{n}(t)X_{n},$ (18)

one has

$B_{n}(t)=\mathrm{e}^{-\mathrm{ad}_{\Omega_{n}}}(B_{n-1}-G_{n}),\qquad\mbox{
with }\qquad G_{n}(t)=d\exp_{\Omega_{n}(t)}(\dot{\Omega}_{n}(t))$ (19)

so that the solution of equation (1) is expressed as

$x(t)=\mathrm{e}^{\Omega_{1}(t)}\,\mathrm{e}^{\Omega_{2}(t)}\cdots\mathrm{e}^{\Omega_{n}(t)}X_{n}(t).$
(20)

Alternatively, we can write $B_{n}(t)$ in Eq. (19) as follows. Since it is
also true that [8]

$\displaystyle d\exp_{\Omega_{n}(t)}(\dot{\Omega}_{n}(t))$ $\displaystyle=$
$\displaystyle\left(\frac{d}{dt}\mathrm{e}^{\Omega_{n}}\right)\,\mathrm{e}^{-\Omega_{n}}=\left(\int_{0}^{1}\mathrm{e}^{x\Omega_{n}}\dot{\Omega}_{n}\mathrm{e}^{(1-x)\Omega_{n}}dx\right)\mathrm{e}^{-\Omega_{n}}$
(21) $\displaystyle=$
$\displaystyle\int_{0}^{1}\mathrm{e}^{x\,\Omega_{n}}\dot{\Omega}_{n}\mathrm{e}^{-x\,\Omega_{n}}dx,$

we have

$\displaystyle B_{n}=$
$\displaystyle\mathrm{e}^{-\Omega_{n}}B_{n-1}\mathrm{e}^{\Omega_{n}}-\int_{0}^{1}\mathrm{e}^{-u\,\Omega_{n}}\dot{\Omega}_{n}\mathrm{e}^{u\,\Omega_{n}}du$
(22) $\displaystyle=$ $\displaystyle\sum_{k\geq
0}\frac{(-1)^{k}}{(k+1)!}\left((k+1)\mathrm{ad}_{\Omega_{n}}^{k}B_{n-1}-\mathrm{ad}_{\Omega_{n}}^{k}\dot{\Omega}_{n}\right),\qquad
n\geq 1.$

The important point is of course how to choose $B_{n}$, or alternatively
$\Omega_{n}$, i.e., the specific requirements each transformation has to
satisfy in order to be useful to approximately solve Eq. (1). There are
obviously many possibilities and in the following we analyze two of them,
leading to two different and well known exponential perturbation
factorizations mentioned in the Introduction, namely the Wilcox [1] and Fer
[11] expansions.

##  3 Wilcox expansion

###  3.1 Recurrences

Let us introduce the (dummy) parameter $\lambda$ in Eq. (1) and replace $A$
with $\lambda A$. This is helpful when collecting coefficients, and at the end
we can always take $\lambda=1$.

Since the solution of Eq. (1) when $A$ is constant, or more generally when
$A(t_{1})A(t_{2})=A(t_{2})A(t_{1})$ for all $t_{1}\neq t_{2}$, is
$x(t)=\exp(\int_{0}^{t}A(u)du)$, it makes sense to take the generator for the
first transformation as

$\Omega_{1}(t)=\int_{0}^{t}B_{0}(u)\,du=\lambda\int_{0}^{t}A(u)\,du\equiv\lambda\,W_{1}(t).$
(23)

Then, according to Eqs. (12)-(13), we have

$C_{1}\equiv B_{0}-G_{1}=\sum_{k\geq 1}\lambda^{k}(b_{0,k}-g_{1,k})$ (24)

where

$\displaystyle b_{0,1}=A(t),\qquad\qquad b_{0,l}=0,\;\;l>1$ (25)
$\displaystyle
g_{1,1}=\dot{W}_{1},\qquad\qquad\;g_{1,l}=\frac{1}{l!}\mathrm{ad}_{W_{1}}^{l-1}\dot{W}_{1},\;\;l>1.$

With the choice $\dot{W}_{1}=A(t)$, it turns out that $B_{1}$ is a power
series in $\lambda$ starting with $\lambda^{2}$,

$B_{1}(t)=\mathrm{e}^{-\mathrm{ad}_{\Omega_{1}}}(B_{0}-G_{1})=\mathrm{e}^{-\lambda\,\mathrm{ad}_{W_{1}}}C_{1}=\sum_{l\geq
2}\lambda^{l}\,b_{1,l},$ (26)

where

$b_{1,l}=\sum_{k=0}^{l-2}\frac{(-1)^{k}}{k!}\mathrm{ad}_{W_{1}}^{k}\big{(}b_{0,l-k}-g_{1,l-k}\big{)}.$
(27)

We can analogously choose the second transformation proportional to
$\lambda^{2}$, i.e., as $\Omega_{2}(t)\equiv\lambda^{2}W_{2}(t)$ for a given
$W_{2}$ to be determined. Then a straightforward calculation shows that

$C_{2}\equiv B_{1}-G_{2}=\sum_{l\geq 2}\lambda^{l}(b_{1,l}-g_{2,l})$ (28)

with

$g_{2,2}=\dot{W}_{2},\qquad
g_{2,2l}=\frac{1}{l!}\mathrm{ad}_{W_{2}}^{l-1}\dot{W}_{2},\qquad
g_{2,r}=0,\quad r\neq 2l.$ (29)

The generator $W_{2}$ is then obtained by imposing that $b_{1,2}-g_{2,2}=0$ in
Eq. (28), i.e.,

$g_{2,2}=\dot{W}_{2}=b_{1,2}=b_{0,2}-g_{1,2}=-\frac{1}{2}\mathrm{ad}_{W_{1}}\dot{W}_{1}.$
(30)

In this way $B_{2}(t)$ is a power series in $\lambda$ starting with
$\lambda^{3}$,

$B_{2}(t)=\mathrm{e}^{-\lambda^{2}\,\mathrm{ad}_{W_{2}}}C_{2}=\sum_{l\geq
3}\lambda^{l}\,b_{2,l}$ (31)

with

$b_{2,l}=\sum_{k=0}^{[(l-1)/2]-1}\frac{(-1)^{k}}{k!}\mathrm{ad}_{W_{2}}^{k}\big{(}b_{1,l-2k}-g_{2,l-2k}\big{)},$
(32)

where $[\cdot]$ stands for the integer part of the argument. In general, the
$n$-th transformation $\Omega_{n}(t)\equiv\lambda^{n}\,W_{n}(t)$ is determined
in such a way that the power series of $B_{n}(t)$ starts with $\lambda^{n+1}$.
This can be done as follows: from $B_{n-1}=\sum_{l\geq
n}\lambda^{l}b_{n-1,l}$, we compute

$C_{n}\equiv
B_{n-1}-G_{n}=\sum_{r=n}^{\infty}\lambda^{r}c_{n,r}=\sum_{r=n}^{\infty}\lambda^{r}(b_{n-1,r}-g_{n,r})$
(33)

with

$g_{n,r}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{l!}\mathrm{ad}_{W_{n}}^{l-1}\dot{W}_{n},&r=n\,l\\\
0,&r\neq n\,l\end{array}\right.\qquad l=1,2,\ldots$ (34)

Then, $W_{n}$ is obtained by taking $b_{n-1,n}-g_{n,n}=0$, i.e.,

$\dot{W}_{n}=b_{n-1,n}$ (35)

and finally $B_{n}$ is determined as

$B_{n}=\mathrm{e}^{-\lambda^{n}\,\mathrm{ad}_{W_{n}}}C_{n}=\sum_{l=n+1}^{\infty}\lambda^{l}\,b_{n,l}$
(36)

with

$b_{n,l}=\sum_{k=0}^{[\frac{l-1}{n}]-1}\,\frac{(-1)^{k}}{k!}\mathrm{ad}_{W_{n}}^{k}c_{n,l-n\,k}.$
(37)

Notice that, in view of Eq. (34) and Eq. (37), equation (35) simplifies to

$\dot{W}_{n}=b_{n-2,n}\qquad\mbox{ for }\qquad n\geq 3.$ (38)

The solution of Eq. (1) is expressed, after $n$ such transformations, as

$x(t)=\mathrm{e}^{\lambda
W_{1}(t)}\,\mathrm{e}^{\lambda^{2}W_{2}(t)}\,\cdots\,\mathrm{e}^{\lambda^{n}W_{n}(t)}X_{n}(t).$
(39)

An approximation to the exact solution containing all the dependence up to
$\mathcal{O}(\lambda^{n})$ is obtained by taking $X_{n}=x_{0}$ in Eq. (39).

In spite of this, the truncated factorization $\mathrm{e}^{\lambda
W_{1}(t)}\,\mathrm{e}^{\lambda^{2}W_{2}(t)}\,\cdots\,\mathrm{e}^{\lambda^{n}W_{n}(t)}$
still shares relevant qualitative properties, such as orthogonality,
unitarity, etc.

Recursion (33)-(38) allows one to construct _any_ generator $W_{k}$ of the
expansion in terms of $W_{1},\ldots,W_{k-1}$. For illustration, we next
collect the first terms:

$\displaystyle\dot{W}_{1}=A(t),\qquad\quad\dot{W}_{2}=-\frac{1}{2}\mathrm{ad}_{W_{1}}\dot{W}_{1},\qquad\quad\dot{W}_{3}=\frac{1}{3}\mathrm{ad}_{W_{1}}^{2}\dot{W}_{1}$
(40)
$\displaystyle\dot{W}_{4}=-\frac{1}{8}\mathrm{ad}_{W_{1}}^{3}\dot{W}_{1}-\frac{1}{2}\mathrm{ad}_{W_{2}}\dot{W}_{2},\qquad\quad\dot{W}_{5}=\frac{1}{30}\mathrm{ad}_{W_{1}}^{4}\dot{W}_{1}-\mathrm{ad}_{W_{2}}\dot{W}_{3}$

This is the way the Wilcox expansion is built up.

###  3.2 Explicit expressions for $W_{n}(t)$

Although the recursive procedure (33)-(38) turns out to be very
computationally efficient to construct the exponents $W_{n}(t)$ for a given
$A(t)$ in practice, it is clear that much insight about the expansion can be
gained if an explicit expression for _any_ $W_{n}$ can be constructed, thus
generalizing the treatment done originally by Wilcox up to $n=4$ [1].

Such an expression could be obtained, in principle, by working out the
recurrence (33)-(38), but a more direct approach consists in comparing the
Dyson perturbation series of $U(t)$ [22] in the associated initial value
problem

$\dot{U}=\lambda\,A(t)U,\qquad U(0)=I,$ (41)

i.e.,

$U(t)=I+\sum_{n=1}^{\infty}\lambda^{k}P_{k}(t),\qquad\mbox{ with }\qquad
P_{k}(t)=\int_{0}^{t}dt_{1}\cdots\int_{0}^{t_{k-1}}dt_{k}\,A(t_{1})\cdots
A(t_{k})$ (42)

with the expansion in $\lambda$ of the factorization

$U(t)=\mathrm{e}^{\lambda
W_{1}(t)}\,\mathrm{e}^{\lambda^{2}W_{2}(t)}\,\cdots\,\mathrm{e}^{\lambda^{n}W_{n}(t)}\cdots.$
(43)

Thus, for the first terms one has

$\displaystyle P_{1}$ $\displaystyle=W_{1},\qquad\quad
P_{2}=W_{2}+\dfrac{1}{2}W_{1}^{2},\qquad\quad
P_{3}=W_{3}+W_{1}W_{2}+\frac{1}{3!}W_{1}^{3},$ (44) $\displaystyle P_{4}$
$\displaystyle=W_{4}+W_{1}W_{3}+\frac{1}{2}W_{1}^{2}W_{2}+\dfrac{1}{2}W_{2}^{2}+\frac{1}{4!}W_{1}^{4}.$

In general, we can write

$P_{n}=\sum_{p(n)\atop i_{1}\leq i_{2}\leq\cdots\leq
i_{k}}\frac{1}{r_{k}!}\,W_{i_{1}}W_{i_{2}}\cdots W_{i_{k}},\qquad i_{1}\leq
i_{2}\leq\cdots\leq i_{k},$ (45)

where the sum is extended over the total number of partitions $p(n)$ of the
integer $n$. We recall that a partition $p(n)$ of the integer $n$ is an
$n$-tuple $(i_{1},i_{2},\ldots,i_{k})$ such that $i_{1}+i_{2}+\cdots+i_{k}=n$,
with ordering $i_{1}\leq i_{2}\leq\cdots\leq i_{k}$. Thus, the seven
partitions of $n=5$ with the chosen ordering are $(5)$, $(1,4)$, $(2,3)$,
$(1,1,3)$, $(1,2,2)$, $(1,1,1,2)$ and $(1,1,1,1,1)$. In Eq. (45), $r_{k}$ is
the number of repeated indices in the partition considered.

By working out equation (45), one can invert the relations and express $W_{n}$
in terms of $P_{1},\ldots,P_{n-1}$ for any $n\geq 1$. Thus, i one obtains

$\displaystyle W_{1}$ $\displaystyle=P_{1},\qquad\quad
W_{2}=P_{2}-\frac{1}{2}P_{1}^{2},\qquad\quad
W_{3}=P_{3}-P_{1}P_{2}+\frac{1}{3}P_{1}^{3},$ (46) $\displaystyle W_{4}$
$\displaystyle=P_{4}-P_{1}P_{3}+\frac{3}{4}P_{1}^{2}P_{2}+\frac{1}{4}P_{2}P_{1}^{2}-\frac{1}{2}P_{2}^{2}-\frac{1}{4}P_{1}^{4}.$

Notice that $W_{n}$, $n\geq 2$, is expressed in terms of products of iterated
integrals $P_{i_{1}}\cdots P_{i_{j}}$. Interestingly, it is possible to
express these products as proper time-ordered integrals by using a procedure
developed in [23]. If we denote

$A(i_{1}i_{2}\ldots
i_{n})\equiv\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\cdots\int_{0}^{t_{n-1}}dt_{n}\,A(t_{i_{1}})A(t_{i_{2}})\cdots
A(t_{t_{n}}),$ (47)

so that $P_{n}(t)=A(12\ldots n)$, then

$\displaystyle W_{2}=A(12)-\frac{1}{2}A(1)\cdot A(1)$ (48) $\displaystyle
W_{3}=A(123)-A(1)\cdot A(12)+\frac{1}{3}A(1)\cdot A(1)\cdot A(1),$

etc. Taking into account Fubini’s theorem,

$\int_{0}^{\alpha}dy\int_{y}^{\alpha}f(x,y)\,dx=\int_{0}^{\alpha}dx\int_{0}^{x}f(x,y)\,dy,$
(49)

it is clear that $A(1)\cdot A(1)=A(12)+A(21)$, and thus

$W_{2}=A(12)-\frac{1}{2}\big{(}A(12)+A(21)\big{)}=\frac{1}{2}\big{(}A(12)-A(21)\big{)}.$
(50)

We can proceed analogously with the following products

$\displaystyle A(1)\cdot A(12)=A(123)+A(213)+A(312)$ (51) $\displaystyle
A(1)\cdot A(1)\cdot A(1)=A(123)+A(132)+A(213)+A(231)+A(312)+A(321),$

so that

$W_{3}=\frac{1}{3}\,A(123)+\frac{1}{3}\,A(132)-\frac{2}{3}\,A(213)+\frac{1}{3}\,A(231)-\frac{2}{3}\,A(312)+\frac{1}{3}\,A(321).$
(52)

Carrying out this argument to any order, we can expand all the products of
integrals appearing in $W_{n}$. As a result, each product is replaced by the
sum of all possible permutations of time ordering consistent with the time
ordering in the factors of this product [24].

At this point, it is illustrative to consider in detail some examples. Thus,
the product $A(1)\cdot A(12)$ gives the sum of all permutations of three
elements such that the second index is less than the third one. With respect
to $A(1)\cdot A(1)\cdot A(1)$, since there is no special ordering, then all
possible permutations have to be taken into account. Finally, for the product
$P_{1}P_{3}$ appearing in $W_{4}$ one has

$P_{1}P_{3}=A(1)\cdot A(123)=A(4123)+A(3124)+A(2134)+A(1234).$ (53)

Proceeding in a similar way, one can show that any product of iterated
integrals can be expressed as a sum of iterated integrals. This property is in
fact related with a much deeper characterization of the group of permutations
[23]. If $\mathfrak{S}Symm$ denotes the graded $\mathbb{Q}$-vector space with
fundamental basis given by the disjoint union of the symmetric groups $S_{n}$
for all $n\geq 0$, then it is possible to define a product $\ast$ of
permutations and a coproduct $\delta$ in $\mathfrak{S}Symm$ so that there
exists a one-to-one correspondence between iterated integrals and
permutations:

$A(\sigma)\cdot A(\tau)=A(\sigma\ast\tau).$ (54)

The product $\ast$ was introduced in [25], and together with the coproduct
$\delta$ endow $\mathfrak{S}Symm$ with a structure of Hopf algebra [26], the
so-called Malvenuto–Reutenauer Hopf algebra of permutations [27].

In sum, the _general_ structure of the Wilcox expansion terms reads

$W_{n}=\sum_{k=1}^{n!}\omega_{k}^{(n)}A(\sigma_{k}),$ (55)

where the summation extends over all the $n!$ permutations $\sigma_{k}\in
S_{n}$ of $\\{1,2,\ldots,n\\}$. The weights $\omega_{k}^{(n)}$ are given by
rational numbers that can be determined algorithmically for any $n$, although,
it seems not obvious to find a general expression for them, in contrast with
the Magnus expansion, for which such a closed formula exists [28]. This can be
then considered as an open problem.

Moreover, if one is interested in getting a compact expression for $W_{n}$ in
terms of independent nested commutators of $A(t)$, as is done in [1] up to
$n=4$, one can use the class of bases proposed by Dragt & Forest in [24] for
the Lie algebra generated by the operators
$A(t_{1}),A(t_{2}),\ldots,A(t_{n})$. The same procedure carried out in [23]
for the Magnus expansion can be applied here, so that one gets the general
formula

$\displaystyle
W_{n}(t)=\sum_{\tau_{k}}c_{\tau_{k}}^{(n)}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\cdots\int_{0}^{t_{n-1}}dt_{n}\,$
(56)
$\displaystyle\qquad\qquad\qquad[A(t_{\tau(2)}),[A(t_{\tau(3)})\cdots[A(t_{\tau(n)}),A(t_{1})]\cdots]].$

Here the sum extends over the $(n-1)!$ permutations $\tau_{k}$ of the elements
$\\{2,3,\ldots,n\\}$ and $c_{\tau_{k}}^{(n)}$ is a rational number that
depends on the particular permutation. For a given permutation, say
$\tau_{k}$, its value coincides with the prefactor in Eq. (55) of the
particular term $A(\sigma_{k})$ corresponding to the permutation
$\sigma_{k}\in S_{n}$ such that

$\\{\sigma_{k}(1),\sigma_{k}(2),\ldots,\sigma_{k}(n-1),\sigma_{k}(n)\\}=\\{\tau_{k}(2),\tau_{k}(3),\ldots,\tau_{k}(n),1\\}.$
(57)

Thus, if we denote

$A[i_{1}i_{2}\ldots
i_{n}]\equiv\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\cdots\int_{0}^{t_{n-1}}dt_{n}[A(t_{i_{1}}),[A(t_{i_{2}}),\cdots[A(t_{i_{n-1}}),A(t_{i_{n}})]\cdots]]$
(58)

we get for the first terms

$\displaystyle W_{2}=-\frac{1}{2}A[21],$ (59) $\displaystyle
W_{3}=\frac{1}{3}A[231]+\frac{1}{3}A[321],$ $\displaystyle
W_{4}=-\frac{1}{4}A[3241]-\frac{1}{4}A[4231]-\frac{1}{4}A[4321],$
$\displaystyle
W_{5}=-\frac{2}{15}A[23451]-\frac{2}{15}A[23541]-\frac{2}{15}A[24351]-\frac{2}{15}A[24531]$
$\displaystyle\qquad-\frac{2}{15}A[25341]-\frac{2}{15}A[25431]+\frac{1}{5}A[32451]+\frac{1}{5}A[32541]-\frac{2}{15}A[34251]$
$\displaystyle\qquad-\frac{2}{15}A[34521]-\frac{2}{15}A[35241]-\frac{2}{15}A[35421]+\frac{1}{5}A[42351]+\frac{1}{5}A[42531]$
$\displaystyle\qquad+\frac{1}{5}A[43251]+\frac{1}{5}A[43521]-\frac{2}{15}A[45231]-\frac{2}{15}A[45321]+\frac{1}{5}A[52341]$
$\displaystyle\qquad+\frac{1}{5}A[52431]+\frac{1}{5}A[53241]+\frac{1}{5}A[53421]+\frac{1}{5}A[54231]+\frac{1}{5}A[54321]$

In sum, the general structure of the Wilcox expansion terms via commutators
uses the same weights as in Eq. (55) and reads

$W_{n}=\sum_{k=1}^{(n-1)!}{\vphantom{\sum}}^{\prime}\omega_{k}^{(n)}A[\sigma_{k}],$
(60)

where the primed sum requires the rightmost element in the permutation
$\sigma(k)$ to be invariant (as in Eq. (56)). This element may be chosen at
will and, whatever it be that value, the permutations are build up with the
remaining $n-1$ elements. Different, but equivalent, expressions for $W_{n}$
in terms of commutators are obtained depending on the value fixed at the
rightmost position. We stress once again that, although only the first terms
have been collected here for simplicity, the whole procedure is algorithmic in
nature and has been implemented in a computer algebra system furnishing to
evaluate explicitly $W_{n}(t)$ for any $n$ [29]. Notice that $W_{n}(t)$
involves a linear combination of $(n-1)!$ iterated integrals.

###  3.3 Convergence of Wilcox expansion

Recursion (33)-(38) is also very useful to provide estimates for the radius of
convergence of the Wilcox expansion when Eq. (1) is defined in a Banach
algebra $\mathcal{A}$, i.e., an algebra that is also a complete normed linear
space with a sub-multiplicative norm,

$\|X\,Y\|\leq\|X\|\,\|Y\|.$ (61)

If this is the case, then $\|\mathrm{ad}_{X}Y\|\leq 2\,\|X\|\,\|Y\|$ and, in
general, $\|\mathrm{ad}_{X}^{n}Y\|\leq 2^{n}\|X\|^{n}\,\|Y\|$.

As shown in [30, 31], if the series

$M(\lambda;t)=\sum_{j=1}^{\infty}\lambda^{j}\,\|W_{j}(t)\|$ (62)

has a certain radius of convergence $r_{c}$ for a given $t$, then, for
$\lambda<r<r_{c}$, the sequence of functions

$\Psi_{n}\equiv\mathrm{e}^{\lambda
W_{1}(t)}\,\mathrm{e}^{\lambda^{2}W_{2}(t)}\,\cdots\,\mathrm{e}^{\lambda^{n}W_{n}(t)}$
(63)

converges uniformly on any compact subset of the ball $B(0,r_{c})$. Thus,
studying the convergence of the Wilcox expansion reduces to analyze the series
$M(\lambda;t)$ and in particular its radius of convergence $r_{c}$.

Let $k(t)$ be a function such that $\|A(t)\|\leq k(t)$ and denote
$K(t)=\int_{0}^{t}k(s)\,ds$. Then, clearly

$\|b_{0,1}\|\leq k(t);\qquad\quad\|b_{0,l}\|=0,\quad l=2,3,\ldots$ (64)

and

$\|W_{1}\|\leq K(t),\qquad\|W_{2}\|\leq\frac{1}{2}K^{2}(t).$ (65)

In general, the following bounds can be established by induction:

$\displaystyle\|g_{n,r}(t)\|$ $\displaystyle\leq$
$\displaystyle\beta_{n,r}\,K^{r-1}(t)\,k(t),$ $\displaystyle\|b_{n,l}(t)\|$
$\displaystyle\leq$ $\displaystyle\alpha_{n,l}\,K^{l-1}(t)\,k(t),\qquad\qquad
n=1,2,\ldots;\quad l>n$ $\displaystyle\|W_{n}(t)\|$ $\displaystyle\leq$
$\displaystyle c_{n}\,K^{n}(t),$ (66)

where

$\displaystyle\alpha_{0,1}$ $\displaystyle=$ $\displaystyle
1;\qquad\alpha_{0,l}=0,\;\;l>1$ $\displaystyle\alpha_{n,l}$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{[\frac{l-1}{n}]-1}\frac{1}{j!}\,2^{j}\,c_{n}^{j}\,(\alpha_{n-1,l-nj}+\beta_{n,l-nj})\quad\qquad
n\geq 1,\;l>n$ $\displaystyle\beta_{n,r}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{l!}\,2^{l-1}\,n\,c_{n}^{l},&\left[\frac{r}{n}\right]=\frac{r}{n}=l\\\
0,&\left[\frac{r}{n}\right]\neq\frac{r}{n}\end{array}\right.$ (69)
$\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle 1,\quad
c_{2}=\frac{1}{2},\qquad c_{n}=\frac{1}{n}\,\alpha_{n-2,n},\qquad n>2.$

It is clear that if the series $\sum_{j\geq 1}c_{j}\,K^{j}(t)$ converges, so
does $M(\lambda=1;t)$. Therefore, a sufficient condition for convergence of
the Wilcox expansion is obtained by imposing

$\lim_{n\rightarrow\infty}\frac{c_{n+1}K^{n+1}(t)}{c_{n}K^{n}(t)}=K(t)\lim_{n\rightarrow\infty}D_{n}<1,$
(70)

where

$D_{n}\equiv\frac{n}{n+1}\frac{\alpha_{n-1,n+1}}{\alpha_{n-2,n}}.$ (71)

We have computed this quantity up to $n=2000$ and then extrapolated to the
limit $(1/n)\rightarrow 0$. Then $D_{n}\rightarrow D_{\infty}=1.51868$, as
seen in Figure 1, and thus the convergence of the Wilcox expansion is ensured
at least for values of time $t$ such that

$\int_{0}^{t}\|A(s)\|ds\leq K(t)<\xi_{W}=\frac{1}{D_{\infty}}\approx 0.65846$
(72)

This type of extrapolation has also been used to estimate the convergence
radius of the Magnus expansion [32]. Although the estimate is not completely
analytic, the same type of computation has provided accurate results in other
settings. In particular, for the Magnus expansion such an estimate fully
agrees with a purely theoretically deduced bound [10, 8, 32].

Figure 1: $D_{n}$ as a function of $1/n$, and linear extrapolation (red line).

##  4 Fer–like expansions

###  4.1 Standard Fer expansion

In forming the Wilcox expansion the first transformation is chosen in such a
way that $\dot{\Omega}_{1}=B_{0}(t)$, whereas $\dot{\Omega}_{n}\neq B_{n-1}$
for $n\geq 2$. It makes sense, then, to analyze what happens if we impose this
condition at each step of the procedure

$\dot{\Omega}_{n}(t)=B_{n-1}(t)\qquad\mbox{ or, equivalently,
}\qquad\Omega_{n}(t)=\int_{0}^{t}B_{n-1}(u)du$ (73)

for all $n\geq 1$. In that way, expression (22) for $B_{n}$ clearly simplifies
to

$B_{n}(t)=\sum_{k\geq
1}\frac{(-1)^{k}k}{(k+1)!}\,\mathrm{ad}_{\Omega_{n}}^{k}B_{n-1},\qquad n\geq
1.$ (74)

In doing so we recover precisely the Fer expansion, see [11, 10]. Again, after
$n$ transformations, we get

$x(t)=\mathrm{e}^{\Omega_{1}(t)}\,\mathrm{e}^{\Omega_{2}(t)}\cdots\mathrm{e}^{\Omega_{n}(t)}X_{n}(t)$
(75)

so that if we impose $X_{n}(t)=x_{0}$ we are left with another approximation
to the exact solution. Notice that this approximation clearly differs from the
previous Wilcox expansion for $n\geq 2$, as can be seen by analyzing the
dependence on $\lambda$ of each transformation. Whereas
$\Omega_{n}=\lambda^{n}W_{n}(t)$ for the Wilcox expansion, now $\Omega_{n}$
contains terms of order $\lambda^{2^{n-1}}$ and higher. This can be easily
shown by induction: $\Omega_{1}$ is proportional to $\lambda$, so that
$B_{1}$, according to Eq. (74) contains terms of order $\lambda^{2}$ (coming
from $[\Omega_{1},B_{0}]$) and higher. In general, $B_{n-2}$ and
$\Omega_{n-1}$ contain terms of order $\lambda^{2^{n-2}}$ and higher, so the
first term in the series (74) for $B_{n-1}$, i.e., the commutator
$[\Omega_{n-1},B_{n-2}]$, produces a term of order
$(\lambda^{2})^{2^{n-2}}=\lambda^{2^{n-1}}$ in $\Omega_{n}$.

Alternatively, expressing Eq. (74) as

$B_{n}(t)=\int_{0}^{1}dx\int_{0}^{x}du\,\mathrm{e}^{-(1-u)\Omega_{n}}\,[B_{n-1},\Omega_{n}]\,\mathrm{e}^{(1-u)\Omega_{n}}$
(76)

and taking norms, it is then possible to show that the Fer expansion converges
for values of $t$ such that [10]

$\int_{0}^{t}\|A(s)\|ds<0.8604065.$ (77)

###  4.2 Intermediate Fer-like expansions

Notice that the $\lambda$-power series of $\Omega_{n}$ in the Fer expansion
contains infinite terms starting with $\lambda^{2^{n-1}}$, but the
corresponding truncated factorization obtained from Eq. (75) by taking
$X_{n}(t)=x_{0}$ is correct only up to terms of order $\lambda^{2^{n-1}}$. One
might then consider yet another sequence of transformations so that each
$\Omega_{k}$ contains only those relevant terms leading to a correct
approximation up to this order. Of course, both factorizations would be
different, but nevertheless they would produce the correct power series up to
order $\lambda^{2^{n-1}}$. The corresponding factorization can be properly
called a modified Fer expansion.

Our starting point is once again equation (22). Clearly, the first
transformation is the same as in Fer (and Wilcox), i.e.,

$\dot{\Omega}_{1}=B_{0}=\lambda A(t),$ (78)

and thus

$B_{1}=\sum_{k=0}^{\infty}\frac{(-1)^{k}k}{(k+1)!}\mathrm{ad}_{\Omega_{1}}^{k}B_{0}=\mathcal{O}(\lambda^{2}),$
(79)

where the rightmost term points out the lowest $\lambda$ contribution in the
sum.

Next, to reproduce the same dependence on $\lambda$ as the Fer expansion, we
need to enforce that $B_{2}=\mathcal{O}(\lambda^{4})$, and the question is how
to choose $\Omega_{2}$ guaranteeing this feature. An analysis of Eq. (22) with
$n=2$ reveals that this is achieved by taking $\dot{\Omega}_{2}$ as the sum of
terms in $B_{1}$ in Eq. (79) contributing to $\lambda^{2}$ and $\lambda^{3}$,
i.e.,

$\dot{\Omega}_{2}=-\frac{1}{2}\mathrm{ad}_{\Omega_{1}}B_{0}+\frac{1}{3}\mathrm{ad}_{\Omega_{1}}^{2}B_{0},$
(80)

since the next term appearing in the expression of $B_{2}$ involves the
computation of $\mathrm{ad}_{\Omega_{1}}^{3}B_{0}=\mathcal{O}(\lambda^{4})$.
Thus

$B_{2}=\sum_{k=1}^{\infty}\frac{(-1)^{k}}{(k+1)!}\left((k+1)\mathrm{ad}_{\Omega_{2}}^{k}B_{1}-\mathrm{ad}_{\Omega_{2}}^{k}\dot{\Omega}_{2}\right)=\mathcal{O}(\lambda^{4}).$
(81)

Likewise, $\Omega_{3}$ is to be designed so that

$B_{3}=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(k+1)!}\left((k+1)\mathrm{ad}_{\Omega_{3}}^{k}B_{2}-\mathrm{ad}_{\Omega_{3}}^{k}\dot{\Omega}_{3}\right)=\mathcal{O}(\lambda^{8}),$
(82)

and this is guaranteed by taking $\Omega_{3}$ as the sum all the terms in
$B_{2}$ contributing to powers from $\lambda^{4}$ up to $\lambda^{7}$. From
Eq. (81) it is clear that

$B_{2}=\underbrace{-\frac{1}{2}\left(2\mathrm{ad}_{\Omega_{2}}B_{1}-\mathrm{ad}_{\Omega_{2}}\dot{\Omega}_{2}\right)}_{\mathcal{O}(\lambda^{4})}+\underbrace{\frac{1}{3!}\left(3\mathrm{ad}_{\Omega_{2}}^{2}B_{1}-\mathrm{ad}_{\Omega_{2}}^{2}\dot{\Omega}_{2}\right)}_{\mathcal{O}(\lambda^{6})}+\mathcal{O}(\lambda^{8}),$
(83)

where only the relevant terms in the expansion in $B_{1}$ have to be taken
into account. In this way, we can take

$\dot{\Omega}_{3}=-\mathrm{ad}_{\Omega_{2}}B_{1}^{[4]}+\frac{1}{2}\mathrm{ad}_{\Omega_{2}}\dot{\Omega}_{2}+\frac{1}{2}\mathrm{ad}_{\Omega_{2}}^{2}B_{1}^{[2]}-\frac{1}{6}\mathrm{ad}_{\Omega_{2}}^{2}\dot{\Omega}_{2}$
(84)

with

$B_{1}^{[j]}\equiv\sum_{k=1}^{j}\frac{(-1)^{k}}{(k+1)!}k\,\mathrm{ad}_{\Omega_{1}}^{k}B_{0}.$
(85)

Notice, however, that since the second term in $\Omega_{2}$ in Eq. (80) is
$\mathcal{O}(\lambda^{3})$, the expression (84) does contain some
contributions in $\lambda^{8}$ and $\lambda^{9}$ that in principle could be
removed. We prefer, however, to maintain them just to have a more compact
expression.

For this modified Fer expansion, $\Omega_{n}$ is chosen in general so that
$\dot{\Omega}_{n}$ is precisely the sum of all terms of $B_{n-1}$ containing
terms of powers from $\lambda^{2^{n-1}}$ up to $\lambda^{2^{n}-1}$ and then
truncating appropriately the series of $B_{n-2},\ldots,B_{1}$.

Other possibilities for choosing $B_{k}$ at the successive stages clearly
exist, and according with the particular election, different intermediate Fer-
like expansions result. In practice, one of those combinations of commutators
could be more easily computed for a specific problem.

##  5 Applications

###  5.1 Wilcox expansion as the continuous analogue of Zassenhaus formula

The Zassenhaus formula may be considered as the dual of the
Baker–Campbell–Hausdorff (BCH) formula [33] in the sense that it relates the
exponential of the sum of two non-commuting operators $X$ and $Y$ with an
infinite product of exponentials of these operators and their nested
commutators. More specifically,

$\mathrm{e}^{X+Y}=\mathrm{e}^{X}\,\mathrm{e}^{Y}\,\prod_{n=2}^{\infty}\mathrm{e}^{C_{n}(X,Y)}=\mathrm{e}^{X}\,\mathrm{e}^{Y}\,\mathrm{e}^{C_{2}(X,Y)}\,\mathrm{e}^{C_{3}(X,Y)}\,\cdots\,\mathrm{e}^{C_{k}(X,Y)}\,\cdots,$
(86)

where $C_{k}(X,Y)$ is a homogeneous Lie polynomial in $X$ and $Y$ of degree
$k$ [21, 7, 34, 35, 1]. A very efficient procedure to generate all the terms
in Eq. (86) is presented in [21] and allows to construct $C_{n}$ up to a
prescribed value of $n$ directly in terms of the minimum number of independent
commutators involving $n$ operators $X$ and $Y$.

In view of the formal similarity between Eq. (39) and Eq. (86), Wilcox
expansion has been also described as the “continuous analogue of the
Zassenhaus formula” [10], just as the Magnus expansion is sometimes called the
continuous version of the BCH formula. To substantiate this claim, we next
reproduce the Zassenhaus formula (86) by applying the procedure of Section 3
to a particular initial value problem, namely the abstract equation

$\dot{U}=\lambda(X+Y)\,U,\qquad U(0)=I,$ (87)

where $X$ and $Y$ are two non-commuting constant operators.

The formal solution is of course $U(t)=\mathrm{e}^{t\lambda(X+Y)}$, but we can
also solve Eq. (87) by first integrating $\dot{U}_{0}=\lambda XU_{0}$ and
factorizing $U(t)$ as $U(t)=U_{0}\,U_{I}=\mathrm{e}^{t\lambda X}\,U_{I}$,
where $U_{I}$ obeys the equation

$\dot{U}_{I}=\lambda\,\mathrm{e}^{-t\lambda X}Y\mathrm{e}^{t\lambda
X}\,U_{I}\equiv A_{\lambda}(t)\,U_{I},$ (88)

and finally apply to Eq. (88) the sequence of transformations leading to the
Wilcox expansion. Notice, however, that now the coefficient matrix
$A_{\lambda}(t)$ is an infinite series in $\lambda$,

$A_{\lambda}(t)=\lambda\,\mathrm{e}^{-t\lambda\mathrm{ad}_{X}}Y=\sum_{j\geq
0}\frac{(-1)^{j}}{j!}t^{j}\lambda^{j+1}\mathrm{ad}_{X}^{j}Y,$ (89)

so that, when applying the recursion (33)-(38), $\dot{\Omega}_{1}$ is no
longer $B_{0}\equiv A_{\lambda}(t)$, but the term in $A_{\lambda}(t)$ which is
proportional to $\lambda$. In other words,

$\dot{W}_{1}=Y,\qquad\mbox{ and thus }\qquad W_{1}(t)=t\,Y.$ (90)

After some computation, one arrives at

$b_{0,l}=\frac{(-1)^{l-1}}{(l-1)!}t^{l-1}\mathrm{ad}_{X}^{l-1}Y.$ (91)

Since $\dot{W}_{1}=b_{0,1}=Y$, then clearly $g_{1,l}=0$ for all $l>1$ and

$b_{1,l}=\sum_{k=0}^{l-2}\frac{(-1)^{k}}{k!}\mathrm{ad}_{W_{1}}^{k}b_{0,l-k}.$
(92)

By imposing $b_{1,2}=g_{2,2}=\dot{W}_{2}$, we get

$\dot{W}_{2}=-t\,\mathrm{ad}_{X}Y,\qquad\mbox{ so that }\qquad
W_{2}(t)=-\frac{1}{2}\,t^{2}\,\mathrm{ad}_{X}Y.$ (93)

In general, $g_{n,n}=\dot{W}_{n}$, $g_{n,r}=0$ when $r\neq n$ and the
recurrence (33)-(38) reads now

$\displaystyle\dot{W}_{1}$ $\displaystyle=$ $\displaystyle b_{0,1}$
$\displaystyle b_{n,l}$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{[\frac{l-1}{n}]-1}\,\frac{(-1)^{k}}{k!}\mathrm{ad}_{W_{n}}^{k}b_{n-1,l-nk},\qquad
n=1,2,\ldots$ (94) $\displaystyle\dot{W}_{n}$ $\displaystyle=$ $\displaystyle
b_{n-2,n}\qquad n\geq 2$

together with Eq. (91). Working out this recursion we obtain, for the first
terms

$\displaystyle W_{3}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{6}t^{3}\mathrm{ad}_{X}^{2}Y+\frac{1}{3}t^{3}\mathrm{ad}_{Y}\mathrm{ad}_{X}Y$
$\displaystyle W_{4}(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{24}t^{4}\mathrm{ad}_{X}^{3}Y-\frac{1}{8}t^{4}\mathrm{ad}_{Y}\mathrm{ad}_{X}^{2}Y-\frac{1}{8}t^{4}\mathrm{ad}_{Y}^{2}\mathrm{ad}_{X}Y$
$\displaystyle W_{5}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{120}t^{5}\mathrm{ad}_{X}^{4}Y+\frac{1}{30}t^{5}\mathrm{ad}_{Y}\mathrm{ad}_{X}^{3}Y+\frac{1}{20}t^{5}\mathrm{ad}_{Y}^{2}\mathrm{ad}_{X}^{2}Y$
$\displaystyle+\frac{1}{30}t^{5}\mathrm{ad}_{Y}^{3}\mathrm{ad}_{X}Y+\frac{1}{20}t^{5}\mathrm{ad}_{[X,Y]}\mathrm{ad}_{X}^{2}Y+\frac{1}{10}t^{5}\mathrm{ad}_{[X,Y]}\mathrm{ad}_{Y}\mathrm{ad}_{X}Y$

One can see that this procedure agrees with the algorithm presented in [21]
for every term $W_{n}$, $n\geq 1$, in

$U(t)=\mathrm{e}^{t\lambda(X+Y)}=\mathrm{e}^{t\lambda X}\,\mathrm{e}^{\lambda
W_{1}(t)}\,\mathrm{e}^{\lambda^{2}W_{2}(t)}\,\mathrm{e}^{\lambda^{3}W_{3}(t)}\,\cdots.$
(96)

The Zassenhaus formula is recovered by taking $t=1$, i.e.,
$C_{n}(X,Y)=W_{n}(t=1)$.

###  5.2 Expanding the exponential $\exp(A+\varepsilon B)$

Bellman, in his classic book [36], states that “one of the great challenges of
modern physics is that of obtaining useful approximate relations for
$\mathrm{e}^{(A+\varepsilon B)t}$ in the case where $AB\neq BA$”. One such
approximation was proposed and left undisclosed in [12, p. 175]. Assuming that
$\mathrm{e}^{A+\varepsilon B}$ can be written in the form

$\mathrm{e}^{A+\varepsilon B}=\mathrm{e}^{A}\,\mathrm{e}^{\varepsilon
C_{1}}\,\mathrm{e}^{\varepsilon^{2}C_{2}}\,\mathrm{e}^{\varepsilon^{3}C_{3}}\cdots,$
(97)

Bellman proposed to determine the first three terms $C_{1}$, $C_{2}$, $C_{3}$,
and pointed out that, contrary to other expansions, the product expansion (97)
is unitary if $A$ and $B$ are skew-Hermitian.

It turns out that the Wilcox expansion can be used to provide explicit
expressions for $C_{n}$ for any two indeterminates $A$ and $B$, as we will see
in the sequel.

Before proceeding, it is important to remark that this problem differs from
the Zassenhauss formula, in the sense that the expansion parameter affects
only one of the operators in the exponential. The solution goes as follows. We
write

$U(t)\equiv\mathrm{e}^{t(A+\varepsilon B)}=\mathrm{e}^{tA}\,V,$ (98)

and solve the differential equation satisfied by $V$

$\frac{dV}{dt}=\varepsilon
e^{-tA}Be^{tA}\,V\equiv\varepsilon\tilde{B}(t)V,\qquad V(0)=I$ (99)

with the Wilcox expansion, so that

$V(t)=\mathrm{e}^{\varepsilon
W_{1}(t)}\,\mathrm{e}^{\varepsilon^{2}W_{2}(t)}\,\mathrm{e}^{\varepsilon^{3}W_{3}(t)}\cdots.$
(100)

The operators $C_{i}$ in Eq. (97) are then obtained by taking $t=1$.

The successive terms $W_{j}(t)$ in Eq. (100) can be determined either by the
recursion (33)-(38) or the explicit expression (56). For the first term we get

$W_{1}(t)=\int_{0}^{t}\tilde{B}(s)dt=\int_{0}^{t}\mathrm{e}^{-s\,\mathrm{ad}_{A}}B\,ds=\sum_{k=0}^{\infty}\frac{(-1)^{k}\,t^{k+1}}{(k+1)!}\mathrm{ad}_{A}^{k}B$
(101)

and so the expression for $C_{1}$ is given by

$C_{1}=\frac{1-\mathrm{e}^{\mathrm{ad}_{A}}}{\mathrm{ad}_{A}}B=B-\frac{1}{2}[A,B]+\frac{1}{3!}[A,[A,B]]-\frac{1}{4!}[A,[A,[A,B]]]+\cdots.$
(102)

Although it is possible in principle to construct explicit expressions for
$W_{2}$, $W_{3}$, etc., it is perhaps more convenient to apply the recursion
(33)-(38) for each particular application.

###  5.3 Illustrative examples

We next particularize the Bellman problem (97) to matrices where closed
expressions for $C_{1},C_{2},\ldots$ can be obtained. The idea is to
illustrate the behaviour of the product expansion by computing explicitly high
order terms with matrices in the $\mathrm{SU(2)}$ and the $\mathrm{SO(3)}$ Lie
algebras.

#### Matrices $X$ and $Y$ in $\mathrm{SU(2)}$.

In the first example we choose $A=i\,a\sigma_{z}$ and $B=i\,\sigma_{x}$, where
$i=\sqrt{-1}$, $a$ is a real parameter and

$\sigma_{x}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),\quad\sigma_{y}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right),\quad\sigma_{z}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)$ (103)

are Pauli matrices. This instance is borrowed from Quantum Mechanics, where
$\exp[i(a\sigma_{z}+\varepsilon\sigma_{x})]$ is a matrix that transforms the
$\frac{1}{2}$-spin wave function in a Hilbert space.

Using the scalar vector notation to write down a linear combination of Pauli
matrices:
$\vec{v}\cdot\vec{\sigma}=v_{x}\sigma_{x}+v_{y}\sigma_{y}+v_{z}\sigma_{z}$,
the matrix exponential reads

$\exp(i\vec{v}\cdot\vec{\sigma})=\cos v\,I+i\,\frac{\sin
v}{v}\,\vec{v}\cdot\vec{\sigma}.$ (104)

In the sequel we work out the expansion

$\mathrm{e}^{i(a\sigma_{z}+\varepsilon\sigma_{x})}=\mathrm{e}^{ia\sigma_{z}}\,\mathrm{e}^{i\varepsilon
W_{1}}\,\mathrm{e}^{i\varepsilon^{2}W_{2}}\,\mathrm{e}^{i\varepsilon^{3}W_{3}}\ldots$
(105)

up to order eleven in $\varepsilon$ and analyze the increasing accuracy of the
product expansion as far as more terms are considered. The lhs in Eq. (105)
may be thought as a transformation involving $\sigma_{x}$ and $\sigma_{z}$. In
turn, the rhs is a pure $\sigma_{z}$ transformation, i.e.
$\exp(ia\sigma_{z})$, followed by an infinite succession of transformations,
$\exp(i\varepsilon^{k}W_{k})$, whose effect should decrease with $k$. The
truncated product expansion is expected to be accurate as far as
$\varepsilon\ll a$.

$k$ | $W_{k}(t)$
---|---
1 | $\dfrac{1}{2a}[S\sigma_{x}+(1-C)\sigma_{y}]$
2 | $\dfrac{1}{4a^{2}}(2at-S)\sigma_{z}$
3 | $\dfrac{1}{12a^{3}}\left\\{[6at+(C-4)S]\sigma_{x}-(1-C)^{2}\sigma_{y}\right\\}$
4 | $-\dfrac{1}{16a^{4}}[6at+(C-4)S]\sigma_{z}$
5 | $\dfrac{1}{240a^{5}}\left\\{[56S-(4C+7)SC-10at(7+4C-2C^{2})]\sigma_{x}+\right.$
|
$\phantom{\dfrac{i}{240a^{5}}}\left.[4C^{3}-7C^{2}-28C+31+2at(2SC-4S+3at)]\sigma_{y}\right\\}$
Table 1: First five orders in Bellman problem for generic $t$. The operators
$C_{k}$ in Eq. (97) are obtained by taking $t=1$, i.e., $C_{k}=i\,W_{k}(1)$.
We have defined $S\equiv\sin(2at)$ and $C\equiv\cos(2at)$.

In Table 1 we write down explicitly the first five contributions for a generic
$t$ (expressions for $k>5$ are too involved to be collected here).
Wilcox–Bellman’s formula (105) corresponds then to $t=1$. All the terms have
been obtained with the recurrences of section 3, starting from

$\dot{W}_{1}(t)=\mathrm{e}^{-iat\sigma_{z}}\,\sigma_{x}\,\mathrm{e}^{iat\sigma_{z}}=C\sigma_{x}+S\sigma_{y},$
(106)

where $C\equiv\cos(2at)$ and $S\equiv\sin(2at)$.

The formulas in Table 1 show that $\varepsilon/a$ may be considered as an
effective expansion parameter. In Figure 2 we illustrate, for $a=1$, the
accuracy of the Wilcox–Bellman product expansion in the example at hand as a
function of $\varepsilon$. We plot the squared modulus of the non–diagonal
matrix element, say $|U_{1,2}|^{2}$, of Eq. (105) for every analytic
approximation up to order eleven in $\varepsilon$, as well as the exact
result. Even orders do not contribute in this test because $W_{2k}$ is always
proportional to $\sigma_{z}$ and therefore $\exp(i\varepsilon^{2k}W_{2k})$ is
a diagonal matrix.

As regards convergence of the product expansion, the lower bound of Eq. (72)
leads to

$\int_{0}^{1}\|\mathrm{e}^{-iat\sigma_{z}}\,\varepsilon\sigma_{x}\,\mathrm{e}^{iat\sigma_{z}}\|\,\mathrm{d}t=\varepsilon<0.658.$
(107)

In turn, the behaviour of the curves in Figure 2 points out that convergence
of the product expansion extends well beyond that lower bound for this
particular example.

Figure 2: Accuracy of Wilcox–Bellman product expansion up to order eleven as a
function of the ratio $\varepsilon/a$, with $a=1$. The quantity plotted is the
squared modulus of the non–diagonal element of the matrix. The vertical grey
line stands for the convergence lower bound $\varepsilon=0.658$. Figure 3:
Absolute error of the approximations given by curves in Figure 2 with $a=1$.
The vertical grey line is located at the value of the convergence lower bound
$\varepsilon=0.658$.

Eventually, Figure 3 shows the logarithm of the absolute error in the
approximations given by curves in Figure 2.

#### Matrices $X$ and $Y$ in $\mathrm{SO(3)}$.

The second example refers to the matrix that describes a rotation in three
dimensions defined by the vector $\vec{a}=\alpha\hat{a}$. Here $\alpha$ stands
for the rotation angle around the axis given by the unitary vector $\hat{a}$.
A generic 3D rotation matrix can be written as $\exp(\vec{a}\cdot\vec{\rho})$,
where the components of $\vec{\rho}$ are the three fundamental rotation
matrices

$\rho_{x}=\left(\begin{matrix}0&0&0\\\\[-11.38109pt] 0&0&-1\\\\[-11.38109pt]
0&1&0\end{matrix}\right),\qquad\rho_{y}=\left(\begin{matrix}0&0&1\\\\[-11.38109pt]
0&0&0\\\\[-11.38109pt]
-1&0&0\end{matrix}\right),\qquad\rho_{z}=\left(\begin{matrix}0&-1&0\\\\[-11.38109pt]
1&0&0\\\\[-11.38109pt] 0&0&0\end{matrix}\right).$ (108)

We study the particular case
$\vec{a}\cdot\vec{\rho}=\alpha\big{(}\cos\theta\,\rho_{z}+\sin\theta\,\rho_{x}\big{)}$,
and compare the rotation of angle $\alpha$ around the unit vector
$(\sin\theta,0,\cos\theta)$

$\mathrm{e}^{\alpha\big{(}\cos\theta\,\rho_{z}+\sin\theta\,\rho_{x}\big{)}}$
(109)

with the sequence of transformations

$\mathrm{e}^{\alpha\cos\theta\,\rho_{z}}\,\mathrm{e}^{\alpha\sin\theta\,W_{1}}\,\mathrm{e}^{(\alpha\sin\theta)^{2}W_{2}}\,\mathrm{e}^{(\alpha\sin\theta)^{3}W_{3}}\ldots$
(110)

In other words, the question we address is how the pure $z$–axis rotation in
Eq. (110), $\exp(\alpha\cos\theta\,\rho_{z})$, has to be corrected by an
infinite composition of further rotations to reproduce the one defined by
$\vec{a}\cdot\vec{\rho}$. When $\alpha\sin\theta$ is small enough the approach
is expected to converge since the expansion convergence lower bound reads in
this case $|\alpha\sin\theta|<0.658$.

Here the accuracy of the product expansion will depend on both the rotation
angle $\alpha$ and the relative orientation of the rotation axis, determined
by the angle $\theta$. This is illustrated in Figures 4 and 5 for the first
five orders of approximation, with $\alpha=\pi/4,\pi/2,3\pi/4$ and $\pi$.

In order to test the expansion, we have computed the matrix trace of the
successive approximants and compared with the exact result:

$\mathrm{tr}\left(\exp\left(\alpha\big{(}\cos\theta\,\rho_{z}+\sin\theta\,\rho_{x}\big{)}\right)\right)=1+2\cos\alpha.$
(111)

The first two approximants are simple enough to be written down:

$\displaystyle\mathrm{tr}\left(\mathrm{e}^{\alpha\cos\theta\,\rho_{z}}\,\mathrm{e}^{\alpha\sin\theta\,W_{1}}\right)=\Big{(}1+\cos(\alpha\cos\theta)\Big{)}\cos\left(2\tan\theta\,\sin\frac{\alpha\cos\theta}{2}\right)$
(112)
$\displaystyle\mathrm{tr}\left(\mathrm{e}^{\alpha\cos\theta\,\rho_{z}}\,\mathrm{e}^{\alpha\sin\theta\,W_{1}}\,\mathrm{e}^{(\alpha\sin\theta)^{2}W_{2}}\right)=\cos\left(\frac{1}{2}\tan^{2}\theta\,(\sin(\alpha\cos\theta)-\alpha\cos\theta)\right)\cdot$
$\displaystyle\phantom{\mathrm{tr}\left(\mathrm{e}^{\alpha\cos\theta\,\rho_{z}}\,\mathrm{e}^{\alpha\sin\theta\,W_{1}}\,\mathrm{e}^{(\alpha\sin\theta)^{2}W_{2}}\right)=}\cos\left(\alpha\cos\theta-\frac{1}{2}\tan\theta\sin\frac{\alpha\cos\theta}{2}\right)$

Interestingly, the third approximation order is worse than the second one in
all four cases. Also in the case $\alpha=\pi$ the fourth order is better than
the fifth one.

Figure 4: Error in the approximations to the matrix trace as a function of the
rotation angle $\theta$, for two values $\alpha=\pi/2$ and $\pi/4$. Figure 5:
Error in the approximations to the matrix trace as a function of the rotation
angle $\theta$, for two values $\alpha=3\pi/4$ and $\pi$.

##  6 Conclusions

When a linear system of differential equations, defined by the coefficient
matrix $\lambda A(t)$, is transformed under
$\exp(\lambda\int_{0}^{t}A(t^{\prime})\mathrm{d}t^{\prime})$, the coefficient
matrix in the new representation becomes an infinite power series in
$\lambda$, say $\tilde{A}(t)=\sum_{k\geq 1}\lambda^{k}\,a_{k}(t)$. That is the
first step of all matrix exponential methods to approximate the time-evolution
operator. In the framework that we have introduced it is the first move in a
sequence of exponential transformations that change the linear system from a
representation to another with the goal that dynamics will become less and
less relevant. Choosing the transformation
$\exp(\lambda\int_{0}^{t}\tilde{A}(x)\mathrm{d}x)$ as the second move and
iterating this procedure afterwards yields the Fer expansion. Instead,
choosing the transformation $\exp(\lambda\int_{0}^{t}a_{1}(x)\mathrm{d}x)$,
i.e. the leading term of the new coefficient matrix, opens up Wilcox
expansion. The framework allows intermediate expansions taking
$\exp(\sum_{k=1}^{n}\lambda^{k}\int_{0}^{t}a_{k}(x)\mathrm{d}x)$ as initiator,
as well as jumping between schemes, according with the particular requirements
of the problem at hand.

We have seen that the theory of linear transformations (or changes of picture
in the language of Quantum Mechanics) provides a unified framework to deal
with many different exponential perturbative expansions: whereas only one
linear transformation reducing the dynamics to the trivial equation (4) or to
a system with a constant matrix renders the Magnus [9] and the Floquet-Magnus
[37] expansion, respectively, a sequence of such transformations with
different choices of the new matrices lead to Wilcox and Fer factorizations.
From this perspective, other factorizations are possible depending on the
particular problem at hand: one only has to select appropriately the
successive transformations.

In the case of Wilcox expansion we have provided an efficient recursive
procedure to compute it. In addition, we have developed a method to build up
an explicit expression for any $W_{n}$ in terms of commutators. This is
possible by using similar tools as in the case of the Magnus expansion, namely
by relating products of iterated integrals with the structure of the Hopf
algebra of permutations, and by using special bases of nested commutators. A
sufficient condition for the expansion convergence has also been obtained.

We have presented some application examples of the results about Wilcox
expansion. Firstly, we have shown how to obtain Zassenhaus formula from Wilcox
expansion which, in turn, may be interpreted as its continuous analogue.
Secondly, we point out that Wilcox expansion solves the problem of expanding
the exponential $\exp(A+\epsilon B)$ when $A$ and $B$ are non-commuting
operators. We refer to it as Wilcox–Bellman expansion. Two practical cases on
this respect have been analyzed up to high order. Interestingly, in one of
them the convergence seems not to be uniform. For convenience, the interested
reader can find in [29] a _Mathematica_ code generating general explicit
expressions and recurrences for the Wilcox expansion.

While a full assessment of the Wilcox expansion in comparison with Fer
expansion is not the main purpose of this work, we can still mention here some
of their most distinctive features. Both types of expansions construct the
solution of Eq. (1) as an infinite exponential factorization, but in Wilcox
the exponent of each factor is proportional to successive powers of the
expansion parameter $\lambda$, whereas in Fer each exponent contains an
infinite sum of powers of $\lambda$. This means that, when truncated after a
given number of transformations, say $n$, the Wilcox expansion differs from
the exact solution in the power $\lambda^{n+1}$. In other words, each term
$W_{k}(t)$ in the Wilcox expansion collects the effect of the perturbation at
order $k$. On the other hand, the Fer expansion, when truncated after $n$
transformations, provides a much more accurate approximation. This is true, of
course, if the infinite sums involved in each transformation are exactly
computed, an almost impossible task unless the time dependence of $A(t)$ is
simple enough. By contrast, we have explicit expressions for each exponent
$W_{k}(t)$ in the Wilcox expansion for a generic $A(t)$ and, by using the same
techniques as in the Magnus expansion, we can construct appropriate
approximations to the iterated integrals if necessary. As the examples
collected here and some other studies show [20], Wilcox expansion can provide
accurate results after only a few such transformations.

The first three authors has been funded by Ministerio de Ciencia e Innovación
(Spain) through projects MTM2016-77660-P and PID2019-104927GB-C21 (AEI/FEDER,
UE), and by Universitat Jaume I (grants UJI-B2019-17 and GACUJI/2020/05). The
work of J.A.O. has been partially supported by the Spanish MINECO (grant
numbers AYA2016-81065-C2-2 and PID2019-109592GB-100).

## References

  * [1] Wilcox, R.  Exponential operators and parameter differentiation in quantum physics.  J. Math. Phys. 1967, 8, 962–982.
  * [2] Oteo, J.; Ros, J.  From time-ordered products to Magnus expansion.  J. Math. Phys. 2000, 41, 3268–3277.
  * [3] Hale, J.  Ordinary Differential Equations; Krieger Publishing Company, 1980.
  * [4] Maricq, M.M.  Application of average Hamiltonian theory to the NMR of solids.  Phys. Rev. B 1982, 25, 6662–6632.
  * [5] Casas, F.; Oteo, J.A.; Ros, J.  Floquet theory: exponential perturbative treatment.  J. Phys. A: Math. Gen. 2001, 34, 3379–3388.
  * [6] Mananga, E.; Charpentier, T.  On the Floquet–Magnus expansion: Applications in solid-state nuclear magnetic resonance and physics.  Phys. Rep. 2016, 609, 1–50.
  * [7] Magnus, W.  On the Exponential Solution of Differential Equations for a Linear Operator.  Comm. Pure and Appl. Math. 1954, VII, 649–673.
  * [8] Blanes, S.; Casas, F.; Oteo, J.; Ros, J.  The Magnus expansion and some of its applications.  Phys. Rep. 2009, 470, 151–238.
  * [9] Casas, F.; Chartier, P.; Murua, A.  Continuous changes of variables and the Magnus expansion.  J. Phys. Commun. 2019, 3, 095014.
  * [10] Blanes, S.; Casas, F.; Oteo, J.; Ros, J.  Magnus and Fer expansions for matrix differential equations: the convergence problem.  J. Phys. A: Math. Gen. 1998, 22, 259–268.
  * [11] Fer, F.  Résolution de l’equation matricielle $\dot{U}=p{U}$ par produit infini d’exponentielles matricielles.  Bull. Classe Sci. Acad. Roy. Bel. 1958, 44, 818–829.
  * [12] Bellman, R.  Introduction to Matrix Analysis, 2nd ed.; McGraw-Hill, 1970\.
  * [13] Iserles, A.  Solving Linear Ordinary Differential Equations by Exponentials of Iterated Commutators.  Numer. Math. 1984, 45, 183–199.
  * [14] Klarsfeld, S.; Oteo, J.  Exponential infinite-product representations of the time displacement operator.  J. Phys. A: Math. Gen. 1989, 22, 2687–2694.
  * [15] Mananga, E.  On the Fer expansion: applications in solid-state nuclear magnetic resonance and physics.  Phys. Rep. 2016, 608, 1–41.
  * [16] Casas, F.  Fer’s Factorization as a Symplectic Integrator.  Numer. Math. 1996, 74, 283–303.
  * [17] Zanna, A.  Collocation and Relaxed Collocation for the Fer and the Magnus Expansions.  SIAM J. Numer. Anal. 1999, 36, 1145–1182.
  * [18] Iserles, A.; Munthe-Kaas, H.; Nørsett, S.; Zanna, A.  Lie-Group Methods.  Acta Numerica 2000, 9, 215–365.
  * [19] Huillet, T.; Monin, A.; Salut, G.  Lie algebraic canonical representation in nonlinear control systems.  Math. Systems Theory 1987, 20, 193–213.
  * [20] Zagury, N.; Aragão, A.; Casanova, J.; Solano, E.  Unitary expansion of the time evolution operator.  Phys. Rev. A 2010, 82, 042110.
  * [21] Casas, F.; Murua, A.; Nadinic, M.  Efficient computation of the Zassenhaus formula.  Comput. Phys. Comm. 2012, 183, 2386–2391.
  * [22] Galindo, A.; Pascual, P.  Quantum Mechanics; Springer, 1990.
  * [23] Arnal, A.; Casas, F.; Chiralt, C.  A general formula for the Magnus expansion in terms of iterated integrals of right-nested commutators.  J. Phys. Commun. 2018, 2, 035024.
  * [24] Dragt, A.; Forest, E.  Computation of nonlinear behavior of Hamiltonian systems using Lie algebraic methods.  J. Math. Phys. 1983, 24, 2734–2744.
  * [25] Agrachev, A.; Gamkrelidze, R.  The shuffle product and symmetric groups. In Differential Equations, Dynamical Systems, and Control Science; Elworthy, K.; Everitt, W.; Lee, E., Eds.; Marcel Dekker, 1994; pp. 365–382.
  * [26] Hazewinkel, M.; Gubareni, N.; Kirichenko, V.  Algebras, Rings and Modules. Lie Algebras and Hopf Algebras; AMS, Providence RI, 2010.
  * [27] Malvenuto, C.; Reutenauer, C.  Duality between quasi-symmetric functions and the Solomon descent algebra.  J. Algebra 1995, 177, 967–982.
  * [28] Strichartz, R.S.  The Campbell–Baker–Hausdorff–Dynkin Formula and Solutions of Differential Equations.  J. Funct. Anal. 1987, 72, 320–345.
  * [29] Geometric Integration Research Group. Available online: http://www.gicas.uji.es/Research/Wilcox.html, 2021.
  * [30] Bayen, F.  On the convergence of the Zassenhaus formula.  Lett. Math. Phys. 1979, 3, 161–167.
  * [31] Arnal, A.; Casas, F.; Chiralt, C.  On the structure and convergence of the symmetric Zassenhaus formula.  Comput. Phys. Comm. 2017, 217, 58–65.
  * [32] Moan, P.; Oteo, J.  Convergence of the exponential Lie series.  J. Math. Phys. 2001, 42, 501–508.
  * [33] Casas, F.; Murua, A.  An efficient algorithm for computing the Baker–Campbell–Hausdorff series and some of its applications.  J. Math. Phys. 2009, 50, 033513.
  * [34] Suzuki, M.  On the convergence of exponential operators—the Zassenhaus formula, BCH formula and systematic approximants.  Commun. Math. Phys. 1977, 57, 193–200.
  * [35] Weyrauch, M.; Scholz, D.  Computing the Baker–Campbell–Hausdorff series and the Zassenhaus product.  Comput. Phys. Comm. 2009, 180, 1558–1565.
  * [36] Bellman, R.  Perturbation techniques in Mathematics, Engineering & Physics; Dover, 1972.
  * [37] Arnal, A.; Casas, F.; Chiralt, C.  Exponential perturbative expansions and coordinate transformations.  Math. Comput. Appl. 2020, 25, 50.