Higher Order Decompositions of Ordered Operator Exponentials

I Introduction

Decompositions of operator exponentials are widely used to approximate operator exponentials that arise in both physics and applied mathematics. These approximations are used because it is often difficult to exponentiate an operator directly, even if the operator is a sum of a sequence of operators that can be easily exponentiated individually. The goal in a decomposition method is to approximate an exponential of a sum of operators as a product of operator exponentials. Particular examples of decompositions include the Trotter formula, and the related Baker-Campbell-Hausdorff formula as well as Suzuki decompositions suzuki:physa ; suzuki:mathphys . These approximations have found use in many fields including quantum Monte-Carlo calculations suzuki:physa , quantum computing nielsenchuang and classical dynamics chin:integrator to name a few.

The problem is to solve for the operator $U$ given by

where $H$ is a linear operator and $\lambda$ and $\mu$ are real numbers. In general, computing $U$ can be difficult regardless of whether $H$ is dependent on $\lambda$. The case of $\lambda$-dependence makes the problem significantly more complicated and is the focus of this work, whereas the case of $\lambda$-independent $H$ has been well studied. The Lie-Trotter formula gives a simple solution, whereas more efficient higher-order solutions are given by the Lie-Trotter-Suzuki (LTS) product formulae suzuki:mathphys ; sanderssim .

A direct approach to solving $U$ in the $\lambda$-independent case is first to diagonalize $H$, then solve the differential equation (1) by direct exponentiation. The scenario we consider is that this is not possible, and instead we are given a set of $m$ operators $\{H_{j}:j=1,\ldots,m\}$, where it is possible to exponentiate of each of these operators, and

The evolution operator can then be approximated by a product of exponentials of $H_{j}$ for some sequence of $\{j_{i}\}$ and intervals $\{\Delta\lambda_{i}\}$,

The goal is to make an intractable calculation of $U$ tractable by approximating $U$ as a finite-length product of efficiently calculated exponentials. The complexity of the calculation can then be quantified by the number of exponentials $N$, and the scaling of $N$ in terms of the parameter difference $\Delta\lambda=\lambda-\mu$.

The case of $\lambda$-dependence makes the problem harder because, rather than the usual operator exponential of $H$, the solution is an ordered exponential. Diagonalization techniques are not directly applicable to solving ordered exponentials, and methods such as using the ordered product of exponentials are needed. That is, $U$ is approximated by an expression of the form

We consider the case where $H$ may be $\lambda$-dependent and where $H$ is a sum as in (2). One approach would be to replace Eq. (3) with a product formula of ordered exponentials. There would still remain the problem of evaluating the ordered exponentials, which would require an approach such as (4). A simpler approach is to use (3), but choose appropriate values of $\lambda$ at which to evaluate the $H_{j_{i}}$. The approximation is then

Suzuki provides an efficient method for approximating ordered exponentials in this way suzuki:mathphys . The method given by Suzuki is in terms of a time-displacement operator, but is equivalent. For many applications it is desirable to be able to place upper bounds on the error that can be obtained. Suzuki derives an order scaling, but not an upper bound on the error. Berry et al. find an upper bound on the error, but only in the case without $\lambda$-dependence, and for $H$ antihermitian (corresponding to Hamiltonian evolution) sanderssim .

Here we prove upper bounds on the error in the general case where $H$ can depend on $\lambda$, and is not restricted to be antihermitian. Whereas Suzuki uses a time-displacement operator, we provide a proof entirely without the use of this operator. The time-displacement operator is problematic, because it is unclear how it acts for $\lambda$-dependence that is non-analytic. We find that, provided all derivatives of the $H_{j}(\lambda)$ exist, Suzuki’s result holds. If there are derivatives that do not exist, then Suzuki’s result does not necessarily hold; we demonstrate this via a counterexample. We solve the case where derivatives may not exist, and find that Suzuki’s approach can still give better scaling of $N$ with $\Delta\lambda$, although the scaling that can be obtained is limited by how many times the $H_{j}(\lambda)$ are differentiable.

In Sec. II we give the background for Trotter product formulae in detail. In Sec. III we review Suzuki’s decomposition methods, and provide our form of Suzuki’s recursive method. In Sec. IV we introduce our terminology and present our main result. Then we rigorously prove the scaling of the error in Sec. V, and place an upper bound on the error in Sec. VI. We then use the error bounds in Sec. VII to find the appropriate order of the integrator to use.

II Trotter formulae

Typically there are two different scenarios that may be considered. First, one may consider a short interval $\Delta\lambda$; the goal is then to obtain error that decreases rapidly as $\Delta\lambda\to 0$. Alternatively the interval $\Delta\lambda$ may be long, and the goal is to obtain an approximation to within a certain error with as few exponentials as possible. For example, given $\lambda$-independent operators $A$ and $B$, it holds that

This gives an accurate approximation for small $\Delta\lambda$. For large $\Delta\lambda$, we may use Eq. (6) to derive the Trotter formula

This order of error is obtained because the error for interval $\Delta\lambda/n$ is $O((\Delta\lambda/n)^{2})$. Taking the power of $n$ then gives $n$ times this error if the norm of $\exp[(A+B)\Delta\lambda]$ is at most one for any $\Delta\lambda>0$, resulting in the error shown in Eq. (7). To obtain a given error $\epsilon$, the value of $n$ must then scale as $O(\Delta\lambda^{2}/\epsilon)$. The goal is to make the value of $n$ needed to achieve a given accuracy as small as possible ($n$ is proportional to the total number of exponentials).

More generally, for a sum of an arbitrary number of operators $H_{j}$, similar formulae give the same scaling. To obtain better scaling, one can use a different product of exponentials. The Lie-Trotter-Suzuki product formulae suzuki:mathphys replace the product for short $\Delta\lambda$ with another that gives error scaling as $O(\Delta\lambda^{p+1})$.

It can be seen that splitting large $\Delta\lambda$ into $n$ intervals as in Eq. (7) yields an error scaling as $O(\Delta\lambda^{p+1}/n^{p})$ if the norm of $U$ is at most one for any $\lambda$. It may at first appear that this gives worse results for large $\Delta\lambda$ due to the higher power. In fact, there is an advantage due to the fact that a higher power of $n$ is obtained. The value of $n$ required to achieve a given error then scales as $O(\Delta\lambda^{1+1/p}/\epsilon^{1/p})$. Therefore, for large $\Delta\lambda$, increasing $p$ gives scaling of $N$ that is close to linear in $\Delta\lambda$.

Similar considerations hold for the case of ordered exponentials (i.e. with $\lambda$-dependence). Huyghebaert and De Raedt showed how to generalize the Trotter formula to apply to ordered operator exponentials huyghebaert:timedeptrotter . Their formula has a decomposition error that is $O(\Delta\lambda^{2})$, but requires that the integrals of $A(u)$ and $B(u)$ are known. Subsequently Suzuki developed a method to achieve error that scales as $O(\Delta\lambda^{p+1})$ for some ordered exponentials suzuki:jacad , and does not require the integrals of $A$ and $B$ to be known. We find that, in contrast to the $\lambda$-independent case, it is not necessarily possible to obtain scaling as $O(\Delta\lambda^{p+1})$ for arbitrarily large $p$. It is possible if derivatives of all orders exist. If there are higher-order derivatives that do not exist, then it is still possible to use Suzuki’s method to obtain error scaling as $O(\Delta\lambda^{p+1})$ for some values of $p$, but the maximum value of $p$ for which this scaling can be proven depends on what orders of derivatives exist.

III Suzuki Decompositions

In this section we explain Suzuki decompositions in more detail. In general, decompositions are of the form, as in (5),

Here we write the final parameter $\lambda$ as $\mu+\Delta\lambda$, to emphasize the dependence on $\Delta\lambda$ ($=\lambda-\mu$). There are many different types of decompositions, but the type that we focus on in this paper is symmetric decompositions because all Suzuki decompositions are symmetric.

An important method for generating symmetric decompositions is due to Suzuki suzuki:timeorderbook ; suzuki:jacad , which we call Suzuki’s recursive method due to its similarity to the method presented by Suzuki in suzuki:mathphys . Furthermore we call any decomposition formula that is found using this method a Suzuki decomposition.

Suzuki’s recursive method takes a symmetric decomposition formula $U_{p}(\mu+\Delta\lambda,\mu)$, that approximates an ordered operator exponential $U(\mu+\Delta\lambda,\mu)$ with an approximation error that is at most proportional to $\Delta\lambda^{2p+1}$ as input, and outputs a symmetric approximation formula $U_{p+1}(\mu+\Delta\lambda,\mu)$ with an error that is often proportional to $\Delta\lambda^{2p+3}$. The approximation $U_{p+1}(\mu+\Delta\lambda,\mu)$ is found using the following recursion relations,

with $s_{p}\equiv\left(4-4^{1/(2p+1)}\right)^{-1}$.

Suzuki’s recursive method does not actually approximate $U(\mu+\Delta\lambda,\mu)$ but rather it builds a higher order approximation formula out of a lower order one. Therefore this method can only be used to approximate $U(\mu+\Delta\lambda,\mu)$ if it is seeded with an appropriate initial approximation. A convenient approximation formula based on Suzuki’s recursive method is the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula which is defined as follows.

The $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula for the operator $H(u)=\sum_{j=1}^{m}H_{j}(u)$ and the interval $[\mu,\mu+\Delta\lambda]$ is defined to be $U_{k}(\mu+\Delta\lambda,\mu)$, which is found by using

as an initial approximation and by applying Suzuki’s recursive method to it $k-1$ times.

Based on Suzuki’s analysis suzuki:timeorderbook ; suzuki:jacad $U_{k}$ should have approximation error that is proportional to $\Delta\lambda^{2k+1}$. Hence if $\Delta\lambda$ is sufficiently small, then the formula should be highly accurate. One might think that it would be advantageous to increase $k$ without limit, in order to obtain increasingly accurate approximation formulae. This is not the case, because the number of terms in the formula increases exponentially with $k$. The best value of $k$ to use can be expected to depend on the desired accuracy, as well as a range of other parameters sanderssim .

IV Sufficiency Criterion for Decomposition

Suzuki’s recursive method is a powerful technique for generating high-order decomposition formulae for ordered operator exponentials. The $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula in particular seems to be well suited for approximating ordered operator exponentials that appear in quantum mechanics and in other fields; furthermore it appears that these formulae should be applicable to approximating the ordered exponentials of any finite dimensional operator $H$. However it turns out that Suzuki’s recursive method does not always generate a higher order decomposition formula from a lower order one.

We show this using the example of the operator $H_{2}(u)=u^{3}\sin(1/u)\openone$. For this operator the second order Lie-Trotter-Suzuki product formula is not an approximation whose error as measured by the 2-norm is $O(\Delta\lambda^{5})$. In Figure 1 we see that the error is proportional to $\Delta\lambda^{4}$ for the operator $H_{a}(u)$, rather than the $\Delta\lambda^{5}$ scaling that we expect and observe for the analytic operator $H_{b}(u)=\cos(u)$. This shows that the second order Lie-Trotter-Suzuki product formula is not as accurate as may be expected for some non-analytic operators.

Our analysis will show that this discrepancy arises from the fact that $H_{2}(u)=u^{3}\sin(1/u)$ is not smooth enough for the second order Lie-Trotter-Suzuki formula to have an error which is $O(\Delta\lambda^{5})$. In the subsequent discussion we will need to classify the smoothness of the operators that arise in decompositions. We use the smoothness criteria $2k$-smooth and $\Lambda$-$2k$-smooth, which we define below.

The set of operators $\{H_{j}:j=1,\dots,m\}$ is $P$-smooth on the interval $[\mu,\lambda]$ if for each $H_{j}$ the quantity $\|\partial_{u}^{P}H_{j}(u)\|$ is finite on the interval $[\mu,\lambda]$.

Here, and throughout this paper, we define $\|\cdot\|$ to be the 2-norm. Also if $\{H_{j}\}$ is $P$-smooth for every positive integer $P$, we call $\{H_{j}\}$ $\infty$-smooth.

This condition is not precise enough for all of our purposes. For our error bounds we need to introduce the more precise condition of $\Lambda$-$P$-smoothness. This condition is useful because it guarantees that if the set $\{H_{j}\}$ is $\Lambda$-$P$-smooth and $p\leq P$ then $\|H^{(p)}(u)\|\leq\Lambda^{p}$. This property allows us to write our error bounds in a form that does not contain any of the derivatives of $H$ individually, but rather in terms of $\Lambda$ which upper bounds the magnitude of any of these derivatives. We formally define this condition below.

The set of operators $\{H_{j}:j=1,\dots,m\}$ is $\Lambda$-$P$-smooth on the interval $[\mu,\lambda]$ if $\{H_{j}\}$ is $P$-smooth and $\Lambda\geq\sup_{p=0,1,...,P}\left(\sup_{u\in[\mu,\lambda]}\left(\sum_{j=1}^{m}\|H_{j}^{(p)}(u)\|\right)^{1/(p+1)}\right)$.

For example if $\{H(u)\}=\{\sin(2u)\openone\}$, where $\openone$ is the identity operator, then using Definition 4 $\{H(u)\}$ is $2^{2/3}$-$2$-smooth on the interval $[0,\pi]$ because the largest value $\|H(u)^{(p)}\|^{1/(p+1)}$ takes is $2^{2/3}$, for $p=0,1,2$. It is also $2$-$2$-smooth because $2^{2/3}\leq 2$, furthermore since $\|H(u)^{(p)}\|^{1/(p+1)}<{2}$ for all positive integers $p$ then $\{H(u)\}$ is also ${2}$-$\infty$-smooth.

Using this measure of smoothness we can then state the following theorem, which is also the main theorem in this paper.

We prove Theorem 1 in several steps, the details of which are spread over Secs. V and VI. In Sec. V we construct the Taylor series for an ordered operator exponential, and use this series to prove that the Lie-Trotter-Suzuki product formula can generate an approximation whose error is $O(\Delta\lambda^{2k+1})$, if $\{H_{j}\}$ is $2k$-smooth on the interval $[\mu,\mu+\Delta\lambda]$. In Sec. VI we use the order estimates in Sec. V to obtain upper bounds on the error. The result of Theorem 1 then follows by counting the number of exponentials needed to make the error bound less than $\epsilon$. Finally in Sec. VII we show that if $k$ is chosen appropriately, then $N$ scales almost optimally with $\Delta\lambda$ if there exists a value of $\Lambda$ such that $\{H_{j}\}$ is $\Lambda$-$\infty$-smooth on $[\mu,\mu+\Delta\lambda]$ for every $\Delta\lambda>0$.

V Decomposing Ordered Exponentials

In this section we present a new derivation of Suzuki’s recursive method. Our derivation has the advantage that it can be rigorously proven that if $\{H_{j}\}$ is $2k$-smooth, where $H(u)=\sum_{j=1}^{m}H_{j}(u)$, then the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula will have an error of $O(\Delta\lambda^{2k+1})$. We show this in three steps. We first give an expression for the Taylor series expansion of a ordered exponential $U(\mu+\Delta\lambda,\mu)$. Then using this expression for the Taylor series, we show in Theorem 2 that Suzuki’s recursive method can be used to generate approximations to $U(\mu+\Delta\lambda,\mu)$ that invoke an error that is $O(\Delta\lambda^{2k+1})$ if $\{H_{j}\}$ is $2k$-smooth. Finally we show in Corollary 1 that if $\{H_{j}\}$ is $2k$-smooth then the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula has an error that is at most proportional to $\Delta\lambda^{2k+1}$.

It is convenient to expand $U$ in a Taylor series of the form

If $H$ is not analytic, then this Taylor series must be truncated, and the error can be bounded by the following lemma.

where $T_{p}(u)$ is defined by the recursion relation, $T_{p+1}(u)\equiv T_{p}(u)H(u)+\partial_{t}T_{p}(u)$, with $T_{0}\equiv\openone$ chosen to be the initial condition.

We then use (14) and Taylor’s Theorem to conclude that

If $H=\sum_{j=1}^{m}H_{j}$ where the set $\{H_{j}\}$ is for a fixed $p$, $2(p+1)$-smooth on the interval $[\mu,\mu+\Delta\lambda]$ and $U_{p}(\mu+\Delta\lambda,\mu)$ is a symmetric approximation formula such that $\|U_{p}(\mu+\Delta\lambda,\mu)-U(\mu+\Delta\lambda,\mu)\|\in O(\Delta\lambda^{2p+1})$, and $U_{p+1}(\mu+\Delta\lambda,\mu)$ is found by applying Suzuki’s recursive method on $U_{p}(\mu+\Delta\lambda,\mu)$, then

In this proof we compare the Taylor series of $U$ to that of $U_{p}$ and show that by choosing $s_{p}$ appropriately will cause both the terms proportional to $\Delta\lambda^{2p+2}$ and $\Delta\lambda^{2p+3}$ to vanish.

By expanding the recursive formula in Lemma 1 we see that a Taylor polynomial can be constructed for $U$ whose difference from $U$ is $O(\Delta\lambda^{2p+3})$ because $\{H_{j}\}$ is $2(p+1)$-smooth on $[\mu,\mu+\Delta\lambda]$. A similar polynomial can be constructed for $U_{p}$ by Taylor expanding each $H_{j}$ that appears in the exponentials in $U_{p}$, and then expanding each of these exponentials. Then because $\{H_{j}\}$ is $2(p+1)$-smooth, Taylor’s Theorem implies that this polynomial can be constructed such that the difference between it and $U_{p}$ is $O(\Delta\lambda^{2p+3})$. Therefore since $\|U-U_{p}\|\in O(\Delta\lambda^{2p+1})$ there exist operators $C$ and $E$ that are independent of $\Delta\lambda$, such that

We then use the above equation to write $U_{p+1}$ as

Then we see that if $s_{p}=(4-4^{1/(2p+1)})^{-1}$, then the terms of order ${2p+1}$ in the above equation vanish. Hence the error invoked using $U_{p+1}$ instead of $U$ is $O(\Delta\lambda^{2p+2})$ with this choice of $s_{p}$.

This implies that the norm of the difference between $U$ and $U_{p+1}$ is proportional to $\Delta\lambda^{2p+3}$, which concludes our proof of Theorem 2. $\sqcap$$\sqcup$

Now that we have proved that Suzuki’s recursive method will generate a higher order decomposition formula from a lower order one if $H$ is the sum of the elements from a sufficiently smooth set $\{H_{j}\}$, we now show that using the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula invokes an error that is proportional to $\Delta\lambda^{2k+1}$ if $\{H_{j}\}$ is $2k$-smooth.

Let $H(u)=\sum_{j=1}^{m}H_{j}(u)$ where the set $\{H_{j}\}$ is $2k$-smooth on the interval $[\mu,\mu+\Delta\lambda]$ and let $U(\mu+\Delta\lambda,\mu)$ be the ordered operator exponential generated by $H$. Then if $U_{k}(\mu,\mu+\Delta\lambda)$ is the $k^{\text{th}}$ order Lie-Trotter-Suzuki formula, then

Our proof of the corollary follows from an inductive argument on $k$. The validity of the base case can be verified by using Lemma 1. More specifically, since $\{H_{j}\}$ is $2k$-smooth on $[\mu,\mu+\Delta\lambda]$ and since $k\geq 1$, then $H$ is at least three times differentiable on that interval. This means that we can use Lemma 1 to say that

This expansion is also obtained by Taylor expanding $\exp(H(\mu+\Delta\lambda/2)\Delta\lambda)$ to third order, so

Since $U_{1}(\mu+\Delta\lambda,\mu)$ is the Lie-Trotter formula for a constant $H$ equal to $H(\mu+\Delta\lambda/2)$ it follows that,

It follows from the above equations and from the triangle inequality that the norm of the difference between $U_{1}$ and $U$ is at most proportional to $\Delta\lambda^{3}$.

Since we have shown that $U_{1}(\mu+\Delta\lambda,\mu)$ is a symmetric approximation formula whose error is $O(\Delta\lambda^{3})$, it then follows from Theorem 2 and induction, that if $\{H_{j}\}$ is $2k$-smooth then a symmetric approximation formula whose error is $O(\Delta\lambda^{2k+1})$ can be constructed from $U_{1}(\mu+\Delta\lambda,\mu)$ by applying Suzuki’s recursive method to it $k-1$ times. $\sqcap$$\sqcup$

We have shown in this section that if $\{H_{j}\}$ is $2k$-smooth on the interval $[\mu,\mu+\Delta\lambda]$ and if $p\leq k$ then Suzuki’s recursive method can be used to create a symmetric decomposition whose error is $O(\Delta\lambda)^{2p+1}$ out of a symmetric decomposition whose error is $O(\Delta\lambda)^{2p-1}$. Then we have used this fact to show that the norm of the difference between $U(\mu+\Delta\lambda,\mu)$ and the $k^{\text{th}}$ order Lie-Trotter-Suzuki formula is $O(\Delta\lambda^{2k+1})$. In the following section we strengthen this result by providing an upper bound on the error invoked by using the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula.

VI Error Bounds and Convergence for Decomposition

We showed in Sec. V that if the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula is used in the place of the ordered operator exponential of $H$, then an error is incurred that is at most proportional to $\Delta\lambda^{2k+1}$ if $H=\sum_{j=1}^{m}H_{j}$ and the set of operators $\{H_{j}\}$ is sufficiently smooth. We also showed that a sufficient condition for smoothness of the set $\{H_{j}\}$ is a condition that we called $2k$-smooth, where this condition is defined is Definition 3. In this section we extend that result by finding upper bounds on the error invoked in using the Lie-Trotter-Suzuki product formula to approximate ordered operator exponentials if $\{H_{j}\}$ is $\Lambda$-$2k$-smooth. Unlike the previous section, here we assume that $\max_{x>y}\|U(x,y)\|$ is at most one. This assumption is important because it ensures that our error bounds are not exponentially large. Our work can be made applicable to the case where this norm is greater than one by re-normalizing $U$. We discuss the implications of this in Appendix B.

We first provide in this section an upper bound on the error invoked in using the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula to approximate the ordered operator exponential $U(\mu+\Delta\lambda,\mu)$ if $\Delta\lambda$ is sufficiently short. We then use this result to upper bound the error if $\Delta\lambda$ is not short. More specifically, we show that for every $1\leq\epsilon>0$ and $\Delta\lambda>0$ there exists an integer $r$ such that

if $\{H_{j}(u)\}$ is $2k$-smooth on the interval $[\mu,\mu+\Delta\lambda]$. Finally by multiplying the number of exponentials in each $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula by $r$, we find the number of exponentials used in the product in (28). We then use this result to prove Theorem 1.

Our upper bound on the error invoked by using a single $U_{k}$ to approximate the ordered operator exponential $U(\mu+\Delta\lambda,\mu)$ is given in Theorem 3. Before stating Theorem 3 we first define the following terms. Since the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula is a product of $2m5^{k-1}$ exponentials, we can express this product as $\prod_{c=1}^{2m5^{k-1}}\exp(H_{j_{c}}(\mu_{c})\Delta\lambda_{c})$. We then use this expansion to define the following two useful quantities.

We define $q_{c,2k}\equiv\frac{\Delta\lambda_{c}}{\Delta\lambda}$ and also define $Q_{k}\equiv\max_{c}|q_{c,2k}|$.

It can be shown that $Q_{1}=1/2$ and that if $p>1$ then $Q_{p}=|1-4s_{1}|\cdots|1-4s_{p-1}|$. We show in Appendix A that for any integer $p$ that $Q_{p}\leq 2p/3^{p}$, implying that $Q_{p}$ decreases exponentially with $p$. We use this definition of $Q_{k}$ in the following Theorem, that gives an upper bound on the difference between the ordered operator exponential $U(\mu+\Delta\lambda,\mu)$, and the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula $U_{k}(\mu+\Delta\lambda,\mu)$.

Let $H(u)=\sum_{j=1}^{m}H_{j}(u)$, let $\{H_{j}(u)\}$ be $\Lambda$-$2k$-Suzuki-smooth on the interval $[\mu,\mu+\Delta\lambda]$, and let $\max_{x>y}\|U(x,y)\|\leq 1$. Then if $2\sqrt{2}(5)^{k-1}Q_{k}\Lambda\Delta\lambda\leq 1/2$, it follows that

The proof of Theorem 3 requires us to first prove two Lemmas before we can conclude that the theorem is valid. We now introduce some notation to state these lemmas concisely. Since we have assumed that $\{H_{j}\}$ is $2k$-smooth, Theorem 2 implies that the difference between $U(\mu+\Delta\lambda,\mu)$ and $U_{k}(\mu+\Delta\lambda,\mu)$ is $O(\Delta\lambda)^{2k+1}$. Then using this fact, we know that we only need to compare the terms of $O(\Delta\lambda^{2k+1})$ to bound the difference between $U$ and $U_{k}$. We introduce the following notation to denote only those terms that do not necessarily cancel.

If the operator $A(\Delta\lambda)$ can be written as $A(\Delta\lambda)=\sum_{p=0}^{2k}A_{p}\Delta\lambda^{p}+R(\Delta\lambda)$ where the norm of $R(\Delta\lambda)$ is $O(\Delta\lambda)^{2k+1}$, then we define $\mathbf{R}_{2k}[A(\Delta\lambda)]$ to be the norm of $R(\Delta\lambda)$.

This definition simply means that $\mathbf{R}_{2k}$ is the error term for a Taylor expansion to order $2k$. Then using this definition, it follows from the triangle inequality that the norm of the difference between $U$ and $U_{k}$ is at most

Our proof of Theorem 3 then follows from (29) and upper bounds that we place on $\mathbf{R}_{2k}[U(\mu+\Delta\lambda,\mu)]$ and $\mathbf{R}_{2k}[U_{k}(\mu+\Delta\lambda,\mu)]$. Our bound on $\mathbf{R}_{2k}[U(\mu+\Delta\lambda,\mu)]$ follows directly from Lemma 1, but the bound on $\mathbf{R}_{2k}[U_{k}(\mu+\Delta\lambda,\mu)]$ does not. We will provide the latter upper bound in Lemma 3, but first we provide Definition 7 and Lemma 2.

Let $k$ be an integer and let $H=\sum_{j=1}^{m}H_{j}$ then $U_{k}(\mu+\Delta\lambda,\Delta\lambda)$ can be written as a product of the form $\prod_{c=1}^{2m5^{k-1}}\exp(H_{j_{c}}(\mu_{c})\Delta\lambda_{c})$. We then define $X_{p}$ for $p<2k$ to be

Here the quantity $q_{c,2k}$ is given in Definition 5.

Then using this definition our lemma can be expressed as follows.

Let $H(u)=\sum_{j=1}^{m}H_{j}(u)$ and let the set $\{H_{j}\}$ be $2k$-smooth on the interval $[\mu,\mu+\Delta\lambda]$, then the norm of the difference between $U_{k}(\mu+\Delta\lambda,\mu)$ and its Taylor series in powers of $\Delta\lambda$ truncated at order ${2k}$, is bounded above by

We begin our proof of Lemma 2 by writing $U_{k}$ as a product of $2m5^{k-1}$ exponentials and use Taylor’s theorem to write $U_{k}$ as

We introduce the terms $v_{c}=(\mu_{c}-\mu)/\Delta\lambda$ and $q_{c,k}=\frac{\Delta\lambda_{c}}{\Delta\lambda}$ and use them to write $U_{k}(\mu+\Delta\lambda,\mu)$ as

We now prove the lemma by placing an upper bound on $\mathbf{R}_{2k}[U_{k}(\mu+\Delta\lambda,\mu)]$ by expanding this equation in powers of $\Delta\lambda$, while retaining only those terms of order ${2k+1}$ and higher. As mentioned previously, the lower order terms are irrelevant since Theorem 2 guarantees that they cancel.

By expanding the exponentials in (34), taking the norm, using the triangle inequality, upper bounding each of the norms present in the expansion, and collecting terms again, we find that an upper bound on $\mathbf{R}_{2k}[U_{k}(\mu+\Delta\lambda,\mu)]$ is

This equation can be simplified by substituting the constants $X_{p}$ into it. These constants are introduced in Definition 7. After this substitution our upper bound becomes

We use Lemma 2 to provide an upper bound on the sum of the norm of all terms in the Taylor expansion of $U_{k}(\mu+\Delta\lambda,\mu)$ which are of order ${2k+1}$ or higher. This bound is given in the following lemma.

Let $H(u)=\sum_{j=1}^{m}H_{j}(u)$ where the set $\{H_{j}\}$ is $\Lambda$-$2k$-smooth on the interval $[\mu,\mu+\Delta\lambda]$ and let $2\sqrt{2}(5)^{k-1}Q_{k}\Lambda\Delta\lambda\leq\frac{1}{2}$. Then the norm of the difference between $U_{k}(\mu+\Delta\lambda,\mu)$ and its Taylor series in $\Delta\lambda$ truncated at order ${2k}$, is upper bounded by

To simplify the following discussion we introduce $\Gamma_{2k}$, defined by

We first find an upper bound on $\Gamma_{2k}$. Now, by Definition 7,

In Appendix A, we show that $q_{c,k}\leq Q_{k}\leq\frac{2k}{3^{k}}$ for all $c$, where $Q_{k}$ is an upper bound on the $q_{c,k}$ given in Definition 5. For a $2k^{\text{th}}$ order Lie-Trotter-Suzuki product formula, $\frac{\mu_{c}-\mu}{\Delta\lambda}\leq 1$. Plugging these two bounds into the above inequality and using the fact that each element of $\{H_{j}\}$ occurs $2(5^{k-1})$ times in $U_{k}$ yields,

Since we assume that $\{H_{j}\}$ is $\Lambda$-$2k$-smooth, then $\sum_{j=1}^{m}\|H_{j}^{(p)}(\tau)\|\leq\Lambda^{p+1}$, so X_{p}^{1/(p+1)}\leq\big{(}2(5)^{k-1}Q_{k}\big{)}^{1/(p+1)}\Lambda. We also show in Appendix A that $Q_{k}\geq\frac{1}{2}\frac{1}{3^{k-1}}$, implying that $2(5)^{k-1}Q_{k}\geq 1$, and thus that $X_{p}^{1/(p+1)}\leq 2(5)^{k-1}Q_{k}\Lambda$. Since this upper bound holds for all $0\leq p\leq 2k$, we conclude that

We begin the main inequality in Lemma 3, by expanding $\mathbf{R}_{2k}\left[\exp\left(\sum_{p=0}^{2k}\frac{X_{p}}{p!}\Delta\lambda^{p+1}\right)\right]$ in powers of $\Delta\lambda$ and using $\frac{1}{p!}\leq\frac{(\sqrt{2})^{p+1}}{p+1}$ and $X_{p}\leq\Gamma_{2k}^{p+1}$. Writing the resulting expansion as an exponential we find that

It then follows that the right hand side of the above expression is upper bounded by

Using the Taylor expansion of $\ln(1-x)$, we rewrite this as

Provided $\sqrt{2}\Gamma_{2k}\Delta\lambda\leq\frac{1}{2}$, this is upper bounded by $2(\sqrt{2}\Gamma_{2k}\Delta\lambda)^{2k+1}$. Plugging inequality (41) into Eq. (44) then gives that

The lemma follows by applying Lemma 2. $\sqcap$$\sqcup$

Now that we have proven Lemma 3 we have an upper bound on $\mathbf{R}_{2k}[U_{k}(\mu+\Delta\lambda,\mu)]$. We now use this upper bound to prove Theorem 3.

Our proof of Theorem 3 begins by recalling the fact that

We then place an upper bound on $\mathbf{R}_{2k}\left[U(\mu+\Delta\lambda,\mu)\right]$ using Lemma 1. Using the notation of Lemma 1 we write the Taylor series of $U(\mu+\Delta\lambda,\mu)$ as $\sum_{p}T_{p}\Delta\lambda^{p}/p!$. We then use the assumption that $\|U(\mu+\Delta\lambda,\mu)\|$ is less than one, to show from Lemma 1 that $\mathbf{R}_{2k}\left[U(\mu+\Delta\lambda,\mu)\right]$ is at most

Using the recursive relations in Lemma 1 it follows that $T_{2k+1}$ can be written as a sum of $(2k+1)!$ terms that are each products of $H$ and its derivatives. Then since $\{H_{j}:j=1,\dots,m\}$ is $\Lambda$-$2k$-smooth and $H=\sum_{j=1}^{m}H_{j}$, it follows from Definition 4 that $\|H^{(p)}(u)\|\leq\Lambda^{p}$ for all $u$ in the interval $[\mu,\mu+\Delta\lambda]$. It can then be verified that each term in $T_{2k+1}$ must have a norm that is less than $\Lambda^{2k+1}$. Therefore it follows that $\|T_{2k+1}\|\leq(2k+1)!\Lambda^{2k+1}$ and hence

Using Eq. (46) and Lemma 3 we see that if $2\sqrt{2}(5)^{k-1}Q_{k}\Lambda\Delta\lambda\leq 1/2$ then an upper bound on the sum of $\mathbf{R}_{2k}\left[U(\mu+\Delta\lambda,\mu)\right]$ and $\mathbf{R}_{2k}\left[U_{k}(\mu+\Delta\lambda,\mu)\right]$ is

We then replace this upper bound with the following simpler upper bound

This is the claim in Theorem 3, and hence we have proven the theorem. $\sqcap$$\sqcup$

The error bound in Theorem 3 is vital to our remaining work, because it provides us with an upper bound on the error invoked by approximating an ordered operator exponential by $U_{k}(\mu+\Delta\lambda,\mu)$ if $\Delta\lambda$ is short. We will now show a method to devise accurate approximations to the ordered operator exponential $U(\mu+\Delta\lambda,\mu)$ even if $\Delta\lambda$ is not short. Our approach is similar to that used by Berry et al. in sanderssim and that used by Suzuki in suzuki:commmathphys ; we split the ordered exponential into a product of ordinary exponentials, each of which has a short duration. However to do so we need to present a method to relate the error invoked by using one $U_{k}$ to the error invoked by using a product of them. This result is provided in the following lemma.

If $\|A_{p}-B_{p}\|\leq\delta/P$ where $\delta$ is a positive number less than $1/2$ and $\|A_{p}\|\leq 1$ for every $p\in\{1,2,\cdots,P\}$ then the product $\|\prod_{p=1}^{P}A_{p}-\prod_{p=1}^{P}B_{p}\|\leq 2\delta$.

Our proof begins by assuming that there exists some integer $q$ such that

We then prove Lemma 4 by using induction on $q$. The proof of the base case follows from $\|A_{p}-B_{p}\|\leq\delta/P$. We then begin to prove the induction step by noting that from $\|A_{p}-B_{p}\|\leq\delta/P$ there exists an operator $C$ with norm at most one, such that $B_{q+1}=A_{q+1}+(\delta/P)C$. Then by making this substitution and using the triangle inequality it follows that

Then because $\|A_{p}\|\leq 1$ and $\|B_{p}\|\leq 1+\delta/P$ it can be verified using our induction hypothesis that the left hand side of Equation (52) is bounded above by

This proves our induction step, and so it follows that $\|\prod_{p=1}^{P}A_{p}-\prod_{p=1}^{P}B_{p}\|\leq\delta(1+\delta/P)^{P-1}$ by using induction on $q$ until $q=P$. The proof of the Lemma then follows from the fact that if $\delta\leq 1/2$ then $(1+\delta/P)^{P-1}\leq 2$. $\sqcap$$\sqcup$

Using Lemma 4 we can now place an upper bound on the error for decompositions with longer $\Delta\lambda$.

If $H(u)=\sum_{j=1}^{m}H_{j}(u)$ is $2k$-Suzuki-smooth on the interval $[\mu,\mu+\Delta\lambda]$, and $\max_{x>y}\|U(x,y)\|\leq 1$ and $\epsilon\leq 3Q_{k}(5)^{k-1}\Lambda\Delta\lambda$, where $Q_{k}$ is given in Definition 6, and the positive integer $r$ is greater than

where $U_{k}$ is the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula, which we introduced in Definition 2.

Using the bound $\epsilon\leq 3Q_{k}(5)^{k-1}\Lambda\Delta\lambda$, we find using Eq. (54) that $3Q_{k}(5)^{k-1}\Lambda\Delta\lambda/r\leq 1/2$. Hence we can use Theorem 3 to obtain, for each $q=1,\ldots,r$,

Then from (54) we can see that, because $r\geq 2(3Q_{k}(5)^{k-1}\Lambda\Delta\lambda)^{1+1/2k}/{\epsilon^{1/2k}}$ it follows that

Then since both $\epsilon$ and $\max_{x>y}\|U(x,y)\|$ are less than one, the result of this lemma follows from Lemma 4. $\sqcap$$\sqcup$

Lemma 5 shows that if the maximum value of the norm of $U$ is one, then a product of $k^{\text{th}}$ order Lie-Trotter-Suzuki formulae converges to $U$ as $r$ increases. Furthermore we can also use this result to prove Theorem 1 by using the value of $r$ from this Lemma and multiplying it by the number of exponentials in each $U_{k}$ to find a bound on the number of exponentials that are needed to approximate $U(\mu+\Delta\lambda,\mu)$. This proof is presented below.

It can be verified from the definition of $U_{k}(\mu+\Delta\lambda,\mu)$ in Theorem 2 that there are at most $2m5^{k-1}$ exponentials in each $U_{k}$ and since at most $r$ different $U_{k}$ in are needed to approximate $U$ within an error of $\epsilon$ then if $\max_{x>y}\|U(x,y)\|$ and $\epsilon$ are at most one, it follows from Lemma 5 that the number of exponentials used to decompose $U(\mu+\Delta\lambda,\mu)$ is at most

if $\epsilon\leq 3Q_{k}5^{k-1}\Lambda\Delta\lambda$. We then use the upper bound from Appendix A, $Q_{k}\leq 2k/3^{k}$ and the fact that $(2k)^{1/2k}<1.5$ to show that,

Equation (61) is only valid if $\epsilon\leq 3Q_{k}5^{k-1}\Lambda\Delta\lambda$, this bound can be simplified by using the lower bound on $Q_{k}$ in (67). After substituting this lower bound we then find that $\epsilon\leq(9/10)(5/3)^{k}\Lambda\Delta\lambda$ also is sufficient to guarantee that (61) is valid, which proves our theorem. $\sqcap$$\sqcup$

Theorem 1 provides an upper bound on the number of exponentials that are needed to decompose an ordered operator exponential using a product of $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula, while guaranteeing that the approximation error is at most $\epsilon$. In the following section we present a formula that provides a reasonable value of $k$, for a particular set of values for $\epsilon,\Lambda$ and $\Delta\lambda$. Furthermore we show that if $\{H_{j}\}$ is $\Lambda$-$\infty$-smooth and if that formula for $k$ is used, then the number of exponentials used scales near optimally with $\Delta\lambda$.

VII Almost Linear Scaling

Reference sanderssim shows that there exist operator exponentials that, when decomposed into a sequence of $N$ exponentials, require that $N$ scale at least linearly with $\Delta\lambda$ for large $\Delta\lambda$. This implies that any decomposition method that does not use any special properties of the operator being exponentiated, will also require that $N$ scale at least linearly with $\Delta\lambda$. We now show that if there exists a $\Lambda$ such that the set of operators $\{H_{j}\}$ is $\Lambda$-$\infty$-smooth on $[\mu,\mu+\Delta\lambda]$ for every $\Delta\lambda>0$, then we can choose $k$ such that $N$ scales almost linearly in $\Delta\lambda$. Specifically, we show that $N/\Delta\lambda$ is sub-polynomial in $\Delta\lambda$, i.e., that $\lim_{\Delta\lambda\rightarrow\infty}\frac{N}{\Delta\lambda^{1+d}}=0$ for all constants $d>0$, provided that $\max_{x>y}\|U(x,y)\|\leq 1$.

It follows from Theorem 1 and property $\lceil a\rceil\leq a+1$ for any positive $a$, that if $\Lambda\Delta\lambda>1$ and $\epsilon\leq 1$ then

so that $\left(\frac{25}{3}\right)^{k_{0}}\geq\left(\frac{\Lambda\Delta\lambda}{\epsilon}\right)^{1/2k_{0}}$. Then

which is sub-polynomial in $\Delta\lambda$, though not poly-logarithmic in $\Delta\lambda$. In conclusion, if $H$ is $\Lambda$-$\infty$-smooth and if we choose $k=k_{0}$, then the number of exponentials needed to decompose a ordered operator exponential of $H$ using the $k^{\text{th}}$ order Lie-Trotter-Suzuki formula scales almost linearly in $\Delta\lambda$.

This choice of $k_{0}$ will cause $N$ to scale nearly linearly with $\Delta\lambda$ if $H$ is $\infty$-smooth; however if $\{H_{j}\}$ is only $2P$-smooth for some positive integer $P$, then we do not expect this because Theorem 1 cannot be used if $k_{0}>P$. Hence a reasonable choice of $k_{0}$ is

The choice of $k_{0}$ in (65) does not allow for near linear scaling of $N$ with $\Delta\lambda$, but it does cause $N$ to be proportional to $\Delta\lambda^{1+1/(2P)}$ in the limit of large $\Delta\lambda$, and causes $N$ to have the same scaling with $\Delta\lambda$ that a $\Lambda$-$\infty$-smooth $\{H_{j}\}$ would have if $\Delta\lambda$ is sufficiently short.

VIII Conclusions

We have presented in this paper a rigorous derivation of Suzuki’s recursive method for generating higher order approximations to ordered operator exponentials from a lower order formula. We have also shown that if $H(u)=\sum_{j=1}^{m}H_{j}(u)$ and $\{H_{j}\}$ is $2k$-smooth, a condition which we define in Definition 3, then Suzuki’s recursive method can be used to build approximation formulae for the ordered exponential of $H$ while invoking an error that is at most proportional to $\Delta\lambda^{2k+1}$. Furthermore we have shown that the $k^{\text{th}}$ order Lie-Trotter-Suzuki product formula has an error at most proportional to $\Delta\lambda^{2k+1}$ if $H$ is $2k$-smooth.

We have also shown that if $H$ is $\Lambda$-$2k$-smooth, which is a condition that we give in Definition 4, and if $\max_{x>y}\|U(x,y)\|\leq 1$ then the number of exponentials needed to approximate an ordered operator exponential of $H$ using a product of $k^{\text{th}}$ order Lie-Trotter-Suzuki formulae while invoking a total error of at most $\epsilon$ is bounded above by

Finally we have shown that if $\{H_{j}\}$ is $\Lambda$-$\infty$-smooth then for long simulations, the number of exponentials needed to approximate ordered operator exponentials of $H$ with an error of at most $\epsilon$ scales close to linearly with $\Delta\lambda$. Linear scaling has been shown to be optimal sanderssim , so in this sense our scheme is nearly optimal.

An important extension of this work will be to provide upper bounds on the number of exponentials used to decompose an ordered operator exponential when the user chooses step size adaptively, rather than the constant sized steps that we use in our derivation. Choosing these steps adaptively could lead to substantial improvements in the performance of our decomposition for certain ordered exponentials.

When proving bounds on the error introduced by our decomposition of an ordered exponential in Section VI, we use the quantities $Q_{k}$ defined by

for $k>1$ and $Q_{1}=\frac{1}{2}$, where $s_{k}=\frac{1}{4-4^{1/(2k+1)}}$ for all $k\geq 1$. We now show that $Q_{k}$ decreases exponentially in $k$,

The lower bound follows directly by noting that $s_{k}\geq\frac{1}{3}$ for all $k\geq 1$.

Set $a=2\ln(4)$, which is approximately $2.7726$. Using that $-x\geq\ln(1-x)$ for $0\leq x<1$, we then have that for $k\geq 1$,

Taking exponentials on either side yields,

Multiplying by 4 and subtracting 1 on either side gives,

Noting that $4s_{k}-1\geq s_{k}$ since $s_{k}\geq\frac{1}{3}$, and using that $s_{1}\leq\frac{2}{3}$, we conclude that

for $k\geq 2$, and by inspection that the inequality $Q_{k}\leq\frac{2k}{3^{k}}$ also holds for $k=1$.

Appendix B Norms larger than 1

In this work we have restricted the norm of $U(\lambda,\mu)$ to not exceed 1. This means that the eigenvalues of $H(\lambda)$ can have no positive real part. In the case where they do, then the analysis can be performed in the following way. Simply define the new operators

for some $\kappa(\lambda)$ such that the eigenvalues of $H^{\prime}(\lambda)$ have no positive real part. Then the result we have given in Theorem 1 will hold for $H^{\prime}$ and $\{H_{j}^{\prime}\}$ (provided we also define $\Lambda$ in terms of these operators). The difference between $H$ and $H^{\prime}$ simply corresponds to a normalization factor; i.e.

To approximate $U(\lambda,\mu)$ by a series of exponentials, we can simply use the series to approximate $U^{\prime}(\lambda,\mu)$, which gives

Thus the same series of exponentials can be used, except for a normalization factor. There is a difference in the final error that can be obtained, because

It might be imagined that the relative error can be kept below $\epsilon$ with similar scaling of $N$. That is, that $e^{\int_{\mu}^{\lambda}\kappa(x)dx}$ can be replaced with $\|U(\lambda,\mu)\|$ in (79). Unfortunately, that is not the case. The reason is that, due to the submultiplicativity of the operator norm, $\|U(\lambda,\mu)\|$ can be much smaller than $e^{\int_{\mu}^{\lambda}\kappa(x)dx}$.