#  Splitting techniques for approximating the exponential of commutators

F. Casas1, A. Escorihuela-Tomàs2
1 Departament de Matemàtiques and IMAC, Universitat Jaume I, 12071-Castellón,
Spain
Email<EMAIL_ADDRESS>
2 Departament de Matemàtiques and IMAC, Universitat Jaume I, 12071-Castellón,
Spain
Email<EMAIL_ADDRESS>

###### Abstract

We construct product formulas of orders 3 to 6 approximating the exponential
of a commutator of two arbitrary operators in terms of the exponentials of the
operators involved. The new schemes require a reduced number of exponentials
and thus provide more efficient approximations than other previously published
alternatives, whereas they can be still used as a starting methods of
recursive procedures to increase the order of approximation.

##  I Introduction

The problem of approximating the exponential of the sum of two non-commuting
operators $A$ and $B$ in terms of the exponential of each operator has a long
history, with applications in many areas of physics, chemistry and applied
mathematics [3]. The simplest approximation corresponds of course to the so-
called _Lie–Trotter_ formula
$\chi^{[1]}_{t}=\mathrm{e}^{tA}\,\mathrm{e}^{tB}$, verifying

$\mathrm{e}^{t(A+B)}-\chi^{[1]}_{t}=\mathcal{O}(t^{2})\qquad\mbox{ as
}\;t\rightarrow 0.$ (1)

A higher-order approximation can be achieved by considering a symmetrized
version of (1), also called the _Strang_ (splitting) _scheme_ ,

$S^{[2]}_{t}=\mathrm{e}^{\frac{1}{2}tA}\,\mathrm{e}^{tB}\,\mathrm{e}^{\frac{1}{2}tA}.$
(2)

Specifically, $\mathrm{e}^{t(A+B)}-S^{[2]}_{t}=\mathcal{O}(t^{3})$. In fact,
$S^{[2]}_{t}$ can be used as the starting point of the recursion ($k\geq 2$)
[8, 16]

$\displaystyle S^{[2k]}_{t}$
$\displaystyle=S^{[2k-2]}_{\gamma_{k}t}\,S^{[2k-2]}_{(1-2\gamma_{k})t}\,S^{[2k-2]}_{\gamma_{k}t},$
(3) $\displaystyle\gamma_{k}=\frac{1}{2-2^{1/(2k-1)}},$

leading to approximations of arbitrarily high order $2k$, although the number
of exponentials increases rapidly with the order, and large positive and
negative coefficients are present in (3). The alternative corresponding to the
five maps composition [15]

$\displaystyle S^{[2k]}_{t}$
$\displaystyle=S^{[2k-2]}_{\alpha_{k}t}\,S^{[2k-2]}_{\alpha_{k}t}\,S^{[2k-2]}_{(1-4\alpha_{k})t}\,S^{[2k-2]}_{\alpha_{k}t}\,S^{[2k-2]}_{\alpha_{k}t},$
(4) $\displaystyle\qquad\alpha_{k}=\frac{1}{4-4^{1/(2k-1)}},$

although leading to smaller error terms, involves an even larger number of
exponentials. Approximation schemes of high order involving a reduced number
of exponentials have been designed for the numerical integration of
differential equations, in which case $A$ and $B$ correspond to Lie
derivatives associated to different parts of the defining vector field [2].
These _splitting methods_ are of the form

$S_{t}=\mathrm{e}^{a_{1}tA}\,\mathrm{e}^{b_{1}tB}\,\cdots\,\mathrm{e}^{a_{s}tA}\,\mathrm{e}^{b_{s}tB}$
(5)

and the coefficients $a_{j}$, $b_{j}$ are chosen so that $S_{t}$ provides an
approximation to $\mathrm{e}^{t(A+B)}$ of a given order $r$. This is achieved
by first obtaining and then solving the so-called order conditions, polynomial
equations in $a_{j}$, $b_{j}$ whose number and complexity grow rapidly with
$r$. Thus, for $r=1,\ldots,10$, the number of order conditions is,
respectively, $2$, $1$, $2$, $3$, $6$, $9$, $18$, $30$, $56$ and $99$ [12]. In
any case, the number of exponentials involved in (5) is much smaller than in
compositions (3) and (4). Splitting methods are widely used in the context of
_Geometric Numerical Integration_ , where one is interested in obtaining
numerical solutions still possessing the same qualitative features as the
continuous system they approximate [14, 10, 2]. A more detailed treatment can
be found in [12] and the recent review [3].

Important as it is, not only the object $\exp(t(A+B))$ can be approximated in
this way. In fact, as shown in particular in [11, 7, 6], similar compositions
can be designed to approximate exponentials of commutators and more general
Lie polynomials. The simplest compositions correspond to

$\displaystyle\Psi_{2,1}(t)$
$\displaystyle=\mathrm{e}^{tA}\,\mathrm{e}^{tB}\,\mathrm{e}^{-tA}\,\mathrm{e}^{-tB}\qquad\mbox{
and }\qquad$ (6) $\displaystyle\Psi_{2,2}(t)$
$\displaystyle=\mathrm{e}^{-tB}\,\mathrm{e}^{tA}\,\mathrm{e}^{tB}\,\mathrm{e}^{-tA},$

both verifying that
$\Psi_{2,j}(t)=\mathrm{e}^{t^{2}[A,B]}+\mathcal{O}(t^{3})$. The problem of
constructing product formulas of a given order $r\geq 3$, namely, determining
coefficients $c_{0},c_{1},\ldots$ so that

$\Psi_{r}(t)=\mathrm{e}^{c_{0}tA}\,\mathrm{e}^{c_{1}tB}\,\mathrm{e}^{c_{2}tA}\,\mathrm{e}^{c_{3}tB}\,\cdots=\mathrm{e}^{t^{2}[A,B]}+\mathcal{O}(t^{r+1})$
(7)

has been analyzed in [7, 6], where general recursions are presented. They
allow one, starting from a given approximation of order $m$, to construct a
composition of order $m+1$ or $m+2$, for any $m\geq 2$. In particular,
starting with a basic formula of order 3, different formulas of order 5 are
designed and illustrated in practice. This approach, however, involves an
exceedingly large number of exponentials: the most efficient procedure (called
$\sqrt{10}$-copy recursive formula in [6]) requires $(10N_{m}-4)$ exponentials
to generate a formula of order $m+2$ starting from a formula of (odd) order
$m$ with $N_{m}$ exponentials, and $(10N_{m}-2)$ if the basic scheme is of
even order. Thus, from a composition of order 3 with 6 exponentials, the most
efficient product formula of order 5 involves 56 exponentials.

In this work we address the general problem (7) by applying the same strategy
usually pursued to construct splitting methods to approximate
$\mathrm{e}^{t(A+B)}$: first, we obtain the relevant order conditions required
by such composition and then determine the coefficients by solving these order
conditions. Although different procedures can be used to get the order
conditions (see for instance [3]), here we successively apply the
Baker–Campbell–Hausdorff formula to the product (7) to express $\Psi_{r}(t)$
as the exponential of only one operator (a linear combination of $A$, $B$ and
their nested commutators). Comparison with $\mathrm{e}^{t^{2}[A,B]}$ then
gives the relevant equations to be satisfied by the coefficients $c_{j}$. The
relative efficiency of different solutions (i.e., different sets of
coefficients) is then determined by defining a function related with the
effective error, so that the most efficient scheme is taken as the one
minimizing this function. By following this procedure, we construct schemes of
orders 4, 5 and 6 involving 10, 16 and 26 exponentials, respectively. This
reduction in the number of exponentials contributes a good deal to the
improvement in efficiency of the new methods. Furthermore, by considering
additional exponentials (and therefore additional free parameters) it is
possible to get even more efficient approximations. Of course, the methods we
obtain here can be safely used as basic formulas for the recursive procedures
presented in [6] if one is interested in achieving even higher orders.

As shown in [7] and especially in [6], having efficient product formulas
approximating $\mathrm{e}^{t^{2}[A,B]}$ is essential for some quantum
simulations, in which case the exponentials $\mathrm{e}^{As}$ and
$\mathrm{e}^{Bs^{\prime}}$ are native gates in a quantum circuit. We thus
expect the present results (possibly in combination with the recursions of
[6]) can be applied in practical quantum computations where commutators appear
in a natural way, such as those enumerated there.

The rest of this work is organized as follows. In section II we consider in
detail the case of order 3, thus illustrating the technique we use and the
different types of composition that arise, whereas in section III we analyze
the existing invariances and present new schemes of orders 4, 5 and 6. These
methods are illustrated in section IV on some simple examples. Product
formulas designed to approximate the exponential of other elements in the Lie
algebra generated by $A$ and $B$ are presented in section V. Finally, section
VI contains some concluding remarks.

##  II Third-order schemes

Just as it is the case for order 2, equation (6), two types of compositions
may in principle lead to higher order approximations, namely

$\mathrm{AB}:\qquad\mathrm{e}^{c_{0}tA}\,\mathrm{e}^{c_{1}tB}\,\cdots\,\mathrm{e}^{c_{2n}tA}\,\mathrm{e}^{c_{2n+1}tB},$
(8)

and

$\displaystyle\mathrm{BA}:\qquad\Psi(t)$
$\displaystyle=\prod_{j=0}^{n}\mathrm{e}^{c_{2j}tB}\,\mathrm{e}^{c_{2j+1}tA}$
(9)
$\displaystyle=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\cdots\,\mathrm{e}^{c_{2n}tB}\,\mathrm{e}^{c_{2n+1}tA}.$

For our analysis next we consider compositions of type BA, although, as we
will see, the results are also valid for (8).

If we apply sequentially the Baker–Campbell–Hausdorff formula to (9), we end
up with $\Psi(t)$ expressed as one exponential of a series of operators in
powers of $t$, namely

$\Psi(t)=\exp(Y(t)),\qquad\mbox{ with }\qquad Y(t)=\sum_{j\geq 1}t^{j}Y_{j},$
(10)

and $Y_{j}\in\mathcal{L_{j}}(A,B)$ is an element of the homogeneous subspace
$\mathcal{L}_{j}(A,B)$ of degree $j$ of the graded free Lie algebra
$\mathcal{L}(A,B)=\bigoplus_{j\geq 1}\mathcal{L}_{j}(A,B)$ [13]. One can think
of $\mathcal{L}_{j}$ as the vector subspace formed by linear combinations of
nested commutators involving $j$ operators $A$ and $B$, verifying
$[\mathcal{L}_{n},\mathcal{L}_{m}]\subset\mathcal{L}_{n+m}$. For the
particular basis of $\mathcal{L}_{j}$, $1\leq j\leq 6$, collected in Table 1,
the operator $Y(t)$ reads

$\displaystyle Y(t)=t(w_{1,1}A+w_{1,2}B)+t^{2}w_{2,1}[A,B]$ (11)
$\displaystyle\qquad+t^{3}(w_{3,1}E_{3,1}+w_{3,2}E_{3,2})$
$\displaystyle\qquad+t^{4}\sum_{\ell=1}^{3}w_{4,\ell}E_{4,\ell}+t^{5}\sum_{\ell=1}^{6}w_{5,\ell}E_{5,\ell}$
$\displaystyle\qquad+t^{6}\sum_{\ell=1}^{9}w_{6,\ell}E_{6,\ell}+\mathcal{O}(t^{7}),$

where $w_{n,k}$ are polynomials of homogeneous degree $n$ in the parameters
$c_{j}$. In particular,

$\displaystyle w_{1,1}=\sum_{i=0}^{2n-1}c_{2i+1},\qquad
w_{1,2}=\sum_{i=0}^{2n-1}c_{2i},\qquad$ $\displaystyle
w_{12}=\frac{1}{2}w_{1,1}w_{1,2}-\sum_{i=0}^{n-1}c_{2i}\sum_{j=i}^{n-1}c_{2j+1}.$

From (11), it is clear that an approximation to $\mathrm{e}^{t^{2}[A,B]}$ of
order 3 (i.e., with an error $\mathcal{O}(t^{4})$) is obtained by requiring
that

$w_{1,1}=w_{1,2}=0,\quad w_{2,1}=1,\quad w_{3,1}=w_{3,2}=0.$ (12)

Table 1: Particular basis of ${\mathcal{L}}_{j}$, $1\leq j\leq 6$, taken in this work.  ${\mathcal{L}}_{j}$ |  Basis of ${\mathcal{L}}_{j}$
---|---
${\mathcal{L}}_{1}$ | $E_{1,1}=A$ | $E_{1,2}=B$ |
${\mathcal{L}}_{2}$ | $E_{2,1}=[A,B]$ |  |
${\mathcal{L}}_{3}$ | $E_{3,1}=[A,E_{2,1}]$ | $E_{3,2}=[B,E_{2,1}]$ |
${\mathcal{L}}_{4}$ | $E_{4,1}=[A,E_{3,1}]$ | $E_{4,2}=[B,E_{3,1}]$ | $E_{4,3}=-[B,E_{3,2}]$
${\mathcal{L}}_{5}$ | $E_{5,1}=[A,E_{4,1}]$ | $E_{5,2}=[B,E_{4,1}]$ | $E_{5,3}=[A,E_{4,2}]$
| $E_{5,4}=[B,E_{4,2}]$ | $E_{5,5}=[A,E_{4,3}]$ | $E_{5,6}=[B,E_{4,3}]$
${\mathcal{L}}_{6}$ | $E_{6,1}=[A,E_{5,1}]$ | $E_{6,2}=[B,E_{5,1}]$ | $E_{6,3}=[A,E_{5,2}]$
| $E_{6,4}=[A,E_{5,4}]$ | $E_{6,5}=[B,E_{5,2}]$ | $E_{6,6}=[A,E_{5,5}]$
| $E_{6,7}=[B,E_{5,5}]$ | $E_{6,8}=[A,E_{5,6}]$ | $E_{6,9}=[B,E_{5,6}]$

In consequence, at least 5 exponentials in the composition (9) are required.
It turns out, however, that equations (12) do not have solutions in that case.
It is therefore necessary to include an additional exponential, i.e., to
consider the scheme

$\Psi_{3}(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\mathrm{e}^{c_{2}tB}\,\mathrm{e}^{c_{3}tA}\,\mathrm{e}^{c_{4}tB}\,\mathrm{e}^{c_{5}tA},$
(13)

in which case one has a free parameter, say $c_{5}$. The general solution is
then given by ($c_{5}\neq 0$)

$\displaystyle
c_{0}=\frac{1\mp\sqrt{5}}{2c_{5}},\;c_{1}=\frac{c_{5}(-1\pm\sqrt{5})}{2},\;c_{2}=\frac{1}{c_{5}},\;$
(14) $\displaystyle
c_{3}=\frac{c_{5}(-1\mp\sqrt{5})}{2},\;c_{4}=\frac{-3\pm\sqrt{5}}{2c_{5}}.$

It is a common practice when designing splitting methods to choose the free
parameter so as to minimize some objective function typically related with the
leading error term in the asymptotic expansion of $Y(t)$ in (11), in this case
$\sum_{\ell=1}^{3}w_{4,\ell}E_{4,\ell}$. In general, given a composition (9)
involving $s$ exponentials that approximates $\mathrm{e}^{t^{2}[A,B]}$ up to
order $r$, we define its effective error as

$\mathcal{E}^{(r+1)}=s\
\left(\sqrt{\sum_{\ell=1}^{\dim{\mathcal{L}}_{r+1}}w_{r+1,\ell}E_{r+1,\ell}}\right)^{1/r},$
(15)

to take also into account the computational effort required by the scheme [3].
A method of order $r$ is typically more efficient in practice than another one
when it leads to a smaller $\mathcal{E}^{(r+1)}$, so it makes sense to take
the effective error as the objective function to minimize. One should be
aware, however, that the value of $\mathcal{E}^{(r+1)}$ may differ depending
on the basis of ${\mathcal{L}}$ one is using.

In the particular case of the 3rd-order scheme (13)-(14), the minimum value of
$\mathcal{E}^{(4)}$ is obtained with the solution

$\displaystyle c_{0}=\pm\sqrt{\frac{2}{\sqrt{5}+1}},\quad
c_{1}=\mp\sqrt{\sqrt{5}-2},\quad c_{2}=\mp\sqrt{\frac{2}{\sqrt{5}-1}},$ (16)
$\displaystyle c_{3}=-c_{2},\qquad c_{4}=-c_{1},\qquad c_{5}=-c_{0}.$

Notice the peculiar arrangement in the distribution of the coefficients in
this case:

$(c_{0},c_{1},c_{2},-c_{2},-c_{1},-c_{0}).$ (17)

In fact, the same minimum is also obtained by the sequence ($i=\sqrt{-1}$)

$\displaystyle(\tilde{c}_{0},\tilde{c}_{1},\tilde{c}_{2},\tilde{c}_{2},\tilde{c}_{1},\tilde{c}_{0}),\qquad\qquad$
(18)
$\displaystyle\tilde{c}_{0}=ic_{0},\qquad\tilde{c}_{1}=-ic_{1},\qquad\tilde{c}_{2}=ic_{2},$

with $c_{j}$ given in (17), although in this case all the coefficients are
pure imaginary. Sequences such as (17) and (18) turn out to be very convenient
when designing higher order methods.

##  III Higher-order compositions

Before proceeding to the actual construction of higher-order schemes, it is
worth to briefly review some of the invariances that a composition of type (9)
has when approximating $\mathrm{e}^{t^{2}[A,B]}$. This, on the one hand,
provides more insight into the solutions obtained at order 3 and, on the other
hand, may help in our exploration.

###  III.1 Invariances

The first and most obvious invariance takes place when one reverses the sign
of $t$ in (9). If $\Psi_{r}(t)$ is a method of order $r$, i.e.,

$\displaystyle\Psi_{r}(t)$
$\displaystyle=\prod_{j=0}^{n}\mathrm{e}^{c_{2j}tB}\,\mathrm{e}^{c_{2j+1}tA}$
$\displaystyle=\exp\big{(}t^{2}[A,B]+D_{r+1}t^{r+1}+\mathcal{O}(t^{r+2})\big{)},$

then

$\displaystyle\Psi_{r}(-t)$
$\displaystyle=\prod_{j=0}^{n}\mathrm{e}^{-c_{2j}tB}\,\mathrm{e}^{-c_{2j+1}tA}$
$\displaystyle=\exp\big{(}t^{2}[A,B]+(-1)^{r+1}D_{r+1}t^{r+1}+\mathcal{O}(t^{r+2})\big{)}.$

In consequence, and according with the definition of $\mathcal{E}^{(r+1)}$,
both $\Psi_{r}(t)$ and $\Psi_{r}(-t)$ lead to the same effective error. This
is in agreement with the two families of solutions (17) differing by a sign.
In other words, when we propose a method with some specific coefficients
$c_{j}$, the scheme obtained by considering the same composition with
coefficients $-c_{j}$ has the same effective error.

Consider now the effect of the transformation

$t\longmapsto i\,\tilde{t},\qquad A\longmapsto-\tilde{A}$

on composition (9). Then

$\displaystyle\Psi_{r}(t)\longmapsto\widetilde{\Psi}_{r}(\tilde{t})$
$\displaystyle=\mathrm{e}^{ic_{0}\tilde{t}B}\,\mathrm{e}^{-ic_{1}\tilde{t}\tilde{A}}\,\cdots\,\mathrm{e}^{ic_{2n}\tilde{t}B}\,\mathrm{e}^{-ic_{2n+1}\tilde{t}\tilde{A}}$
$\displaystyle=\exp\big{(}\tilde{t}^{2}[\tilde{A},B]+(i)^{r+1}\widetilde{D}_{r+1}\,\tilde{t}^{r+1}+\mathcal{O}(\tilde{t}^{r+2})\big{)},$

where $\widetilde{D}_{r+1}$ is the same linear combination of commutators as
$D_{r+1}$, possibly with some coefficients with the reverse sign, thus leading
again to a scheme with the same effective error. This fact accounts for the
second family of solutions (18) found when $r=3$.

Finally, consider the transformation

$A\longmapsto\hat{B},\qquad B\longmapsto-\hat{A}.$ (19)

Then, clearly

$\displaystyle\Psi_{r}(t)\longmapsto\hat{\Psi}_{r}(t)$
$\displaystyle=\mathrm{e}^{-c_{0}t\hat{A}}\,\mathrm{e}^{c_{1}t\hat{B}}\,\cdots\,\mathrm{e}^{-c_{2n}t\hat{A}}\,\mathrm{e}^{c_{2n+1}t\hat{B}}$
$\displaystyle=\exp\big{(}t^{2}[\hat{A},\hat{B}]+\widehat{D}_{r+1}\,t^{r+1}+\mathcal{O}(t^{r+2})\big{)}.$

In other words, a BA method is transformed into an AB scheme, although
possibly with a different value of $\mathcal{E}^{(r+1)}$, since the map (19)
may modify the basis. In the particular case of the basis collected in Table
1, one has $\hat{E}_{6,4}=E_{6,4}+\frac{1}{3}(E_{6,5}+E_{6,6})$. All other
terms $\hat{E}_{j,\ell}$ up to $j=6$ only differ from $E_{j,\ell}$ in a sign.
We can therefore safely consider schemes of type BA, at least up to this order
of approximation. Notice how $\Psi_{2,1}(t)$ can be obtained from
$\Psi_{2,2}(t)$ in (6) in this way.

###  III.2 ‘Counter-palindromic’ patterns

We have seen how the sequences (17) and (18) lead, when $r=3$, to methods with
the minimum value of the effective error, such as is defined in (15). Scheme
$\Psi_{2,2}(t)$ in (6) actually follows the same pattern:
$(c_{0},c_{1},c_{1},c_{0})$. It seems then reasonable to consider similar
sequences of coefficients when the order of approximation $r>3$, i.e., either

$(c_{0},c_{1},c_{2},\ldots,c_{m-1},c_{m},c_{m},c_{m-1},\ldots,c_{2},c_{1},c_{0})$
(20)

or

$(c_{0},c_{1},c_{2},\ldots,c_{m-1},c_{m},-c_{m},-c_{m-1},\ldots,-c_{2},-c_{1},-c_{0}).$
(21)

Schemes of this class will be called (positive or negative) _counter-
palindromic_ compositions, to emphasize the fact that the coefficients
multiply a different operator ($A$ or $B$) the second time they appear in the
sequence (with the same or opposite sign, respectively). Sequences of this
type have been previously considered in [9] when designing symplectic
integrators for the (space discretized) Schrödinger equation.

Moreover, it turns out that sequences (20) and (21) actually lead to a
significant reduction in the number of order conditions, since they already
satisfy by construction many of them. Specifically, the coefficients in the
expansion (11) of the operator $Y(t)$ associated to the scheme verify the
identities (up to order 6):

$\begin{array}[]{lllc}w_{1,1}=\pm w_{1,2}&w_{3,1}=\mp
w_{3,2}&w_{4,1}=-w_{4,3}&\\\ w_{5,1}=\pm w_{5,6}&w_{5,2}=\pm
w_{5,5}&w_{5,3}=\pm w_{5,4}&\\\
w_{6,1}=-w_{6,9}&w_{6,2}=-w_{6,8}&w_{6,3}=-w_{6,7}&\\\
w_{6,4}=3(w_{6,5}+w_{6,6}).&&&\end{array}$

Here the top sign corresponds to the pattern (20) and the bottom sign to (21).

This can be seen as follows. Suppose we have a composition of type (20),

$\Psi(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\cdots\,\mathrm{e}^{c_{m-1}tB}\,\mathrm{e}^{c_{m}tA}\,\mathrm{e}^{c_{m}tB}\,\mathrm{e}^{c_{m-1}tA}\,\cdots\,\mathrm{e}^{c_{1}tB}\,\mathrm{e}^{c_{0}tA},$
(22)

which can be formally expressed as $\Psi(t)=\exp(Y(t))$, with $Y(t)$ given by
(11), and interchange $A$ and $B$ in (22). This results in the new method

$\widetilde{\Psi}(t)=\mathrm{e}^{c_{0}tA}\,\mathrm{e}^{c_{1}tB}\,\cdots\,\mathrm{e}^{c_{m-1}tA}\,\mathrm{e}^{c_{m}tB}\,\mathrm{e}^{c_{m}tA}\,\mathrm{e}^{c_{m-1}tB}\,\cdots\,\mathrm{e}^{c_{1}tA}\,\mathrm{e}^{c_{0}tB},$

whose associated operator was $\widetilde{Y}(t)$ in
$\widetilde{\Psi}(t)=\exp(\widetilde{Y}(t))$ reads

$\displaystyle\widetilde{Y}(t)=t(w_{1,1}B+w_{1,2}A)-t^{2}w_{2,1}[A,B]$
$\displaystyle\qquad+t^{3}(-w_{3,1}E_{3,2}+w_{3,2}E_{3,1})$
$\displaystyle\qquad+t^{4}(w_{4,1}E_{4,3}-w_{4,2}E_{4,2}+w_{4,3}E_{4,1})+\mathcal{O}(t^{5})$

Now, the adjoint of $\widetilde{\Psi}(t)$, defined as
$\big{(}\widetilde{\Psi}(t)\big{)}^{*}=\big{(}\widetilde{\Psi}(-t)\big{)}^{-1}$
[10], is precisely $\Psi(t)$, so that

$-\widetilde{Y}(-t)=Y(t),$

whence, by comparing the corresponding expansions, one arrives at
$w_{1,1}=w_{1,2}$, $w_{3,1}=-w_{3,2}$, $w_{4,1}=-w_{4,3}$, etc. A similar
argument is also valid for compositions (21) by considering the interchange
$A\longleftrightarrow-B$.

As a result, designing counter-palindromic methods of order $r\geq 4$ requires
solving a smaller number of polynomial equations than with more general
compositions. This is clearly visible in Table 2, where we collect the number
of order conditions required to achieve order $3\leq r\leq 6$ for a general
composition (9) (first row) and the corresponding number for a counter-
palindromic composition (second row). For clarity, we have also included the
minimum number $s$ of exponentials required in each case (in bold).

Table 2: Number of independent order conditions to achieve order $r$ for a general composition (9) (first row) and for counter-palindromic compositions (second row). The number $s$ refers to the number of exponentials involved in the composition. Order $r$ |  3 |  4 |  5 |  6
---|---|---|---|---
General |  5 (6) |  8 ($\geq$ 8) |  14 ($\geq$ 14) |  23 ($\geq$ 23)
CP ($s$) |  3 (6) |  5 (10) |  8 (16) |  13 (26)

Although the number of exponentials increases with respect to the general case
(e.g., 16 vs. 14 for $r=5$), the number of equations to solve (and therefore
the complexity of the problem) reduces considerably (14 vs. 8 equations for
$r=5$). Moreover, although the number of order conditions sets up the minimum
number of exponentials involved in a general composition, this by itself does
not guarantee the existence of real solutions. This fact has been illustrated
for $r=3$ in section II. On the other hand, we recall that the most efficient
5th-order composition built in [6] involves 56 exponentials.

###  III.3 New schemes

We have explored both types of sequences (20) and (21) and solved the required
order conditions to achieve up to order 6 with the minimum number of
exponentials in each case. This task is simplified by considering schemes with
patterns (20) and (21). The coefficients of the most efficient methods (i.e.,
with the smallest value of $\mathcal{E}^{(r+1)}$) are collected in Table 3.
The resulting method is denoted as $\mathcal{PCP}_{s}^{[r]}$ (respectively,
$\mathcal{NCP}_{s}^{[r]}$) if it corresponds to a composition of class (20)
(respect., (21)) of order $r$ involving $s$ exponentials. For clarity, the
methods read, respectively

$\mathcal{PCP}:\begin{array}[]{l}\Psi(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\cdots\,\mathrm{e}^{c_{m-1}tB}\,\mathrm{e}^{c_{m}tA}\\\
\qquad\quad\times\mathrm{e}^{c_{m}tB}\,\mathrm{e}^{c_{m-1}tA}\,\cdots\,\mathrm{e}^{c_{1}tB}\,\mathrm{e}^{c_{0}tA}\quad\mbox{
($m$ odd) }\\\
\Psi(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\cdots\,\mathrm{e}^{c_{m-1}tA}\,\mathrm{e}^{c_{m}tB}\\\
\qquad\quad\times\mathrm{e}^{c_{m}tA}\,\mathrm{e}^{c_{m-1}tB}\,\cdots\,\mathrm{e}^{c_{1}tB}\,\mathrm{e}^{c_{0}tA}\quad\mbox{
($m$ even) }\end{array}$ (23)

and

$\mathcal{NCP}:\begin{array}[]{l}\Psi(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\cdots\,\mathrm{e}^{c_{m-1}tB}\,\mathrm{e}^{c_{m}tA}\\\
\qquad\quad\times\mathrm{e}^{-c_{m}tB}\,\mathrm{e}^{-c_{m-1}tA}\,\cdots\,\mathrm{e}^{-c_{1}tB}\,\mathrm{e}^{-c_{0}tA}\quad\mbox{
($m$ odd) }\\\
\Psi(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\cdots\,\mathrm{e}^{c_{m-1}tA}\,\mathrm{e}^{c_{m}tB}\\\
\qquad\quad\times\mathrm{e}^{-c_{m}tA}\,\mathrm{e}^{-c_{m-1}tB}\,\cdots\,\mathrm{e}^{-c_{1}tB}\,\mathrm{e}^{-c_{0}tA}\quad\mbox{
($m$ even) }\end{array}$ (24)

Table 3: Coefficients of the most efficient counter-palindromic schemes of
order $3\leq r\leq 6$ with the minimum number of exponentials.

$\mathcal{NCP}_{6}^{[3]},\qquad\qquad\mathcal{E}^{(4)}/s\approx 0.473$
$c_{0}=c_{1}-c_{2}$ $c_{1}=-\sqrt{\sqrt{5}-2}$
$c_{2}=-\sqrt{\frac{2}{\sqrt{5}-1}}$
$\mathcal{NCP}_{10}^{[4]},\qquad\qquad\mathcal{E}^{(5)}/s\approx 0.606$
$c_{0}=\sum_{j=1}^{4}(-1)^{j+1}c_{j}$ $c_{1}=0.4920434066428167763156$
$c_{2}=-1.569846260451462851779$ $c_{3}=-0.0340560371300231615989$
$c_{4}=3.007307207357765662262$
$\mathcal{PCP}_{16}^{[5]},\qquad\qquad\mathcal{E}^{(6)}/s\approx 0.505$
$c_{0}=-\sum_{j=1}^{7}c_{j}$ $c_{1}=0.2969175443796203417835$
$c_{2}=1.418243492034305431995$ $c_{3}=0.4347212029859471608694$
$c_{4}=-0.127142127469064995044$ $c_{5}=-2.014276365712093993010$
$c_{6}=0.8493401946712687892513$ $c_{7}=-0.305642216160471071886$
$\mathcal{PCP}_{26}^{[6]},\qquad\qquad\mathcal{E}^{(7)}/s\approx 0.447$
$c_{0}=-\sum_{j=1}^{12}c_{j}$ $c_{1}=0.2464427486685065253599$
$c_{2}=0.437855533639627516106$ $c_{3}=-0.6290554972825559401392$
$c_{4}=-1.160402744300525331934$ $c_{5}=-0.5248160600039844378749$
$c_{6}=-0.2264322765760404736976$ $c_{7}=0.1165418804073705040233$
$c_{8}=0.4687839445292851414849$ $c_{9}=1.983312306755703005101$
$c_{10}=-0.9894918460835968618662$ $c_{11}=0.6722571007458945095097$
$c_{12}=-0.2387711966553848135336$

We should remark that, due to the symmetries and invariances explored in
subsection III.1, there are more compositions than those collected in Table 3
leading to the same efficiency.

A useful strategy to improve the efficiency of splitting methods when
approximating $\mathrm{e}^{t(A+B)}$ consists in including additional
exponentials in the composition, and therefore additional parameters, which
are then fixed in the optimization process of the objective function [3]. In
fact, very often the most efficient method does not necessarily correspond to
the scheme requiring the minimum number of exponentials: the extra cost can be
compensated by the additional accuracy obtained, although solving the
polynomial equations with additional exponentials and free parameters is not a
trivial task. By pursuing the same strategy in this context we have obtained
the schemes of order 4 and 5 collected in Table 4. Notice how the effective
error is reduced: the extra cost resulting from an increased number of
exponentials is compensated by the higher accuracy achieved.

Table 4: Coefficients of optimized counter-palindromic schemes of order 4 and
5 involving 12 and 18 exponentials, respectively.

$\mathcal{PCP}_{12}^{[4]},\qquad\qquad\mathcal{E}^{(5)}/s\approx 0.455$
$c_{0}=-\sum_{j=1}^{5}c_{j}$ $c_{1}=0.3263285743794757829237$
$c_{2}=-1.564170317916158642032$ $c_{3}=-0.0234725141740210902965$
$c_{4}=2.920816850699232751348$ $c_{5}=-0.8045459762846959202889$
$\mathcal{NCP}_{18}^{[5]},\qquad\qquad\mathcal{E}^{(6)}/s\approx 0.395$
$c_{0}=\sum_{j=1}^{8}(-1)^{j+1}c_{j}$ $c_{1}=-0.6410115692148225407946$
$c_{2}=0.3165189600901244909982$ $c_{3}=0.2075766074841999769730$
$c_{4}=-1.042459743800714071012$ $c_{5}=1.027769699504593533740$
$c_{6}=1.290831433928573680468$ $c_{7}=0.7061407649397449413288$
$c_{8}=0.253358191085494126186$

##  IV Numerical examples

In this section we illustrate the previous methods for some simple examples.
The first one has also been used as a test bench in [6] and consists in
approximating the exponential of the commutator of Pauli matrices, namely
$\mathrm{e}^{[-i\sigma_{x},-i\sigma_{z}]}$, by products of elementary
exponentials of $A=-i\sigma_{x}$ and $B=-i\sigma_{z}$. Specifically, we
compute

$\big{\|}\Psi(1/\sqrt{n})^{n}-\mathrm{e}^{[A,B]}\big{\|}_{2}$ (25)

for several values of $n$ and plot its error as a function of the total number
of exponentials $ns$ required by each method $\Psi(t)$ to take into account
its overall computational cost. In this way we get the diagram of Figure 1
(left).

Figure 1: Error committed in the approximation of $\mathrm{e}^{[A,B]}$ by
different product formulas vs. the number of elementary exponentials when
$A=-i\sigma_{x}$, $B=-i\sigma_{z}$ (left) and when $A$ and $B$ are $16\times
16$ real matrices with random elements (right). Methods of Table 3 are
compared with schemes designed in [6].

Here $\mathcal{S}_{2}$, $\mathcal{S}_{3}$ and $\mathcal{G}_{5}$ correspond to
schemes of order 2, 3 and 5 proposed in [6], involving 4, 6 and 56
exponentials, respectively. Actually, $\mathcal{S}_{2}$ is the approximation
$\Psi_{2,1}$ of (6), whereas $\mathcal{G}_{5}$ is obtained by applying the
most efficient ($\sqrt{10}$-copy) recursive formula presented in [6] to the
basic scheme $\mathcal{S}_{3}$. We also show in the diagram the results
achieved by the methods of Table 3. Notice that, whereas the improvement
achieved by the new 3-order approximation $\mathcal{NCP}_{6}^{[3]}$ is only
marginal, this is not the case for the new schemes of order 4 and 5,
essentially due to the much reduced number of exponentials involved. The
improvement is still more remarkable in the case of the 6th-order method
$\mathcal{PCP}_{26}^{[6]}$, which in fact provides the most efficient
approximation. This method should be compared with the 6th-order scheme
labelled $\mathcal{G}_{6}$, which is obtained by applying the recursion of [6]
to the basic scheme $\mathcal{NCP}_{10}^{[4]}$ (involving a total of 98
exponentials). In all cases, the leftmost point along the line corresponds to
$n=1$.

To discard the possible effect of the low dimensionality ($2\times 2$) and the
particular structure (skew-Hermitian) of the matrices involved, we repeat the
same experiment, but this time by constructing two $16\times 16$ real matrices
$\tilde{A}$, $\tilde{B}$ with elements randomly generated from a normal
distribution and then taking $A=\tilde{A}/\|\tilde{A}\|$,
$B=\tilde{B}/\|\tilde{B}\|$. As before, we plot in Figure 1 (right) the
efficiency diagram obtained by the methods of Table 3 in comparison with the
most efficient approximation designed in [6]. We observe the same behavior as
for Pauli matrices, with the new 6th-order method $\mathcal{PCP}_{26}^{[6]}$
showing the best performance.

Figure 2 illustrates, on the same examples, the behavior of the optimized
methods of order 4 and 5 whose coefficients are collected in Table 4. Notice
that, even if the computational cost per step increases, they still show a
better overall efficiency in both cases.

Figure 2: Error in the approximation of $\mathrm{e}^{[A,B]}$ obtained by
product formulas of Table 4 when $A=-i\sigma_{x}$, $B=-i\sigma_{z}$ (left) and
when $A$ and $B$ are $16\times 16$ real matrices with random elements (right).

We can also consider values $t>1$ to check how the different schemes perform
in the long run. To this end, we take $t=10$ and compute accordingly

$\big{\|}\Psi(10/\sqrt{n})^{n}-\mathrm{e}^{10[A,B]}\big{\|}_{2}.$

The corresponding efficiency diagrams for Pauli (left) and random $16\times
16$ matrices (right) are shown in Figure 3.

Figure 3: Same as Figure 1, but now when approximating $\mathrm{e}^{t[A,B]}$
with $t=10$.

In our next experiment, we determine the number of elementary exponentials
(gates) $ns$ required for each scheme $\Psi(t)$ so that the error for
different values of $x$ in $\mathrm{e}^{x^{2}[A,B]}$, namely

$\big{\|}\Psi(x/\sqrt{n})^{n}-\mathrm{e}^{x^{2}[A,B]}\big{\|}_{2}$ (26)

is smaller than a prescribed tolerance. We choose again $A=-i\sigma_{x}$,
$B=-i\sigma_{z}$, and $x\in[0.1,0.9]$. The corresponding results are depicted
in Figure 4 for a tolerance $10^{-4}$ (left) and $10^{-7}$ (right). The
improvement with respect to the most efficient scheme presented in [6] is
quite remarkable, even for small values of $x$.

Figure 4: Number of elementary exponentials (gates) required by each product
formula to approximate $\mathrm{e}^{x^{2}[-i\sigma_{x},-i\sigma_{z}]}$ with
error smaller than $10^{-4}$ (left) and $10^{-7}$ (right).

##  V Extensions

Not only the exponential of a commutator can be approximated by products of
exponentials of the elementary operators involved, but in fact the same
technique can be applied to any element of the free Lie algebra generated by
$A$ and $B$. In this section we illustrate the procedure in two particular
cases, already considered previously in the literature.

###  V.1 Approximating the exponential of $F=t(A+B)+t^{2}R[A,B]$

Suppose one is interested in computing an approximation to the exponential of
both the sum and commutator of $A$, $B$, namely, $F=t(A+B)+t^{2}R[A,B]$, for
an arbitrary $R\in\mathbb{R}$. A particular application arises in the
simulation of counter-adiabatic time evolution of two-qubit Hamiltonians [6].

Clearly, the coefficients of an approximation of the form (9) of order $r$,
namely,

$\mathrm{e}^{c_{0}tRB}\,\mathrm{e}^{c_{1}tRA}\,\cdots\,\mathrm{e}^{c_{2n}tRB}\,\mathrm{e}^{c_{2n+1}tRA}=\mathrm{e}^{F}+\mathcal{O}(r+1)$

have to verify

$\displaystyle\qquad w_{1,1}=w_{1,2}=1,\qquad w_{2,1}=R^{2},$ $\displaystyle
w_{j,\ell}=0,\qquad 3\leq j\leq r,\;\ell=1,\ldots,\dim({\mathcal{L}}_{j}).$

Thus, we get an approximation of order 2 with only 3 exponentials:

$\displaystyle\qquad\Phi_{3}^{[2]}(t)=\mathrm{e}^{c_{0}tRB}\,\mathrm{e}^{c_{1}tRA}\,\mathrm{e}^{c_{2}tRB},\qquad$
$\displaystyle c_{0}=-\frac{2R^{2}-1}{2R},\quad c_{1}=\frac{1}{R},\quad
c_{2}=\frac{2R^{2}+1}{2R},$

whereas the inclusion of a free parameter $c_{3}$ allows us to reduce the
effective error:

$\Phi_{4}^{[2]}(t)=\mathrm{e}^{c_{0}tRB}\,\mathrm{e}^{c_{1}tRA}\,\mathrm{e}^{c_{2}tRB}\,\mathrm{e}^{c_{3}tRA},$

with

$\displaystyle c_{0}=\frac{2R^{2}+2Rc_{3}-1}{2R\left(Rc_{3}-1\right)},\quad
c_{1}=-\frac{Rc_{3}-1}{R},$ $\displaystyle
c_{2}=-\frac{2R^{2}+1}{2R\left(Rc_{3}-1\right)}.$

Moreover, in this case we recover the scheme $\Psi_{2,2}(t)$ of (6) in the
limit $R\rightarrow\infty$ with $c_{3}=-1$.

By following the same strategy, the following 3rd-order scheme can be
obtained:

$\Phi_{5}^{[3]}(t)=\mathrm{e}^{c_{0}tRB}\,\mathrm{e}^{c_{1}tRA}\,\mathrm{e}^{c_{2}tRB}\,\mathrm{e}^{c_{3}tRA}\,\mathrm{e}^{c_{4}tRB},$

with coefficients

$\displaystyle
c_{0}=\frac{1}{R}\left(3R^{4}-R^{2}\mp\Delta+\frac{1}{4}\right),\quad
c_{1}=\frac{6R^{4}\pm 2\Delta-\frac{1}{2}}{R(12R^{4}-1)},\quad$ $\displaystyle
c_{2}=\frac{1}{2R}-6R^{3},\quad c_{3}=\frac{6R^{4}\mp
2\Delta-\frac{1}{2}}{R(12R^{4}-1)},$ $\displaystyle
c_{4}=\frac{1}{R}\left(3R^{4}-R^{2}\pm\Delta+\frac{1}{4}\right),$
$\displaystyle\Delta=\frac{\sqrt{3(12R^{4}-1)(36R^{4}+1)}}{12}.$

It is worth noticing that, whereas 6 exponentials are necessary to approximate
$\mathrm{e}^{t^{2}[A,B]}$ up to order 3 (cf. Section III), with only 5 we are
able to reproduce $\mathrm{e}^{F}$ in general, in contrast with the treatment
in [6].

###  V.2 Approximating $\mathrm{e}^{t^{3}[A,[A,B]]}$

In some problems arising in quantum computing it is desirable to suppress the
effect of undesirable interactions in coupled spin systems by using composite
pulse sequences to generate an effective Hamiltonian for which the nested
commutator $[A,[A,B]]$ is the dominant term [4]. This can itself be
approximated by a product of exponentials of $A$ and $B$. Thus, a
straightforward calculation shows that

$\mathrm{e}^{tA}\,\mathrm{e}^{tB}\,\mathrm{e}^{-tA}\,\mathrm{e}^{-tB}\,\mathrm{e}^{-tA}\,\mathrm{e}^{tB}\,\mathrm{e}^{tA}\,\mathrm{e}^{-tB}=\mathrm{e}^{t^{3}[A,[A,B]]}+\mathcal{O}(t^{4}),$
(27)

whereas, if the identity $[A,[B,[B,A]]]=0$ holds, then error of this
approximation is of order $\mathcal{O}(t^{5})$ [7]. By applying the same
strategy as in section III it is possible to construct a left-right
palindromic composition of order 4 involving only 9 exponentials,

$\Psi(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\mathrm{e}^{c_{2}tB}\,\mathrm{e}^{c_{3}tA}\,\mathrm{e}^{c_{4}tB}\,\mathrm{e}^{c_{3}tA}\,\mathrm{e}^{c_{2}tB}\,\mathrm{e}^{c_{1}tA}\,\mathrm{e}^{c_{0}tB}.$
(28)

In fact, this composition has essentially the same computational cost as (27)
when it is used iteratively (since the last exponential can be concatenated
with the first one at the next step). The analysis shows that there is a free
parameter in (28), so that the sequence of coefficients reads

$(c_{0},c_{1},c_{2},c_{3},c_{4})=\left(-\frac{c_{2}}{2},\pm\frac{1}{\sqrt{c_{2}}},c_{2},\mp\frac{1}{\sqrt{c_{2}}},-c_{2}\right)$

and the minimum effective error is achieved with
$c_{2}=\left((\sqrt{1346}-36)/25\right)^{1/3}$. By following the same
approach, an approximation of order 6 can be obtained with 21 exponentials. On
the other hand, the composition

$\Psi(t)=\mathrm{e}^{c_{0}tB}\,\mathrm{e}^{c_{1}tA}\,\mathrm{e}^{c_{2}tB}\,\mathrm{e}^{c_{3}tA}\,\mathrm{e}^{c_{4}tB},$

involving only 5 exponentials verifies

$\Psi(t)=\exp\big{(}t(A+B)+t^{2}[A,B]+t^{3}[A,[A,B]]\big{)}+\mathcal{O}(t^{4})$

when the coefficients are given by

$\displaystyle c_{0}=-\frac{3}{4}+\frac{\alpha}{4},\quad
c_{1}=\frac{1}{2}+\frac{\alpha}{2},\quad c_{2}=\frac{1}{2},$
$\displaystyle\qquad c_{3}=\frac{1}{2}-\frac{\alpha}{2},\quad
c_{4}=\frac{5}{4}-\frac{\alpha}{4}$

and $\alpha=\sqrt{47/3}$.

##  VI Concluding remarks

In this work we have obtained efficient approximations for the exponential of
the commutator $\mathrm{e}^{t^{2}[A,B]}$ of two arbitrary operators $A$ and
$B$ as compositions of exponentials of the operators involved. Previous
analyses of this same problem show that it is in fact possible to get
arbitrarily high-order approximations by considering recursive formulas
starting from a basic low-order scheme, in a similar way as the well known
triple (3) and quintuple jump (4) formulas do for approximating
$\mathrm{e}^{t(A+B)}$. Although these recursive formulas allow one to get
explicit upper bounds for the error committed and on the number of
exponentials required, this number is exceedingly large, even for moderate
orders of approximation, and so their application to solve problems in
practice is questionable at least. This is the case in particular in quantum
computing, where the action of $\mathrm{e}^{tA}$ and $\mathrm{e}^{tB}$ is
simulated by a quantum circuit. In that setting it is clearly relevant to
reduce to a minimum the number of exponentials required to achieve a certain
accuracy in the approximation.

Here, by applying the same techniques as in the construction of splitting
methods for the time integration of differential equations, and considering
particular patterns in the distribution of the exponentials, we have obtained
product formulas of orders 4, 5 and 6 with a much reduced computational cost
(only two more exponentials than the theoretical minimum). The numerical
examples we gather clearly show how this reduction greatly contributes to the
improvement in the efficiency. The new methods could then be applied in
quantum simulations where it is essential to have a number of gates as small
as possible whereas still having good accuracy. Other possible applications
include the computation of the successive terms in the Zassenhaus formula [5]
and its continuous symmetric analogue [1].

Although we have limited our analysis up to order 6, it is clear that similar
calculations can be carried out for higher orders if necessary. In any case,
the new methods can be used as basic schemes for the general recursive
procedures developed in [7, 6].

### Acknowledgements

This work has been funded by Ministerio de Ciencia e Innovación (Spain)
through project PID2022-136585NB-C21, MCIN/AEI/10.13039/501100011033/FEDER,
UE, and also by Generalitat Valenciana (Spain) through project
CIAICO/2021/180. The authors wish to thank S. Blanes for his very useful
comments.

## References

  * [1] P. Bader, A. Iserles, K. Kropielnicka, and P. Singh, Effective approximation for the semiclassical Schrödinger equation, Found. Comput. Math., 14 (2014), pp. 689–720.
  * [2] S. Blanes and F. Casas, A Concise Introduction to Geometric Numerical Integration, CRC Press, 2016.
  * [3] S. Blanes, F. Casas, and A. Murua, Splitting methods for differential equations, tech. rep., arXiv:2401.01722, 2024.
  * [4] T. W. Borneman, C. E. Granade, and D. G. Cory, Parallel information transfer in a multinode quantum information processor, Phys. Rev. Lett., 108 (2012), p. 140502.
  * [5] F. Casas, A. Murua, and M. Nadinic, Efficient computation of the Zassenhaus formula, Comput. Phys. Comm., 183 (2012), pp. 2386–2391.
  * [6] Y.-A. Chen, A. M. Childs, M. Hafezi, Z. Jiang, H. Kim, and Y. Xu, Efficient product formulas for commutators and applications to quantum simulations, Phys. Rev. Res., 4 (2022), p. 013191.
  * [7] A. M. Childs and N. Wiebe, Product formulas for exponentials of commutators, J. Math. Phys., 54 (2013), p. 062202.
  * [8] M. Creutz and A. Gocksch, Higher-order hybrid Monte Carlo algorithms, Phys. Rev. Lett., 63 (1989), pp. 9–12.
  * [9] S. Gray and D. E. Manolopoulos, Symplectic integrators tailored to the time-dependent Schrödinger equation, J. Chem. Phys., 104 (1996), pp. 7099–7112.
  * [10] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, Springer-Verlag, Second ed., 2006.
  * [11] F. Jean and P.-V. Koseleff, Elementary approximation of exponentials of Lie polynomials, in Applied Algebra, Algebraic Algorithms and Error-Correcting Codes, T. Mora and H. Mattson, eds., Springer, 1997, pp. 174–188.
  * [12] R. I. McLachlan and R. Quispel, Splitting methods, Acta Numerica, 11 (2002), pp. 341–434.
  * [13] H. Munthe-Kaas and B. Owren, Computations in a free Lie algebra, Phil. Trans. Royal Soc. A, 357 (1999), pp. 957–981.
  * [14] J. Sanz-Serna, Geometric integration, in The State of the Art in Numerical Analysis (York, 1996), vol. 63 of Inst. Math. Appl. Conf. Ser. New Ser., New York, 1997, Oxford Univ. Press, pp. 121–143.
  * [15] M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations, Phys. Lett. A, 146 (1990), pp. 319–323.
  * [16] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A, 150 (1990), pp. 262–268.