Optimized Lie-Trotter-Suzuki decompositions for two and three non-commuting terms

TL;DR

This paper develops optimized Lie-Trotter-Suzuki decompositions up to sixth order for efficiently approximating operator exponentials in quantum and classical systems, improving accuracy over previous methods.

Contribution

It introduces new optimized decompositions up to order six for multiple non-commuting operators, with minimized error coefficients, enhancing simulation efficiency.

Findings

Optimized decompositions significantly reduce approximation errors.

Results are close to previous optimal decompositions for two terms.

Decompositions are applicable to lattice models with finite-range interactions.

Abstract

Lie-Trotter-Suzuki decompositions are an efficient way to approximate operator exponentials $\exp(t H)$ when $H$ is a sum of $n$ (non-commuting) terms which, individually, can be exponentiated easily. They are employed in time-evolution algorithms for tensor network states, digital quantum simulation protocols, path integral methods like quantum Monte Carlo, and splitting methods for symplectic integrators in classical Hamiltonian systems. We provide optimized decompositions up to order $t^6$. The leading error term is expanded in nested commutators (Hall bases) and we minimize the 1-norm of the coefficients. For $n=2$ terms, several of the optima we find are close to those in McLachlan, SlAM J. Sci. Comput. 16, 151 (1995). Generally, our results substantially improve over unoptimized decompositions by Forest, Ruth, Yoshida, and Suzuki. We explain why these decompositions are sufficient to efficiently simulate any one- or two-dimensional lattice model with finite-range interactions. This follows by solving a partitioning problem for the interaction graph.

Figures

Tables

Table 1. Table 1: Hall bases for the free Lie algebra generated by n = 2 𝑛 2 n=2 and 3 3 3 operators, respectively.

Degree	Hall basis elements
1	$A$,  $B$
2	$[A,B]$
3	$[A,[A,B]]$,  $[B,[A,B]]$
4	$[A,[A,[A,B]]]$,  $[B,[A,[A,B]]]$,  $[B,[B,[A,B]]]$
5	$[A,[A,[A,[A,B]]]]$,  $[B,[A,[A,[A,B]]]]$,  $[B,[B,[A,[A,B]]]]$,
	$[B,[B,[B,[A,B]]]]$,  $[[A,B],[A,[A,B]]]$,  $[[A,B],[B,[A,B]]]$
1	$A$,  $B$,  $C$
2	$[A,B]$,  $[A,C]$,  $[B,C]$
3	$[A,[A,B]]$,  $[A,[A,C]]$,  $[B,[A,B]]$,  $[B,[A,C]]$,
	$[B,[B,C]]$,  $[C,[A,B]]$,  $[C,[A,C]]$,  $[C,[B,C]]$
4	$[A,[A,[A,B]]]$,  $[A,[A,[A,C]]]$,  $[B,[A,[A,B]]]$,  $[B,[A,[A,C]]]$,  $[B,[B,[A,B]]]$,
	$[B,[B,[A,C]]]$,  $[B,[B,[B,C]]]$,  $[C,[A,[A,B]]]$,  $[C,[A,[A,C]]]$,  $[C,[B,[A,B]]]$,
	$[C,[B,[A,C]]]$,  $[C,[B,[B,C]]]$,  $[C,[C,[A,B]]]$,  $[C,[C,[A,C]]]$,  $[C,[C,[B,C]]]$,
	$[[A,B],[A,C]]$,  $[[A,B],[B,C]]$,  $[[A,C],[B,C]]$

Table 2. Table 2: Numbers d k subscript 𝑑 𝑘 d_{k} of Hall bases elements with degree k 𝑘 k for the free Lie algebra generated by n = 2 𝑛 2 n=2 and 3 3 3 operators, respectively.

Degree $k$	1	2	3	4	5	6	7	8	9	10	11	12
$n=2$	2	1	2	3	6	9	18	30	56	99	186	335
$n=3$	3	3	8	18	48	116	312	810	2184	5880	16104	44220

Table 3. Table 3: Hall basis elements of odd degree for the free Lie algebra generated by operators { A 1 , A 3 , A 5 , … } subscript 𝐴 1 subscript 𝐴 3 subscript 𝐴 5 … \{A_{1},A_{3},A_{5},\dotsc\} , where A k subscript 𝐴 𝑘 A_{k} has degree k 𝑘 k .

Degree	Hall basis elements
1	$A_1$
3	$A_3$
5	$A_5$,  $[A_1,[A_1,A_3]]$
7	$A_7$,  $[A_1,[A_1,A_5]]$,  $[A_3,[A_1,A_3]]$,  $[A_1,[A_1,[A_1,[A_1,A_3]]]]$
9	$A_9$,  $[A_1,[A_1,A_7]]$,  $[A_3,[A_1,A_5]]$,  $[A_5,[A_1,A_3]]$,  $[A_1,[A_1,[A_1,[A_1,A_5]]]]$,
	$[A_3,[A_1,[A_1,[A_1,A_3]]]]$,  $[[A_1,A_3],[A_1,[A_1,A_3]]]$,  $[A_1,[A_1,[A_1,[A_1,[A_1,[A_1,A_3]]]]]]$

Table 4. Table 4: Numbers d L , k subscript 𝑑 L 𝑘 d_{\text{L},k} of Hall bases elements with odd degree k 𝑘 k for the free Lie algebra generated by operators { A 1 , A 3 , A 5 , … } subscript 𝐴 1 subscript 𝐴 3 subscript 𝐴 5 … \{A_{1},A_{3},A_{5},\dotsc\} , where A k subscript 𝐴 𝑘 A_{k} has degree k 𝑘 k .

Degree $k$	1	3	5	7	9	11	13
$d_{\text{L},k}$	1	1	2	4	8	18	40

Table 5. Table 6: Hall basis elements of odd degree for the free Lie algebra generated by operators { A 1 , A 2 , A 3 , A 4 , … } subscript 𝐴 1 subscript 𝐴 2 subscript 𝐴 3 subscript 𝐴 4 … \{A_{1},A_{2},A_{3},A_{4},\dotsc\} , where A k subscript 𝐴 𝑘 A_{k} has degree k 𝑘 k .

Degree	Hall basis elements
1	$A_1$
3	$A_3$,  $[A_1,A_2]$
5	$A_5$,  $[A_1,A_4]$,  $[A_2,A_3]$,  $[A_1,[A_1,A_3]]$,  $[A_2,[A_1,A_2]]$,  $[A_1,[A_1,[A_1,A_2]]]$
7	$A_7$,  $[A_1,A_6]$,  $[A_2,A_5]$,  $[A_3,A_4]$,  $[A_1,[A_1,A_5]]$,  $[A_2,[A_1,A_4]]$,
	$[A_2,[A_2,A_3]]$,  $[A_3,[A_1,A_3]]$,  $[A_4,[A_1,A_2]]$,  $[A_1,[A_1,[A_1,A_4]]]$,
	$[A_2,[A_1,[A_1,A_3]]]$,  $[A_2,[A_2,[A_1,A_2]]]$,  $[A_3,[A_1,[A_1,A_2]]]$,  $[[A_1,A_2],[A_1,A_3]]$,
	$[A_1,[A_1,[A_1,[A_1,A_3]]]]$,  $[A_2,[A_1,[A_1,[A_1,A_2]]]]$,  $[[A_1,A_2],[A_1,[A_1,A_2]]]$,
	$[A_1,[A_1,[A_1,[A_1,[A_1,A_2]]]]]$

Table 6. Table 7: Numbers d E , k subscript 𝑑 E 𝑘 d_{{\text{E}},k} of Hall bases elements with odd degree k 𝑘 k for the free Lie algebra generated by operators { A 1 , A 2 , A 3 , A 4 , … } subscript 𝐴 1 subscript 𝐴 2 subscript 𝐴 3 subscript 𝐴 4 … \{A_{1},A_{2},A_{3},A_{4},\dotsc\} , where A k subscript 𝐴 𝑘 A_{k} has degree k 𝑘 k .

Degree $k$	1	3	5	7	9	11	13
$d_{\text{E},k}$	1	2	6	18	56	186	630

Equations

$$e^{tH}=e^{x_{1}tA}{e}^{x_{2}tB}{e}^{x_{3}tA}{e}^{x_{4}tB}\cdots e^{x_{m}tA}+\mathcal{O}(t^{p+1})$$

$$U_{\text{N},m}:=\{\begin{array}{ll}{e}^{a_{1}tA}{e}^{b_{1}tB}{e}^{a_{2}tA}{e}^{b_{2}tB}\cdots e^{b_{m/2}tB}& \text{for even m},\\ e^{a_{1}tA}{e}^{b_{1}tB}{e}^{a_{2}tA}{e}^{b_{2}tB}\cdots e^{a_{(m+1)/2}tA}& \text{for odd m}\end{array}$$

$$U_{\text{S},m}:=\{\begin{array}{ll}{e}^{a_{1}tA}{e}^{b_{1}tB}{e}^{a_{2}tA}\cdots e^{b_{\nu -1}tB}{e}^{a_{\nu }tA}{e}^{b_{\nu -1}tB}\cdots e^{a_{1}tA}& \text{for odd ν},\\ e^{a_{1}tA}{e}^{b_{1}tB}{e}^{a_{2}tA}\cdots e^{a_{\nu -1}tA}{e}^{b_{\nu }tB}{e}^{a_{\nu -1}tA}\cdots e^{a_{1}tA}& \text{for even ν}\end{array}$$

$$U_{\text{L}}(\tau ):=e^{\frac{1}{2}\tau tA}{e}^{\tau B}{e}^{\frac{1}{2}\tau A}$$

$$U_{\text{SL},m}$$

$$=e^{\frac{1}{2}{w}_{1}tA}{e}^{w_{1}tB}{e}^{\frac{1}{2}(w_1+w_2)tA}{e}^{w_{2}tB}{e}^{\frac{1}{2}(w_2+w_3)tA}\cdots e^{w_{\nu }tB}\cdots e^{\frac{1}{2}{w}_{1}tA}$$

$$U=e^{tH+\mathcal{O}(t^{p+1})}=e^{tH}+\mathcal{O}(t^{p+1}).$$

$$\text{SL}\subseteq \text{S}\subseteq \text{N}.$$

$$\log \left(e^X{e}^Y\right)=X+Y+\frac{1}{2}[X,Y]+\frac{1}{12}[X,[X,Y]]-\frac{1}{12}[Y,[X,Y]]-\frac{1}{24}[Y,[X,[X,Y]]]\dots$$

$$Z=tH+\mathcal{O}(t^{p+1})\text{with}Z:=\log U.$$

$$X_1=tA\text{and}{Y}_1=tB\text{to get}{Z}_1=\log (e^{tA}{e}^{tB}).$$

$$X_1=Z_1\text{and}{Y}_2=tC\text{to get}{Z}_2=\log (e^{Z_1}{e}^{tC})=\log (e^{tA}{e}^{tB}{e}^{tC})$$

$$[V,[W,X]]+[W,[X,V]]+[X,[V,W]]=0.$$

$$d_k=\frac{1}{j}\sum _{j|k}\mu (j)\cdot n^{k/j},$$

$$\log U_{\text{L}}(\tau )=\tau Z_{\text{L}}^{(1)}+{\tau }^3{Z}_{\text{L}}^{(3)}+{\tau }^5{Z}_{\text{L}}^{(5)}+\dots \text{with}{Z}_{\text{L}}^{(1)}=H.$$

$$\log U=tH+t^{p+1}{\sum }_i{c}_i{X}_i+\mathcal{O}(t^{p+2})$$

$$\Vert \log U-tH\Vert =t^{p+1}\Vert {\sum }_i{c}_i{X}_i\Vert +\mathcal{O}(t^{p+2})\le t^{p+1}{\sum }_i|c_i|\Vert X_i\Vert +\mathcal{O}(t^{p+2}).$$

$$U=e^{tH}+t^{p+1}{\sum }_i{c}_i{X}_i+\mathcal{O}(t^{p+2})\text{and}\Vert U-e^{tH}\Vert \le t^{p+1}{\sum }_i|c_i|\Vert X_i\Vert +\mathcal{O}(t^{p+2}).$$

$$\epsilon :={\left(\frac{m}{p}\right)}^p\min \left({\sum }_i|c_i^{AB}|,{\sum }_i|c_i^{BA}|\right)$$

$$U_{\text{L}}(wt)U_{\text{L}}((1-2w)t)U_{\text{L}}(wt)=e^{tH}+\mathcal{O}(t^5)\text{with}w=(2-2^{1/3}{)}^{-1}\approx 1.35121$$

$$\begin{array}{rl}& U_{\text{Y},q+1}(t):=U_{\text{Y},q}(y_{q}t)U_{\text{Y},q}((1-2y_q)t)U_{\text{Y},q}(y_{q}t)=e^{tH}+\mathcal{O}(t^{2q+3})\\ & \text{with}{y}_q={\left(2-2^{1/(2q+1)}\right)}^{-1}.\end{array}$$

$$\begin{array}{rl}& U_{\text{Z},q+1}(t):=U_{\text{Z},q}(z_{q}t)U_{\text{Z},q}(z_{q}t)U_{\text{Z},q}((1-4z_q)t)U_{\text{Z},q}(z_{q}t)U_{\text{Z},q}(z_{q}t)=e^{tH}+\mathcal{O}(t^{2q+3})\\ & \text{with}{z}_q={\left(4-4^{1/(2q+1)}\right)}^{-1}.\end{array}$$

$$b_1=\frac{1}{2},a_2=1-2a_1,a_1=\frac{1}{12y}(y^2+6y-2)\approx 0.19318\text{with}y=(2\sqrt{326}-36{)}^{1/3}$$

$$b_1=\frac{1}{2},a_2=1-2a_1,a_1=\frac{1}{6}(3-\sqrt{3})\approx 0.21132$$

$$\begin{array}{rl}& b_2=\frac{1}{2}-b_1,a_3=1-2(a_1+a_2),b_1=\frac{6}{11}\approx 0.54545,\\ & a_1=\frac{1}{3924}\left(642+\sqrt{471}\right)\approx 0.16914,a_2=\frac{121}{3924}\left(12-\sqrt{471}\right)\approx -0.29919\end{array}$$

$$a_1=0.1720865590295143,b_1=0.5915620307551568,a_2=-0.1616217622107222$$

$$\begin{array}{rl}& b_2=\frac{1}{2}-b_1,a_3=1-2(a_1+a_2),b_1=-0.35905925216967795307,\\ & a_1=0.26756486526206148829,a_2=-0.034180403245134195595\end{array}$$

$$\begin{array}{rl}& b_2=\frac{1}{2}-b_1,a_3=1-2(a_1+a_2),b_1=-\frac{1}{3},\\ & a_1=\frac{17}{2}-\frac{5}{2}\sqrt{\frac{65}{6}},a_2=\frac{3}{20}\left(\sqrt{390}-20\right).\end{array}$$

$$\begin{array}{rl}& b_2=\frac{1}{2}-b_1,a_3=1-2(a_1+a_2),b_1=0.60417497648530223585,\\ & a_1=0.17285948240376668244,a_2=-\mathrm{0.14265971252922336963.}\end{array}$$

$$w_3=1-2(w_1+w_2),w_1=\frac{3\pm \sqrt{3}}{6},w_2=\frac{3\mp \sqrt{3}}{6}$$

$$w_3=1-2(w_1+w_2),w_1=0.28,w_2=0.62546642846767004501$$

$$w_3=1-2(w_1+w_2),w_1=0.3221375960817984,w_2=0.5413165481700430$$

$$w_3=1-2(w_1+w_2),w_1=0.25686635900587695859,w_2=0.67762403230558747362$$

Full text

Lie-Trotter-Suzuki decomposition, Lie-Trotter product formula, operator exponential, tensor network states, quantum Hamiltonian simulation, quantum Monte Carlo, symplectic integrators, splitting methods, Baker-Campbell-Hausdorff formula, Hall basis, Gröbner basis, graph partitioning problem

1. I Introduction

2. II Coarse-graining and partitioning of interaction graphs

3. III Lie-Trotter-Suzuki decompositions for $n=2$ terms and general aspects

4. III.1 Considered types of decompositions

5. III.2 Baker-Campbell-Hausdorff formula

6. III.3 Free Lie algebra and Hall bases

7. III.4 Numbers of parameters and constraints, symmetries

8. III.5 Ideal of constraint polynomials and Gröbner bases

9. III.6 Relevance of 1-norm as an error measure

10. III.7 Unoptimized decompositions due to Forest, Ruth, Yoshida, and Suzuki

11. IV Optimized decompositions for $n=2$ terms

12. IV.1 Order $p=2$

13. IV.2 Order $p=4$

14. IV.3 Order $p=6$

15. V Lie-Trotter-Suzuki decompositions for $n=3$ terms

16. V.1 Considered types of decompositions

17. V.2 Numbers of parameters and constraints, symmetries

18. VI Optimized decompositions for $n=3$ terms

19. VI.1 Order $p=1$

20. VI.2 Order $p=2$

21. VI.3 Order $p=4$

22. VI.4 Order $p=6$

23. VII Discussion

I Introduction

In many situations, we need to evaluate or apply operator exponentials $e^{tH}$ where $H$ acts in a huge vector space. A first trick is to decompose the task into small time steps, $e^{(t_{2}-t_{1})H}=(e^{tH})^{N}$ with $Nt=t_{2}-t_{1}$, and we will always consider $t$ to denote that time step in the following. While $H$ refers to a Hamiltonian in many applications, we do not assume $H$ to be Hermitian or skew-Hermitian. For example, it could also be a Liouville super-operator for a Markovian open quantum system. When $H$ is a sum of $n$ (non-commuting) terms which, individually, can be exponentiated easily, one can approximate $e^{tH}$ by a product of these easy exponentials Trotter1959 ; Suzuki1976-51 . This gives Lie-Trotter-Suzuki decompositions like

[TABLE]

for the case of a generator $H=A+B$ with $n=2$ terms. Depending on the number of factors $m$ in the product approximation (1), one can achieve different approximation orders $p$ with respect to $t$ and may still have free parameters to minimize the approximation error. In particular, when there are free parameters left, we can use them to minimize the amplitude of the leading order error term $\propto t^{p+1}$. We discuss and present such optimized Lie-Trotter-Suzuki decompositions for $n=2$ and $n=3$ terms in the generator $H$.

Lie-Trotter-Suzuki decomposition have many applications. For example, they are employed in time evolution algorithms for matrix product states Vidal2003-10 ; White2004 ; Daley2004 ; Orus2008-78 or higher-dimensional tensor network states called projected entangled pair states Verstraete2004-7 ; Verstraete2006-96 ; Niggemann1997-104 ; Nishino2000-575 ; Martin-Delgado2001-64 . Importantly, Lie-Trotter-Suzuki decompositions can also be used for digital quantum simulation of many-body systems on quantum computers Lloyd1996-273 ; Barry2007-270 ; Lanyon2011-334 ; Kliesch2011-107 ; Barreiro2011-470 . Furthermore, they are important tools in path integral methods like worldline quantum Monte Carlo Suzuki1977-58 , diffusion Monte Carlo Kolorenc2011-74 , and for approximate symplectic integrators for the Hamilton equations of classical systems Ruth1983-30 ; Yoshida1990 ; McLachlan1995 . In the latter context, the decompositions are also called splitting or composition methods.

Here, we provide optimized decompositions up to order $p=6$ for $n=2$ and $n=3$ terms in the generator $H$. We minimize a bound on the operator norm of the leading error term by expanding in terms of nested commutators (Hall bases) and optimizing the 1-norm of the coefficients. In the comparison of different decompositions, we take into account that increased numbers $m$ of factors in the decomposition need to be compensated by correspondingly larger time steps $t$ to keep computation costs constant. For $n=2$ terms, with a few exceptions, we find optima that are close to those in Ref. McLachlan1995 , where McLachlan minimized an error measure suitable for classical Hamiltonian systems. Generally, our results improve substantially over unoptimized decompositions due to Forest and Ruth, Yoshida, and Suzuki Forest1990-43 ; Yoshida1990 ; Suzuki1990 . We also compare to decompositions of Kahan and Li, and Omelyan et al. Kahan1997-66 ; Omelyan2002-146 . The optimized decompositions for $n\in\{2,3\}$ can be used to efficiently simulate any one-dimensional (1d) and two-dimensional lattice model with finite-range interactions (geometrically local Hamiltonians) as we explain with a coarse-graining argument. In fact, several of the considered decompositions (those of type SL and SE) are generally applicable for any number of terms $n$, but they are in general not optimal when applied for $n\geq 4$.

In section II, we show that interaction graphs of 1d and 2d lattice models with finite-range interactions can always be partitioned into $n\leq 3$ subgraphs $A$, $B$ (and $C$), each corresponding to sums of local operators with disjoint spatial supports. In section III, we discuss Lie-Trotter-Suzuki decompositions for $n=2$ terms and further general aspects: the considered types of decompositions (Sec. III.1), the Baker-Campbell-Hausdorff formula which is essential to expand decompositions in terms of nested commutators (Sec. III.2), Hall bases for the free Lie algebra generated by the terms in $H$ (used to remove linear dependence of nested commutators; see Sec. III.3), the numbers of parameters and constraints in the different types of $n=2$ decompositions (Sec. III.4), Gröbner bases which can be used to solve the systems of polynomial constraints and to identify free parameters (Sec. III.5), the 1-norm error measure which bounds the operator-norm approximation accuracy and is minimized in optimizations as well as the issue of quasi-locality (Sec. III.6), and generic unoptimized decompositions due to Forest, Ruth, Yoshida, and Suzuki (Sec. III.7). In section IV, we provide various optimized decompositions for $n=2$, compare to previous results, and make recommendations. Section V describes the considered types of decompositions for $n=3$ terms (Sec. V.1) as well as their numbers of parameters and constraints (Sec. V.2). In section VI, we provide various optimized decompositions for $n=3$, compare to previous results, and make recommendations. We close with a short discussion of the results and alternatives to Lie-Trotter-Suzuki decompositions (Sec. VII).

II Coarse-graining and partitioning of interaction graphs

The first step in the design of Lie-Trotter-Suzuki decompositions for $e^{tH}$ is to partition $H$ into a sum of $n$ terms $A,B,\dots$ such that each of these terms can be exponentiated easily, i.e., with low cost. For classical systems, these are typically the kinetic energy term and the potential energy term. In this section, we consider quantum many-body lattice models with finite-range 2-local interactions. As explained in the following, we can always find partitionings into $n\leq 2$ terms for such 1d systems and $n\leq 3$ terms for 2d systems. In particular, every term will be a sum of mutually commuting operators with finite (local) spatial support. As all elementary interactions are two-particle interactions, 2-locality is natural and, in solid state systems, lattice models with finite-range interactions naturally arise due to screening effects and the tight-binding approximation Mahan2000 .

Using coarse graining as in Kadanoff’s block spin transformation, we can always reduce the problem to nearest neighbor interactions:

Lemma 1** (Reduction to nearest neighbor interactions).**

Consider lattice models with finite-range 2-local interactions. In 1d, a finite number of coarse-graining steps is sufficient to map to a model with onsite and nearest neighbor interactions only. In 2d, a finite number of coarse-graining steps is sufficient to map to a model on the triangular lattice with onsite and nearest neighbor interactions only.

The situation for 1d is displayed in Fig. 1(a). When we coarse grain two sites into one new effective site, an interaction between sites of distance $\Delta x$ maps to an interaction between effective sites of distance $\operatorname{floor}(\Delta x/2)$ or $\operatorname{floor}(\Delta x/2)+1$. After a finite number of coarse-graining steps, we arrive at a model with interactions of range $\Delta x\in\{0,1\}$. The situation for 2d is displayed in Fig. 1(b-d). In a first step, any 2d lattice can be deformed into a square lattice. When we then coarse grain a square of four sites into one new effective site, an interaction between sites of distance $(\Delta x,\Delta y)$ maps to an interaction between effective sites with $x$ distance $\operatorname{floor}(\Delta x/2)$ or $\operatorname{floor}(\Delta x/2)+1$ and $y$ distance $\operatorname{floor}(\Delta y/2)$ or $\operatorname{floor}(\Delta y/2)+1$. After a finite number of coarse-graining steps, we arrive at a square lattice with interactions of range $\Delta x,\Delta y\in\{0,1\}$, i.e., nearest and next-nearest neighbor interactions as shown in Fig. 1(c). In a final coarse graining step, we can map two neighboring sites to one effective site to arrive at a triangular lattice with onsite and nearest neighbor interactions only [Fig. 1(c,d)].

Lemma 2** (Partitioning 1d and 2d interaction graphs).**

Consider lattice models with finite-range 2-local interactions. The interactions for any 1d model can be partitioned into $n=2$ terms. The interactions for any 2d model can be partitioned into $n=3$ terms. Each of the resulting terms is a sum of local interaction operators with disjoint spatial supports.

According to Lemma 1, we can use coarse graining to reduce the problem to nearest neighbor interactions on 1d lattices and 2d triangular lattices. For 1d lattices with nearest neighbor interactions, the interaction graph can be partitioned into two terms, the first containing interactions on even bonds (of effective sites after coarse graining) and the second containing interactions on odd bonds as shown in Fig. 2(a). For 2d triangular lattices with nearest neighbor interactions, the interaction graph can be partitioned into three terms. All interactions are grouped into three-site operators acting on triangular plaquettes. As indicated in Fig. 2(e), three terms ($A,B,C$) are required such that the three-site operators in each term are all mutually commuting.

The described coarse-graining to achieve a partitioning into $n=2$ or $3$ terms, is always possible in 1d and 2d but may not always be the best choice. For instance, in a 1d model with nearest and next-nearest neighbor interactions, we would coarse grain once [Fig. 1(a)] to arrive at a partitioning with $n=2$. The resulting nearest neighbor interactions between effective sites correspond to four-site operators for the original lattice. Depending on the specific classical or quantum computation costs, it may be preferable to partition into three terms as in Fig. 2(b) and use an $n=3$ Lie-Trotter-Suzuki decomposition instead of a decomposition for $n=2$. Another example is the square 2d lattice with nearest neighbor interactions. One can partition into $n=2$ terms with four-site operators [Fig. 2(c)] or into $n=3$ terms with three-site operators [Fig. 2(d)].

III Lie-Trotter-Suzuki decompositions for $n=2$ terms and general aspects

For $n=2$, the generator $H=A+B$ consists of two terms, where exponentials $e^{tA}$ and $e^{tB}$ can be computed easily. For example, both can be sums of commuting few-particle interaction terms. Examples for lattice models are given in Fig. 2(a,c,g). Keep in mind that we actually do not assume $H$ to be Hermitian such that our decompositions apply equally well to real- and imaginary-time evolution, the evolution under Liouvillians, and quite generally, for any matrix exponentials.

III.1 Considered types of decompositions

We consider the following types of Lie-Trotter-Suzuki decompositions for $e^{tH}$, where $m$ denotes the total number of operator exponentials:

• •

Type N with $\nu=m$ parameters. This is the most generic type of decompositions

[TABLE]

with $a_{i},b_{i}\in\mathbb{R}$.

• •

Type S with $\nu=(m+1)/2$ parameters and odd $m$. This is the most generic type of symmetric decompositions

[TABLE]

with $a_{i},b_{i}\in\mathbb{R}$.

• •

Type SL with $\nu=\operatorname{ceil}((m-1)/4)$ parameters and odd $m$. This type of symmetric decomposition is a product of “leapfrog” terms 111The name alludes to the similarity to the leapfrog integration method for differential equations. This second order decomposition is also known as the Verlet integrator.

[TABLE]

such that

[TABLE]

with $w_{i}\in\mathbb{R}$. Note that, as shown in the last line, exponentials of $A$ for subsequent leapfrog terms can be contracted into one such that the total number of required operator exponentials is indeed $m$.

Depending on the number $m$ of factors in a decomposition $U$ of type (2)-(5), the parameters can be chosen such that $U$ coincides with the exact $e^{tH}$ up to order $p$ in the sense that

[TABLE]

Depending on $m$, the number of parameters in a decomposition can be larger than the number of constraints due to Eq. (6). We can then use the remaining freedom to minimize the (leading-order) error term of the decomposition in a suitable metric.

Of course, for the same number of factors $m$, type SL is a subclass of type S, and type S is a subclass of type N,

[TABLE]

At first sight, one might consider types S and SL to be superfluous. But we will see in Sec. III.4 that types S and SL do not only have smaller numbers of parameters $\nu$, but also fewer constraints than a generic type-N decomposition for the same order $p$. For certain $m$ and $p$, type-N decompositions solving the constraint (6) may then necessarily have to be of type S or such of type S may have to be of type SL. In this sense, below a certain minimum $m$ only decompositions of type S (or SL) may exist but none of type (SL or) N.

III.2 Baker-Campbell-Hausdorff formula

How can we compare a Lie-Trotter-Suzuki decomposition $U$ and $e^{tH}$ as suggested by Eq. (6)? One option would be to study the Taylor expansions of both operators in terms of operator monomials like $A^{2}BA$, require the difference of the expansions to vanish up to order $t^{p}$ and minimize some metric for the scalar coefficients of operator monomials in the $t^{p+1}$ term. However, in typical applications of many-body physics with finite-range interactions in $H$, norms for operator monomials of order $p+1$ scale as $\mathcal{O}(L^{p+1})$ with the system size $L$ which results in very large error bounds (cf. Sec. III.6).

We can use the Baker-Campbell-Hausdorff (BCH) formula Campbell1897-29 ; Baker1905-s2 ; Hausdorff1906 ; Dynkin1947-57 ; Dynkin1949 ; Varadarajan1984

[TABLE]

to resolve this issue. In particular, we can apply it recursively, to expand $\log U$ for the Lie-Trotter-Suzuki decompositions in terms of nested commutators of $A$ and $B$. Then, the constraints (6) are equivalent to

[TABLE]

In applications with finite-range interactions in $H$, norms (norm bounds) for the nested commutators in $Z$ will all scale linearly with the system size $L$.

To obtain $Z=\log U$ for a decomposition $U=e^{tA}\,e^{tB}\,e^{tC}\,e^{tD}\dotsb$, one can first use the BCH formula with

[TABLE]

and so on until we arrive at $Z=\log U$.

III.3 Free Lie algebra and Hall bases

Campbell, Baker, and Hausdorff Campbell1897-29 ; Baker1905-s2 ; Hausdorff1906 found that $\log\left(e^{X}e^{Y}\right)$ can be expressed in terms of nested commutators, i.e, that it is an element of the Lie algebra generated by $X$ and $Y$ Rossmann2002 . Dynkin finally derived an explicit formula Dynkin1947-57 ; Dynkin1949 .

The iteration (10) of the BCH formula gives an expansion of $Z=\log U$ in terms of nested commutators of $A$ and $B$. The number of constraints, imposed by requiring $U$ to coincide with $e^{tH}$ up to order $t^{p}$, can be determined from the number of terms in the expansion. In general, the nested commutators are not linearly independent, especially, because of the Jacobi identity

[TABLE]

Let us discuss this problem more generally for $n$ operators instead of just two ($A$ and $B$). To resolve the above problem, we need a basis for the Lie algebra generated by operators $A,B,C,\dots$. For a free Lie algebra, we can use Hall bases Hall1950-1 ; Serre1992 ; Reutenauer1993 . A free Lie algebra is fully characterized by the properties of the commutator from which all algebraic relations like the Jacobi identity (11) follow and no further relations exit. For a finite vector space, the Lie algebra generated by $n$ operators will of course close at some point. However, for the large Hilbert spaces relevant in many-body physics and low expansion orders considered for our optimized decompositions, the assumption of a free Lie algebra is generally sufficient as we will see below.

For the construction of a Hall basis, one introduces a total order “$<$” on the generators and nested commutators, e.g.: $A<B<C<\dots$ for the generators, $X<Y$ for two nested commutators if $X$ is of lower degree than $Y$ and $[X,Y]<[V,W]$ if $X<V$, or $X=V$ and $Y<W$. According to Witt’s formula Witt1956-64 ; Reutenauer1993 , the number of Hall basis elements of degree $k$ is given by the necklace polynomial

[TABLE]

where $\mu$ denotes the Möbius function and the sum is over all integers $j$ that divide $k$. For $n=2$ and $3$ generators, Tables 1 and 2 list Hall basis elements and their numbers $d_{k}$. Note that Witt’s formula (12) addresses the case where all generators have degree one. For the decomposition type SL [Eq. (5)] and further decomposition for $n=3$ terms in Sec. V.1, we will also consider generators of higher degree.

III.4 Numbers of parameters and constraints, symmetries

Assuming that all terms that can occur do occur in the expansion of $Z=\log U$ for an order-$p$ Lie-Trotter-Suzuki decomposition $U$, we can determine the number of constraints from the number of relevant Hall basis elements. For example, a type-N decomposition (2) with $m$ factors has $m$ parameters $\{a_{i},b_{i}\}$ and the constraints due to Eq. (9) are the polynomials in $\{a_{i},b_{i}\}$ that appear in $Z$ as coefficients of Hall basis elements with degree $\leq p$. The numbers $\sum_{k=1}^{p}d_{k}$ of these are given in Eq. (12) and Table 2.

After taking into account symmetries etc., we will find that indeed all relevant Hall basis elements occur for the considered decompositions and expansion orders except for some obvious cases that are due to the fact that type SL decompositions are a subclass of those of type S, and the latter a subclass of the type N decompositions. To check and determine how many free parameters we actually have, after the expansion in a Hall basis, we can employ the Gröbner basis of the constraint polynomials as discussed in the following section (Sec. III.5).

Type N. – For the decompositions $U_{\text{N},m}$ in Eq. (2), we have $\nu=m$ parameters $\{a_{i},b_{i}\}$. As the decomposition has no further symmetry or structure, all Hall basis elements should occur in the expansion of $Z$ and the number of constraints to achieve approximation order $p$ is $\sum_{k=1}^{p}d_{k}$. For example, at first order, we have the two constraints $\sum_{i}a_{i}=\sum_{i}b_{i}=1$ to achieve $Z=tH+\mathcal{O}(t^{2})$.

Type S. – For the decompositions $U_{\text{S},m}$ in Eq. (3), we have $\nu=(m+1)/2$ parameters and odd $m$. Due to the symmetry in the factors, the decompositions obey the time reversal symmetry $U_{\text{S},m}(t)U_{\text{S},m}(-t)=\mathbbm{1}$. It follows that $Z(t)$ only contains terms of odd order in $t$ such that no Hall basis elements of even degree occur Yoshida1990 :

Lemma 3** (Time reversal symmetry and order).**

Let $U(t)$ be of the form $U(t)=e^{tA}\,e^{tB}\,e^{tC}\dotsc$ and obey time reversal symmetry $U(t)U(-t)=\mathbbm{1}$. In an expansion $\log U(t)=tZ^{(1)}+t^{2}Z^{(2)}+t^{3}Z^{(3)}+t^{4}Z^{(4)}+\dots$, all even-order terms vanish, i.e., $Z^{(2)}=Z^{(4)}=Z^{(6)}=\dots=0$, and $Z^{(1)}=A+B+C+\dots$

As an immediate consequence, a symmetric order $p$ decomposition (obeys $Z^{(1)}=H$ and $Z^{(2)}=Z^{(3)}=\dots=Z^{(p)}=0$) with odd $p$ is actually also an order $p+1$ decomposition Yoshida1990 ; Suzuki1990 . For the proof of Lemma 3, note that $\log U(-t)=-tZ^{(1)}+t^{2}Z^{(2)}-t^{3}Z^{(3)}+t^{4}Z^{(4)}+\dots$ Applying the BCH formula (8), we find $\log[U(t)U(-t)]=2t^{2}Z^{(2)}+\mathcal{O}(t^{3})$ such that $Z^{(2)}=0$ because of the constraint $U(t)U(-t)=\mathbbm{1}$. With this information, we can reconsider the BCH formula and find $\log[U(t)U(-t)]=2t^{4}Z^{(4)}+\mathcal{O}(t^{5})$, showing that $Z^{(4)}=0$. Continuing in this way, we find that all even-order terms vanish.

For a symmetric decomposition that has no further structure, all Hall basis elements of odd degree should occur in the expansion of $Z$ and the number of constraints to achieve approximation order $p$ is $\sum_{q,2q-1\leq p}d_{2q-1}$. For example, at first order, we have the two constraints that the coefficients $a_{i}$ in the $A$ factors and the coefficients $b_{i}$ in the $B$ factors in Eq. (3) sum to one to achieve $Z=tH+\mathcal{O}(t^{2})$.

Type SL. – For the decompositions $U_{\text{SL},m}$ in Eq. (5), we have $\nu=\operatorname{ceil}((m-1)/4)$ parameters and odd $m$. As these decompositions are symmetric, according to Lemma 3 only terms of odd degree occur in the expansion of $Z=\log U_{\text{SL},m}$. The decomposition is a product of leapfrog terms $U_{\text{L}}(\tau)$ [Eq. (4)] which are symmetric second order decompositions of $e^{\tau H}$ and can hence be expanded in the form

[TABLE]

$Z^{(k)}_{\text{L}}$ is a Lie polynomial containing only nested commutators of degree $k$. To determine the number of constraints for $U_{\text{SL},m}$ to be an order-$p$ decomposition of $e^{tH}$, we can consider the Lie algebra generated by $\{Z^{(1)}_{\text{L}},Z^{(3)}_{\text{L}},Z^{(5)}_{\text{L}},\dots\}$. Again, assuming no further relevant algebraic relations for the generators, we can treat the problem as that of a free Lie algebra. The number of constraints due to Eq. (9) is then given by the number $\sum_{q,2q-1\leq p}d_{\text{L},2q-1}$ of Hall basis elements of odd degree $2q-1\leq p$. Here, for example, the degree of $[Z^{(1)}_{\text{L}},[Z^{(1)}_{\text{L}},Z^{(5)}_{\text{L}}]]$ is 7. Tables 3 and 4 list the Hall basis elements and their numbers $d_{\text{L},k}$. At first order, we only have the constraint that the coefficients $w_{i}$ of all leapfrog factors in Eq. (5) sum to one to achieve $Z=tH+\mathcal{O}(t^{2})$.

To summarize this section, we give the total numbers of constraints by order in Table 5. The table also states $m_{p}^{\text{H}}$, the minimum number of factors needed in the different decompositions to obtain an order-$p$ decomposition [cf. Eq. (6)], under the assumption that all constraint polynomials are independent. This is of course not true in all cases. For example, for sixth order decompositions we have $m_{p}^{\text{H}}=15,19,23$ for types SL, S, and N, respectively. But still, any decomposition of type SL (S) is at the same time of type S (N). How to reconcile these statements? If we want to construct sixth order decompositions of type S or N with $m<19$, we can do so, but they turn out to be of the specific type SL, i.e., consist of leapfrog factors. And if we want to construct a type-N decomposition with $19\leq m<23$, we can do so, but they turn out to be of the specific type S, i.e., symmetric.

Although we have not encountered this situation in our current search for optimized decomposition, it could also occur, that the constraint polynomials (the coefficients of the different Hall basis elements discussed above) are not independent. For example, a polynomial could be a multiple of another one.

III.5 Ideal of constraint polynomials and Gröbner bases

For decompositions $U$ as defined in Sec. III.1, the constraints to obtain a decomposition of order $p$ [Eq. (9)] are polynomials in the $\nu$ parameters of the decomposition ($\{a_{i},b_{i}\}$ or $\{w_{i}\}$) and can be read off as the coefficients of Hall basis elements in an expansion of $Z=\log U$. If the constraints are all independent, the number of free parameters is simply given by $\nu-N^{\text{H}}_{p}$. Here, $N^{\text{H}}_{p}$ denotes the number of Hall basis elements with degree $\leq p$ in an expansion of $Z$ (cf. Table 5).

The independence of the constraints can be checked by constructing a Gröbner basis for the ideal defined by the constraint polynomials 222For a finite set $\mathcal{C}=\{c_{1},\cdots,c_{k}\}$ of polynomials for variables ${\bm{x}}$, the ideal generated by $\mathcal{C}$ is the set of linear combinations of the $c_{i}$ with arbitrary polynomial coefficients.. The Gröbner basis Buchberger1965 ; Adams1994 – a generating set for the ideal – can also be used to determine maximal sets of independent variables or the number of solutions if it is finite. Finally, the Gröbner basis can be used to find solutions for the system of polynomial equations. Algorithms to compute Gröbner bases were given by Buchberger Buchberger1976-10 ; Cox2007 and Faugère Faugere1999-139 ; Faugere2002 .

III.6 Relevance of 1-norm as an error measure

Above, we have discussed how to construct Lie-Trotter-Suzuki decompositions $U$ as products of operator exponentials to approximate $e^{tH}=e^{t(A+B)}$ to order $p$ [Eq. (6)]. The more factors we include in a decomposition, the more free parameters we have which allow us to reduce the amplitude of deviations from $e^{tH}$. We will only consider the leading error term which is of order $t^{p+1}$. Like McLachlan in Ref. McLachlan1995 we will use the free parameters in the decomposition to minimize a measure for the leading error term.

We can expand the leading error term in elements $\{X_{i}\}$ of degree $p+1$ from a Hall basis that is generated by $A$ and $B$ such that

[TABLE]

The relevant error measure is the operator norm distance which we can bound using the triangle inequality such that

[TABLE]

The $X_{i}$ are nested commutators of $A$ and $B$ with degree $p+1$. In our generic optimization we do not want to assume any detailed information about $A$ and $B$. To simplify matters further, we use a uniform norm bound for all $X_{i}$: For a typical situation where $A$ and $B$ are both sums of finite-range interaction terms with disjoint spatial supports, the number of terms in all $X_{i}$ will be linear in the system size $L$. However, the spatial support of a term $[A,[A,[A,[A,B]]]]$ would then be smaller than that of terms like $[B,[A,[B,[A,B]]]]$ and one could use these properties to further improve the error measure.

So, with a uniform norm bound on the $X_{i}$ we can use the 1-norm $\sum_{i}|c_{i}|$ of the expansion coefficients as an error measure. There is one more complication. The $c_{i}$ can and often do depend on the choice of the Hall basis. The latter depends on the order that is chosen for the generators $A$ and $B$. As both $A<B$ and $B<A$ are perfectly fine and result in a valid upper bound on the norm distance, we can minimize with respect to the ordering. Let the two corresponding sets of coefficient polynomials be denoted by $\{c^{AB}_{i}\}$ and $\{c^{BA}_{i}\}$, respectively. We will then use the error measure

[TABLE]

to quantify the magnitude of the deviation of $U$ from $e^{tH}$, where $m$ denotes again the number of factors in the decomposition $U$. Like the coefficient polynomials $c_{i}$, the error $\varepsilon$ is a function of the parameters in the decomposition.

Similar to Ref. McLachlan1995 , the prefactor $m^{p}$ in Eq. (16) has been chosen to allow for a fair comparison of decompositions with the same order $p$, but different numbers $m_{j}$ of factors: We want to compare the accuracies for evolving the system for a time $T$ at constant computation cost. The evolution can be accomplished by applying the decomposition $U$ with time step size $t$ for $T/t$ times. Consider two order-$p$ decompositions with $m_{1}$ and $m_{2}$ factors, respectively. Assuming that the computational cost for the implementation of the exponentials $e^{atA}$ and $e^{btB}$ is uniform or comparable, we should choose time steps of size $t_{j}\propto m_{j}$. This scaling of the time step, the prefactor $t^{p+1}$ of the leading error term in Eq. (15), and the number $T/t$ of time steps motivate the factor $m^{p}$ in Eq. (16). The additional factor $1/p^{p}$ is irrelevant and just added in order to prohibit $\varepsilon$ from increasing too much when increasing $p$. Note that $\varepsilon$ can not be used to compare decompositions of different order $p$.

We have discussed above that norm bounds for the nested commutators in Eqs. (15) are linear in the system size $L$ for the case of finite-range interactions. In fact this is overpessimistic in many situations, in particular, if we are only interested in the evolution of local quantities, i.e., observables which are a sum of operators with finite spatial support. For the unitary time evolution of closed systems and Markovian dynamics of open quantum systems, $L$ can be replaced by $(vT)^{d}$, where $d$ denotes the number of spatial dimensions, $T$ the maximum time for which we want to evolve, and $v$ is a Lieb-Robinson velocity Lieb1972-28 ; Poulin2010-104 ; Nachtergaele2011-552 ; Barthel2012-108b . This is due to quasi-locality Nachtergaele2007-12a ; Barthel2012-108b : In the Heisenberg picture, one can truncate the evolution outside a region of size $\sim(vT)^{d}$ around the spatial support of a given local observable.

III.7 Unoptimized decompositions due to Forest, Ruth, Yoshida, and Suzuki

Based on the symmetric second order (leapfrog) decomposition $U_{\text{L}}$ in Eq. (4), Forest and Ruth Forest1990-43 found the fourth order decomposition

[TABLE]

consisting of three leapfrog factors. It is hence of type SL with $m=7$ in Eq. (5).

Yoshida Yoshida1990 generalized the approach. Let us define $U_{\text{Y},1}(t):=U_{\text{L}}(t)$ for the leapfrog decomposition (4). Yoshida then showed how to obtain decompositions $U_{\text{Y},q}$ of arbitrary (even) order $p=2q$ through the recursion

[TABLE]

This approach has two drawbacks: (a) the number $m$ of factors in the decompositions is $2\cdot 3^{q-1}+1=3,7,19,55,163,487,1459,\dotsc$ and grows considerably faster than the theoretical minimum $m_{p}^{\text{H}}$ for type-SL decompositions given in Table 5, and (b) the parameters $y_{q}=1.35121,\,1.17467,\,1.11618,\dotsc$ are all larger than one. This results in large error values $\varepsilon$ as we will see below and the series of decompositions does not converge in the limit $q\to\infty$.

Suzuki Suzuki1990 resolved the latter issue with what he called fractal decompositions, using more factors than above. Let us define $U_{\text{Z},1}(t):=U_{\text{L}}(t)$ for the leapfrog decomposition. Suzuki then showed how to obtain decompositions $U_{\text{Z},q}$ of arbitrary (even) order $p=2q$ through the recursion

[TABLE]

Now, both $z_{q}=0.414491,\,0.373066,\,0.359585,\,\dotsc$ and $|1-4z_{q}|$ are smaller than one. Specifically, $\frac{1}{3}<z_{q}<\frac{1}{2}$ and $|1-4z_{q}|<\frac{2}{3}$ Suzuki1990 . Still, this approach has drawbacks: (a) the number $m$ of factors in the decompositions is $2\cdot 5^{q-1}+1=3,11,51,251,1251,6251,31251\dotsc$ and grows hence even faster than in Yoshida’s scheme, and (b) the resulting error values $\varepsilon$ are still much larger than those of the optimized decompositions discussed below.

IV Optimized decompositions for $n=2$ terms

Here, we optimize the decompositions of types N, S, SL defined in Sec. III.1 with respect to their parameters to minimize the error measure $\varepsilon$ in Eq. (16). Keep in mind, that $\varepsilon$ is properly scaled to allow for a fair comparison of decompositions with the same order $p$ but different numbers of factors $m$; it takes into account that, for larger $m$, also the time step $t$ should be increased to compare integrators with the same computation costs.

We compare the results to the unoptimized type-SL decompositions of Forest and Ruth Forest1990-43 , Yoshida Yoshida1990 , and Suzuki Suzuki1990 . In most cases, we find that the optimal decompositions are relatively close to those obtained by McLachlan McLachlan1995 who used a different error measure. He expanded in a somewhat different basis suitable for symplectic integrators in classical Hamiltonian systems and, for practical reasons, used the 2-norm of the expansion coefficients instead of the more relevant 1-norm (16). The definition for McLachlan’s “Hamiltonian truncation error” is given in Ref. McLachlan1992-5 .

In the discussion below, $\nu$ denotes the number of parameters as stated in Sec. III.1 and the number of constraints for the different decomposition types and approximation orders $p$ are given in Table 5. We check the applicability of the corresponding counting argument using Gröbner bases for the constraint polynomials as discussed in Sec. III.5. For each order $p$, the recommended decomposition (usually smallest $\varepsilon$ found) is indicated by a star. When the two possible orders $A<B$ and $B<A$ chosen for the construction of the Hall basis are not equivalent, we specify in brackets the order that yields the minimum for $\varepsilon$.

IV.1 Order $p=2$

• Leapfrog ($m=3$, type SL, $\nu=1$). – The only parameter is fixed to $w_{1}=1$ by the constraint that the first order term in $\log U$ is $tH=t(A+B)$. The error is $\varepsilon=\left(\frac{3}{2}\right)^{2}\frac{1}{8}=0.28125$.

• McLachlan ($m=5$, type S, $\nu=3$). – There are two constraints and hence one free parameter. McLachlan states

[TABLE]

which has error $\varepsilon\approx 0.075192$.

$\bigstar$

Optimized ($m=5$, type S, $\nu=3$). – There are two constraints and hence one free parameter; we choose $a_{1}$. The error is minimized for

[TABLE]

which gives $\varepsilon\approx 0.069778$. At this point, the coefficient of the term $[A,[A,B]]$ vanishes.

$\blacksquare$

Discussion. – The best decomposition found here is of type S with $m=5$ as specified in Eq. (21). It improves over the common leapfrog decomposition by a factor $\sim 1/4$. According to Table 5, we can reach order $p=4$ with $m=7$ factors.

IV.2 Order $p=4$

• Forest & Ruth, Yoshida ($m=7$, type SL). – For the decomposition $U_{\text{Y},q=2}$ with $q=2$ in Eq. (18), the error is $\varepsilon\approx 0.38640$ ($A<B$).

• McLachlan, Omelyan et al. ($m=9$, type S, $\nu=5$). – There are four constraints and hence there is one free parameter. McLachlan states

[TABLE]

which has error $\varepsilon\approx 0.072483$ ($B<A$). Optimizing the coefficient 2-norm, Omelyan et al. Omelyan2002-146 find

[TABLE]

which gives $\varepsilon\approx 0.069248$ ($B<A$).

• Optimized ($m=9$, type S, $\nu=5$). – There are four constraints and hence there is one free parameter; according to the Gröbner basis, we can choose $b_{1}$. The error is minimized for

[TABLE]

• Suzuki, Kahan & Li ($m=11$, type SL). – For Suzuki’s decomposition $U_{\text{Z},q=2}$ in Eq. (19), the error is $\varepsilon\approx 0.216883$ ($A<B$). In Ref. Kahan1997-66 , Kahan and Li state the two solutions

[TABLE]

which have the error $\varepsilon\approx 0.17706$ ($A<B$).

• McLachlan, Omelyan et al. ($m=11$, type SL, $\nu=3$). – There are two constraints and hence there is one free parameter. McLachlan states

[TABLE]

which has error $\varepsilon\approx 0.11155$ ($A<B$). Optimizing the coefficient 2-norm, Omelyan et al. Omelyan2002-146 find

[TABLE]

which gives $\varepsilon\approx 0.13365$ ($A<B$).

• Optimized ($m=11$, type SL, $\nu=3$). – There are two constraints and hence there is one free parameter; according to the Gröbner basis, we can choose $w_{2}$. The error is minimized for

[TABLE]

• McLachlan ($m=11$, type S, $\nu=6$). – There are four constraints and hence two free parameters. McLachlan states

[TABLE]

which has error $\varepsilon\approx 0.023685$ ($A<B$).

$\bigstar$

Optimized ($m=11$, type S, $\nu=6$). – There are four constraints and hence two free parameters; according to the Gröbner basis, we can choose $\{b_{1},b_{2}\}$. The error is minimized for

[TABLE]

• Optimized ($m=13$, type SL, $\nu=3$). – There are two constraints and hence there is one free parameter; according to the Gröbner basis, we can choose $w_{2}$. The error is minimized for

[TABLE]

which gives $\varepsilon\approx 0.28728$ ($A<B$).

• Optimized ($m=13$, type S, $\nu=7$). – There are four constraints and hence three free parameters; according to the Gröbner basis, we can choose $\{a_{2},b_{2},a_{3}\}$. It is nontrivial to locate the global minimum. The best solution we found is

[TABLE]

$\blacksquare$

Discussion. – The best decomposition found here is of type S with $m=13$ as specified in Eq. (32). It improves over the decomposition $U_{\text{Y},q=2}$ [Eq. (18)] of Forest, Ruth, and Yoshida Forest1990-43 ; Yoshida1990 by a factor $\sim 1/28$ and by a factor of $\sim 1/16$ over the widely applied decomposition $U_{\text{Z},q=2}$ [Eq. (19)] due to Suzuki Suzuki1990 . Actually, the error of the decomposition (32) with $m=13$ factors does not improve too much over that of the type-S decomposition with $m=11$ factors [Eqs. (30a) and (30b)]. We have applied the latter in many tensor network simulations as in Refs. Barthel2013-15 ; Lake2013-111 ; Cai2013-111 ; Barthel2016-94 ; Barthel2017_08unused ; Binder2018-98 and recommend it generally for fourth order integration. It is not surprising that the optimized type-SL decomposition (31) with $m=13$ factors has a larger error than the SL decomposition (28a) with $m=11$ factors. The number of free parameters does not increase when going to $m=13$, but the increased number of factors is taken account of in the definition (16) of $\varepsilon$ and results in a larger error value. According to Table 5, we can reach order $p=6$ with $m=15$ factors.

IV.3 Order $p=6$

• Yoshida ($m=15$, type SL, $\nu=4$). – There are four constraints and hence no free parameters. According to the Gröbner basis, there are three real solutions which have already been determined numerically by Yoshida Yoshida1990 . Note that this is different from Yoshida’s generic solution (18) which has $m=19$ factors for $p=6$ and is discussed below. The best of the three solutions is

[TABLE]

which has error $\varepsilon\approx 0.44573$ ($A<B$). The other two solutions have much larger errors $\varepsilon\approx 5.7167$ and $\varepsilon\approx 5.8160$, respectively (both for order $A<B$).

• Yoshida, Kahan & Li ($m=19$, type SL). – For Yoshida’s decomposition $U_{\text{Y},q=3}$ with $q=3$ in Eq. (18), the error is $\varepsilon\approx 26.18692$ ($A<B$). In Ref. Kahan1997-66 , Kahan and Li state two similar solutions. The better of the two is

[TABLE]

with error $\varepsilon\approx 0.22167$ ($A<B$).

$\bigstar$

Optimized ($m=19$, type SL, $\nu=5$). – There are four constraints and hence there is one free parameter; according to the Gröbner basis, we can choose $w_{2}$. The error is minimized for

[TABLE]

which gives $\varepsilon\approx 0.17255$ ($B<A$). See Fig. 3(e). Albeit optimizing a different error measure, McLachlan gave almost the same decomposition in Ref. McLachlan1995 .

• Optimized ($m=23$, type SL, $\nu=6$). – There are four constraints and hence there are two free parameters; according to the Gröbner basis, we can choose $\{w_{4},w_{5}\}$. The best solution we found has

[TABLE]

which gives $\varepsilon\approx 0.17204$ ($A<B$). We identified two further good local minima with $\varepsilon\approx 0.19599$ ($B<A$) and $\varepsilon\approx 0.26551$ ($B<A$), respectively. See Fig. 3(f).

Suzuki ($m=51$, type SL). – For the decomposition $U_{\text{Z},q=3}$ in Eq. (19), the error is $\varepsilon\approx 0.84749$ ($B<A$). The order $A<B$ for the Hall basis would give a considerably larger error of $\approx 16.992$.

$\blacksquare$

Discussion. – The best decomposition found here is of type SL with $m=23$ as specified in Eq. (36). It improves over Yoshida’s generic decomposition $U_{\text{Y},q=3}$ [Eq. (18)] by a factor $\sim 1/152$, by a factor of $\sim 1/5$ over Suzuki’s decomposition $U_{\text{Z},q=3}$ [Eq. (19)], and by a factor of $\sim 2/5$ over the solution (33) that Yoshida found numerically Yoshida1990 . Actually, the error of the decomposition (36) with $m=23$ factors is almost identical to that of the decomposition (35) with $m=19$ factors. We hence recommend using the latter for sixth order integration. According to Table 5, true type-S decompositions require $m=19$ factors to reach order $p=6$. We have checked explicitly that type-S solutions with $m=15$ simply reproduce the corresponding type-SL decomposition (33). With $m=19$, type-S decompositions should just reproduce the type-SL decomposition (35) with one free parameter. With $m=23$, type-S decompositions have two free parameters like the corresponding type-SL decompositions (36). According to Table 5, there exist type-SL decompositions of order $p=8$ for $m\geq 31$ factors.

V Lie-Trotter-Suzuki decompositions for $n=3$ terms

For $n=3$, the generator $H=A+B+C$ consists of three terms, where exponentials $e^{tA}$, $e^{tB}$, $e^{tC}$ can be computed easily. Examples for corresponding partitionings of lattice model interaction graphs are given in Fig. 2(b,d,e,f). Much of the treatment for two terms in Sec. III carries over, but there are also some differences. In particular, we have two further relevant types of symmetric decompositions.

V.1 Considered types of decompositions

We consider the following types of Lie-Trotter-Suzuki decompositions for $e^{tH}$, where $m$ denotes the total number of operator exponentials:

• •

Type N with $\nu=m$ parameters. This is the most generic type of decompositions with the building block $ABCB$

[TABLE]

with $a_{i},b_{i},c_{i}\in\mathbb{R}$.

• •

Type S with $\nu=\operatorname{ceil}(m/2)=(m+1)/2$ parameters and even $(m-1)/2$. This is a symmetric decomposition with building block $ABCB$:

[TABLE]

with $a_{i},b_{i},c_{i}\in\mathbb{R}$.

• •

Type S-abc with $\nu=\operatorname{ceil}(m/2)$ parameters and even $(m+1)/3$. This is a symmetric decomposition with building block $ABC$:

[TABLE]

with $a_{i},b_{i},c_{i}\in\mathbb{R}$.

• •

Type SE with $\nu=(m-1)/4$ parameters and even $(m-1)/2$. This type of symmetric decomposition is a product of alternating two types of “Euler” terms 333The name alludes to the similarity to the Euler integration method for differential equations. This first order decomposition is the famous Lie-Trotter product formula Trotter1959 .

[TABLE]

such that

[TABLE]

with $u_{i},v_{i}\in\mathbb{R}$. Note that, as shown in the last line, exponentials of $C$ or $A$ for subsequent Euler terms $U_{\text{E}+}(ut)U_{\text{E}-}(vt)$ and $U_{\text{E}-}(vt)U_{\text{E}+}(ut)$, respectively, can be contracted into one such that the total number of required operator exponentials is indeed $m$. For $n=2$ terms, we did not discuss this type of decompositions because it is in that case simply equivalent to type S in Eq. (3) McLachlan1995 .

• •

Type SL with $\nu=\operatorname{ceil}((m-1)/8)$ parameters and integer $(m-1)/4$. This type of symmetric decomposition is a product of leapfrog terms

[TABLE]

such that

[TABLE]

with $w_{i}\in\mathbb{R}$. Note that, as shown in the last line, exponentials of $A$ for subsequent leapfrog terms can be contracted into one such that the total number of required operator exponentials is indeed $m$.

Depending on the number $m$ of factors in a decomposition $U$ of type (37)-(43), the parameters can be chosen such that $U$ coincides with the exact $e^{tH}$ up to order $p$ in the sense of Eq. (6). If free parameters remain, we can use these to minimize the leading error term of the decomposition.

Of course, for the same number of factors $m$, some types of decompositions are subclasses of others:

[TABLE]

V.2 Numbers of parameters and constraints, symmetries

As in Eqs. (10), we can use the BCH formula recursively, to compute $\log U$ in terms of nested commutators of $A$, $B$, and $C$. The number of constraints, imposed by requiring $U$ to coincide with $e^{tH}$ up to order $t^{p}$ [Eq. (9)], can be determined from the number of terms in the expansion. For the large Hilbert spaces relevant in many-body physics and low expansion orders considered for our optimized decompositions, it is sufficient to treat the Lie algebra generated by $\{A,B,C\}$ as free. We can hence work with Hall bases Hall1950-1 ; Serre1992 ; Reutenauer1993 as discussed in Sec. III.3. The number of Hall basis elements of degree $k$ is given by the necklace polynomial (12). Tables 1 and 2 list Hall basis elements and their numbers $d_{k}$.

Assuming that all terms that can occur do occur in the expansion of $Z=\log U$ for an order-$p$ Lie-Trotter-Suzuki decomposition $U$, we read off constraint polynomials as coefficients of the relevant Hall basis elements. To check and determine how many free parameters we actually have, we can employ the Gröbner basis of the constraint polynomials as discussed in Sec. III.5.

Type N. – For the decompositions $U_{\text{N},m}$ in Eq. (37), we have $\nu=m$ parameters $\{a_{i},b_{i},c_{i}\}$. As the decomposition has no further symmetry or structure, all Hall basis elements should occur in the expansion of $Z$ and the number of constraints to achieve approximation order $p$ is $\sum_{k=1}^{p}d_{k}$. For example, at first order, we have the three constraints $\sum_{i}a_{i}=\sum_{i}b_{i}=\sum_{i}c_{i}=1$ to achieve $Z=tH+\mathcal{O}(t^{2})$.

Type S. – For the decompositions $U_{\text{S},m}$ in Eq. (38), we have $\nu=(m+1)/2$ parameters and even $(m-1)/2$. Due to the symmetry in the factors, the decompositions obey the time reversal symmetry $U_{\text{S},m}(t)U_{\text{S},m}(-t)=\mathbbm{1}$. According to Lemma 3, it follows that $Z(t)$ only contains terms of odd order in $t$ such that no Hall basis elements of even degree occur Yoshida1990 . For a symmetric decomposition that has no further structure, all Hall basis elements of odd degree should occur in the expansion of $Z$ and the number of constraints to achieve approximation order $p$ is $\sum_{q,2q-1\leq p}d_{2q-1}$. For example, at first order, we have the three constraints that the coefficients $a_{i}$ in the $A$ factors, $b_{i}$ in the $B$ factors and, $c_{i}$ in the $C$ factor of Eq. (38) sum to one to achieve $Z=tH+\mathcal{O}(t^{2})$.

Type S-abc. – For the decompositions $U_{\text{S}^{\prime},m}$ in Eq. (39), we have $\nu=\operatorname{ceil}(m/2)$ parameters and even $(m+1)/3$. As these decompositions are symmetric, only terms of odd degree occur in the expansion of $Z=\log U_{\text{S}^{\prime},m}$. As for type S, the number of constraints to achieve approximation order $p$ is $\sum_{q,2q-1\leq p}d_{2q-1}$.

Type SE. – For the decompositions $U_{\text{SE},m}$ in Eq. (41), we have $\nu=(m-1)/4$ parameters and even $(m-1)/2$. As these decompositions are symmetric, only terms of odd degree occur in the expansion of $Z=\log U_{\text{SE},m}$. The decomposition is a product of Euler terms $U_{\text{E}+}(\tau)$ and $U_{\text{E}-}(\tau)$ in Eq. (40) which are non-symmetric first order decompositions of $e^{\tau H}$ Trotter1959 . To determine the number of constraints, we can use the following expansions.

Lemma 4** (Euler term expansions).**

The expansions of the forward and backward Euler terms $U_{{\text{E}}+}(\tau)=e^{\tau A_{1}}\,e^{\tau A_{2}}\dotsb e^{\tau A_{n}}$ and $U_{{\text{E}}-}(\tau)=e^{\tau A_{n}}\,e^{\tau A_{n-1}}\dotsb e^{\tau A_{1}}$ coincide up to sign factors. In particular,

[TABLE]

For the proof, note that $U_{{\text{E}}+}(\tau)U_{{\text{E}}-}(-\tau)=\mathbbm{1}$ and let us define $\log U_{{\text{E}}\pm}(\tau)=:\tau Z^{(1)}_{\text{E}}+\tau^{2}Z^{(2)}_{\pm}+\tau^{3}Z^{(3)}_{\pm}+\tau^{4}Z^{(4)}_{\pm}+\dots$ where $Z^{(1)}_{\text{E}}=\sum_{i=1}^{n}A_{i}$. Applying the BCH formula (8), we find $\log[U_{{\text{E}}+}(\tau)U_{{\text{E}}-}(-\tau)]=\tau^{2}(Z^{(2)}_{+}+Z^{(2)}_{-})+\mathcal{O}(\tau^{3})$ such that $Z^{(2)}_{+}=-Z^{(2)}_{-}=:Z^{(2)}_{{\text{E}}+}$ because of $U_{{\text{E}}+}(\tau)U_{{\text{E}}-}(-\tau)=\mathbbm{1}$. With this information, we can reconsider the BCH formula and find $\log[U_{{\text{E}}+}(\tau)U_{{\text{E}}-}(-\tau)]=\tau^{3}(Z^{(3)}_{+}-Z^{(3)}_{-})+\mathcal{O}(\tau^{4})$, showing that $Z^{(3)}_{+}=Z^{(3)}_{-}=:Z^{(3)}_{{\text{E}}+}$. Continuing in this way, Eq. (45) is established.

$Z^{(k)}_{\text{E}}$ is a Lie polynomial containing only nested commutators of degree $k$. To determine the number of constraints for $U_{\text{SE},m}$ to be an order-$p$ decomposition of $e^{tH}$, we can consider the Lie algebra generated by $\{Z^{(1)}_{\text{E}},Z^{(2)}_{\text{E}},Z^{(3)}_{\text{E}},\dots\}$ as free. The number of constraints due to Eq. (9) is then given by the number $\sum_{q,2q-1\leq p}d_{{\text{E}},2q-1}$ of Hall basis elements of odd degree $2q-1\leq p$. Tables 6 and 7 list the Hall basis elements and their numbers $d_{{\text{E}},k}$. At first order, we only have the constraint that the coefficients $u_{i}$ and $v_{i}$ of all Euler factors in Eq. (41) sum to one to achieve $Z=tH+\mathcal{O}(t^{2})$.

Type SL. – For the decompositions $U_{\text{SL},m}$ in Eq. (43), we have $\nu=\operatorname{ceil}((m-1)/8)$ parameters and integer $(m-1)/4$. As these decompositions are symmetric, only terms of odd degree occur in the expansion of $Z=\log U_{\text{SL},m}$. The decomposition is a product of leapfrog terms $U_{\text{L}}(\tau)$ in Eq. (42) which are symmetric second order decompositions of $e^{\tau H}$ and hence can be expanded in the form $\log U_{\text{L}}=\tau Z^{(1)}_{\text{L}}+\tau^{3}Z^{(3)}_{\text{L}}+\tau^{5}Z^{(5)}_{\text{L}}+\dots$ with $Z^{(1)}_{\text{L}}=H$ [Eq. (13)]. $Z^{(k)}_{\text{L}}$ is a Lie polynomial containing only nested commutators of degree $k$. To determine the number of constraints for $U_{\text{SL},m}$ to be an order-$p$ decomposition of $e^{tH}$, we can consider the Lie algebra generated by $\{Z^{(1)}_{\text{L}},Z^{(3)}_{\text{L}},Z^{(5)}_{\text{L}},\dots\}$ as free. The number of constraints due to Eq. (9) is then given by the number $\sum_{q,2q-1\leq p}d_{\text{L},2q-1}$ of Hall basis elements of odd degree $2q-1\leq p$. Tables 3 and 4 list the Hall basis elements and their numbers $d_{\text{L},k}$. At first order, we only have the constraint that the coefficients $w_{i}$ of all leapfrog factors in Eq. (43) sum to one to achieve $Z=tH+\mathcal{O}(t^{2})$.

To summarize this section, we give the total numbers of constraints by order in Table 8. The table also states $m_{p}^{\text{H}}$, the minimum number of factors needed in the different decompositions to obtain an order-$p$ decomposition [cf. Eq. (6)], under the assumption that all constraint polynomials are independent.

VI Optimized decompositions for $n=3$ terms

As discussed in Sec. III.6 we can quantify the accuracy of an order-$p$ decomposition by the leading error term which is of order $t^{p+1}$. The relevant error measure is the operator norm distance which we bound using the triangle inequality. This leads to Eq. (15) and the 1-norm $\sum_{i}|c_{i}|$ of the Hall basis expansion coefficients $|c_{i}|$ as an error measure. Furthermore, we are free to impose any order for the generators $A,B,C$ in the construction of the Hall basis and the $c_{i}$ depend on that choice. Let the sets of coefficient polynomials for the six possible orders be denoted by $\{c^{ABC}_{i}\}$, $\{c^{BCA}_{i}\}$ etc. We will then use the error measure

[TABLE]

to quantify the magnitude of the deviation of $U$ from $e^{tH}$, where $m$ denotes again the number of factors in the decomposition $U$. Like the coefficient polynomials $c_{i}$, the error $\varepsilon$ is a function of the parameters in the decomposition.

Here, we optimize the decompositions of types N, S, S-abc, SE, and SL defined in Sec. V.1 with respect to their parameters to minimize the error measure (46). We compare the results to the unoptimized type-SL decompositions (18) and (19) of Yoshida Yoshida1990 and Suzuki Suzuki1990 which generalize without modification to arbitrary numbers $n$ of terms in $H$. Similarly, we compare to type-SL decompositions of McLachlan McLachlan1995 , adapted to the $n=3$ case, and those of Kahan and Li Kahan1997-66 .

In the discussion below, $\nu$ denotes the number of parameters as stated in Sec. V.1 and the number of constraints for the different decomposition types and approximation orders $p$ are given in Table 8. We check the applicability of the corresponding counting argument using Gröbner bases for the constraint polynomials as discussed in Sec. III.5. For each order $p$, the best decomposition found (smallest $\varepsilon$) is indicated by a star. When the different possible orders $A<B<C$, $B<C<A$ etc., chosen for the construction of the Hall basis, are not equivalent, we specify in brackets the one that yields the minimum for $\varepsilon$.

VI.1 Order $p=1$

$\bigstar$

Euler ($m=3$, type N, $\nu=1$). – All three parameters are fixed to $a_{1}=b_{1}=c_{1}=1$ by the constraint that the first order term in $\log U$ is $tH=t(A+B+C)$. The error is $\varepsilon=\left(\frac{3}{1}\right)^{1}\frac{3}{2}=\frac{9}{2}$.

VI.2 Order $p=2$

• Leapfrog ($m=5$, type S, $\nu=1$). – The only parameter is fixed to $w_{1}=1$ by the constraint that the first order term in $\log U$ is $tH=t(A+B+C)$. The error is $\varepsilon=\left(\frac{5}{2}\right)^{2}\frac{13}{24}=\frac{325}{96}\approx 3.3854$ ($A<B$ and $A<C$).

$\bigstar$

Optimized ($m=9$, type S, $\nu=5$). – There are three constraints and hence two free parameters; we choose $\{a_{1},b_{1}\}$. The error is minimized for

[TABLE]

which gives $\varepsilon\approx 1.0496$ ($B<A$ or $C<A$). At this point, the coefficients of the terms $[B,[A,B]]$, $[B,[B,C]]$, and $[C,[A,C]]$ vanish. See Fig. 6(a).

• Optimized ($m=11$, type S-abc, $\nu=6$). – There are three constraints and hence three free parameters; we choose $\{a_{1},b_{1},c_{1}\}$. The error is minimized for

[TABLE]

which gives $\varepsilon\approx 2.3391$. At this point, the coefficients of the terms $[B,[A,C]]$ and $[C,[A,B]]$ vanish.

• Optimized ($m=11$, type S, $\nu=6$). – There are three constraints and hence three free parameters; we choose $\{a_{1},b_{1},b_{2}\}$. The error is minimized for

[TABLE]

which gives $\varepsilon\approx 1.3054$ ($B<A$ or $C<A$). At this point, the coefficients of the terms $[A,[A,B]]$, $[B,[B,C]]$, and $[C,[A,C]]$ vanish.

$\blacksquare$

Discussion. – The best decomposition found here is of type S with $m=9$ as specified in Eq. (47). It improves over the common leapfrog decomposition by a factor $\sim 1/3$. According to Table 8, we can reach order $p=4$ with $m=13$ factors.

VI.3 Order $p=4$

• Forest & Ruth, Yoshida ($m=13$, type SL). – For the decomposition $U_{\text{Y},q=2}$ with $q=2$ in Eq. (18), the error is $\varepsilon\approx 65.721$ ($A<C<B$).

• Optimized ($m=17$, type S-abc, $\nu=9$). – There are no free parameters and two real solutions. The first solution reproduces the type-SL decomposition above with $m=13$. The second solution has a very large error $\varepsilon\approx 1243.13$.

• Optimized ($m=17$, type SL, $\nu=2$). – There are two constraints and hence no free parameters. According to the Gröbner basis, there are only two complex solutions.

• Optimized ($m=17$, type SE, $\nu=4$). – There are three constraints and hence there is one free parameter. For practical reasons, we reparametrize according to

[TABLE]

• Suzuki, McLachlan, Kahan & Li, Omelyan et al. ($m=21$, type SL). – For Suzuki’s decomposition $U_{\text{Z},q=2}$ in Eq. (19), the error is $\varepsilon\approx 35.239$. For McLachlan’s decomposition (26) and the decomposition (27) of Omelyan et al., adapted to the $n=3$ case, we find errors $\varepsilon\approx 19.479$ and $\varepsilon\approx 22.827$, respectively. Kahan and Li’s decompositions (25) have the error $\varepsilon\approx 33.346$ (order $A<B<C$ in all cases).

• Optimized ($m=21$, type SL, $\nu=3$). – There are two constraints and hence there is one free parameter; according to the Gröbner basis, we can choose $w_{2}$. The error is minimized for

[TABLE]

$\bigstar$

Optimized ($m=21$, type SE, $\nu=5$). – There are three constraints and hence two free parameters. For practical reasons, we reparametrize according to

[TABLE]

• Optimized ($m=25$, type SL, $\nu=3$). – There are two constraints and hence there is one free parameter; according to the Gröbner basis, we can choose $w_{2}$. The error is minimized for

[TABLE]

which gives $\varepsilon\approx 56.179$ ($A<C<B$).

• Optimized ($m=25$, type SE, $\nu=6$). – There are three constraints and hence there are three free parameters. For practical reasons, we reparametrize according to

[TABLE]

$\blacksquare$

Discussion. – The best decomposition found here is of type SE with $m=25$ as specified in Eq. (54). It improves over the decomposition $U_{\text{Y},q=2}$ [Eq. (18)] of Forest, Ruth, and Yoshida Forest1990-43 ; Yoshida1990 by a factor $\sim 1/20$ and by a factor of $\sim 1/10$ over the decomposition $U_{\text{Z},q=2}$ [Eq. (19)] due to Suzuki Suzuki1990 . Actually, the error of the decomposition (54) with $m=25$ factors does not improve too much over that of the type-SE decomposition with $m=21$ factors [Eq. (52)]. We hence recommend using the latter for fourth order integration. It is not surprising that the optimized type-SL decomposition (53) with $m=25$ factors has a larger error than the SL decomposition (51a) with $m=21$ factors. The number of free parameters does not increase when going to $m=25$, but the increased number of factors is taken account of in the definition (16) of $\varepsilon$ and results in a larger error value. We have checked explicitly that type-S decompositions with $m=13$ factors just reproduce the type-SL decomposition $U_{\text{Y},q=2}$ due to Forest, Ruth, and Yoshida, and that type-S decompositions with $m=17$ and $m=21$ reproduce the type-SE decompositions (50) and (52), respectively. Similarly, all inspected solutions for type-S decompositions with $m=25$ were of type SE. According to Table 8, we can reach order $p=4$ with $m=29$ factors.

VI.4 Order $p=6$

• Optimized ($m=29$, type SL, $\nu=4$). – There are four constraints and hence no free parameters. According to the Gröbner basis, there are three real solutions. The best of the three solutions is

[TABLE]

which has error $\varepsilon\approx 722.85$ ($A<B<C$). It, of course, coincides with the corresponding decomposition (33) for $n=2$ terms. The other two solutions have much larger errors $\varepsilon\approx 15940$ and $\varepsilon\approx 16470$, respectively (both for order $A<B<C$).

• Yoshida, Kahan & Li ($m=37$, type SL). – For Yoshida’s decomposition $U_{\text{Y},q=3}$ with $q=3$ in Eq. (18), the error is $\varepsilon\approx 68024$ ($A<B<C$). Kahan and Li’s decomposition (34) has the error $\varepsilon\approx 687.06$ ($A<B<C$).

$\bigstar$

Optimized ($m=37$, type SL, $\nu=5$). – There are four constraints and hence there is one free parameter; according to the Gröbner basis, we can choose $w_{4}$. The error is minimized for

[TABLE]

• Suzuki ($m=101$, type SL). – For the decomposition $U_{\text{Z},q=3}$ in Eq. (19), the error is $\varepsilon\approx 51034$ ($A<B<C$).

$\blacksquare$

Discussion. – The best decomposition found here is of type SL with $m=37$ as specified in Eq. (56a). It improves over Yoshida’s generic decomposition $U_{\text{Y},q=3}$ [Eq. (18)] by a factor $\sim 1/165$, by a factor of $\sim 1/124$ over Suzuki’s decomposition $U_{\text{Z},q=3}$ [Eq. (19)], and by a factor of $\sim 4/7$ over the simplest sixth order decomposition (55). According to Table 8, true type-SE and type-S decompositions require $m=37$ and $m=117$ factors, respectively, to reach order $p=6$, and there exist type-SL decompositions of order $p=8$ for $m\geq 61$ factors.

VII Discussion

We have determined optimized Lie-Trotter-Suzuki decompositions for $n=2$ and $n=3$ terms up to order $t^{6}$. Using a coarse-graining argument, we have explained why these decompositions are sufficient to simulate any 1d and 2d lattice models with finite-range interactions. Decompositions of different approximation order are constructed by expanding in terms of nested commutators, using Hall bases to remove linear dependencies, and solving systems of polynomial constraints resulting from the comparison with $e^{tH}$. The sizes of Hall bases are also essential to understand the numbers of constraints and free parameters for the different decomposition types. The free parameters are used to minimize the amplitudes of leading error terms. For these optimizations, we employ an error measure that bounds the operator-norm distance and allows for a fair comparison of decompositions with different numbers of factors $m$ in the sense that the time step $t$ should be chosen proportional to $m$ to keep computation costs constant.

For $n=2$ terms, at order $p=2$, we recommend the type-S decomposition (21) with $m=5$ factors, at order $p=4$, the type-S decomposition (30a) with $m=11$ factors and, at order $p=6$, the type-SL decomposition (35) with $m=19$ factors. We have applied the recommended $p=4$ decomposition, in particular, in many precise tensor network simulations as in Refs. Barthel2013-15 ; Lake2013-111 ; Cai2013-111 ; Barthel2016-94 ; Barthel2017_08unused ; Binder2018-98 .

For $n=3$ terms, at order $p=2$, we recommend the type-S decomposition (47) with $m=9$ factors, at order $p=4$, the type-SE decomposition (52) with $m=21$ factors and, at order $p=6$, the type-SL decomposition (56a) with $m=37$ factors.

Ref. Sornborger1999-60 discusses decompositions for an arbitrary number of terms $n$. For $n=2$ and $n=3$ they have errors similar to those of the decomposition $U_{\text{Y},q=2}$ in Eq. (18) due to Forest, Ruth, and Yoshida Forest1990-43 ; Yoshida1990 . Note that the type-SL and type-SE decompositions presented here are generally applicable for any number of terms $n$, but they are in general not optimal when applied for $n\geq 4$.

Of course there are alternatives to using Lie-Trotter-Suzuki decompositions. Tensor network states and matrix product states, in particular, can also be evolved using Runge-Kutta methods Cazalilla2001 ; Feiguin2005 , Krylov subspace methods Schmitteckert2004-70 ; Garcia-Ripoll2006-8 ; Dargel2012-85 ; Wall2012-14 , or the time-dependent variational principle Haegeman2011-107 ; Haegeman2016-96 . Some reviews are given in Refs. Garcia-Ripoll2006-8 ; Schollwoeck2011-326 ; Paeckel2019_01 . For the purpose of digital quantum simulation (a.k.a. Hamiltonian simulation), algorithms with a gate count that is poly-logarithmic in the desired accuracy have been developed Berry2014-283 ; Berry2015-114 ; Low2017-118 ; Haah2018_01 , e.g., by introducing ancillary qubits and implementing truncated Taylor expansions. For classical systems, popular choices are linear multistep methods and Runge-Kutta methods.

We gratefully acknowledge discussions with R. Mosseri and J. Socolar, and support through US Department of Energy grant DE-SC0019449.