Greatly improved higher-order product formulae for quantum simulation

Mauro E. S. Morales Centre for Quantum Computation and Communication Technology, University of Technology Sydney, Sydney, NSW 2007, AU Centre for Quantum Software and Information, University of Technology Sydney, Sydney, NSW 2007, AU Joint Center for Quantum Information and Computer Science (QuICS), University of Maryland, College Park, Maryland 20742, USA Pedro C. S. Costa School of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109, AU Quantum for New South Wales, Sydney, NSW 2000, AU Giacomo Pantaleoni School of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109, AU Daniel K. Burgarth School of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109, AU Department Physik, Friedrich-Alexander-Universität Erlangen-Nürnberg, Staudtstraße 7, 91058 Erlangen, Germany Yuval R. Sanders Centre for Quantum Software and Information, University of Technology Sydney, Sydney, NSW 2007, AU School of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109, AU Dominic W. Berry School of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109, AU

Abstract

Quantum algorithms for simulation of Hamiltonian evolution are often based on product formulae. The fractal method of Suzuki gives a systematic way to find arbitrarily high-order product formulae, but results in a large number of exponentials. On the other hand, product formulae with fewer exponentials can be found by numerical solution of simultaneous nonlinear equations. It is also possible to reduce the cost of long-time simulations by processing, where a kernel is repeated and a processor need only be applied at the beginning and end of the simulation. In this work, we found thousands of new product formulae of both 8th and 10th order, and numerically tested these formulae, together with many formulae from prior literature. We provide methods to fairly compare product formulae of different lengths and different orders. We have found a new 8th order processed product formula with exceptional performance, that outperforms all other tested product formulae for about eight orders of magnitude in system parameters $T$ (time) and $\epsilon$ (allowable error). That includes most reasonable combinations of parameters to be used in quantum algorithms.

I Introduction

The Lie-Trotter product formula is commonly used in quantum algorithms for Hamiltonian simulation, where one can approximate the Hamiltonian evolution of $H=A+B$ in terms of exponentials of $A$ and $B$ when these operators do not commute. For short time, a standard approximation is the second-order symmetric formula $S_{2}(t)=e^{-iAt/2}e^{-iBt}e^{-iAt/2}$ , which satisfies $e^{-iHt}=S_{2}(t)+\mathcal{O}(t^{3})$ . More generally, the error in an order $k$ formula is $\mathcal{O}(t^{k+1})$ . Longer times are simulated by breaking the evolution into many repetitions of a short time. The number of repetitions needed is reduced with the order, motivating the search for higher-order product formulae. A systematic method to produce arbitrarily high order formulae is the fractal method of Suzuki [1, 2], which has found applications in Hamiltonian simulation [3]. The first explicit use of product formulae for quantum simulation was given in [4], applying it for quantum evolution under local Hamiltonians. Later work considered the broader class of sparse Hamiltonians [5] and higher orders [3], as well as methods beyond product formulae [6, 7, 8, 9].

Recent work has shown that despite its simplicity, the Lie-Trotter product formula can compete with other Hamiltonian simulation algorithms due to the low error that it achieves in practice [10]. Methods based on linear combinations of unitaries [11, 12] or quantum signal processing [13] give complexity logarithmic in the inverse error, but the error is not required to be extremely small, meaning those methods do not provide as large an advantage as might be expected. Product formula error bounds can be further improved by considering particular physical systems [14, 15, 16] or leveraging randomisation [17, 18, 19]. Moreover, Trotter formulae are expected to be relevant for both noisy intermediate-scale quantum computation and fault-tolerant computation. It is then of great importance to seek efficient implementations of product formulae as it can have a great impact on the efficiency of Hamiltonian simulation algorithms in practice.

The downside of the Suzuki method to generate higher-order formulae is that the number of exponential operators to implement it grows very rapidly. Suzuki’s product formulae are usually assumed in quantum computing, but they can be greatly improved upon. An alternative method by Yoshida [20] can be used to obtain product formulae with a smaller number of exponentials. Similar to Suzuki’s formulae, they are given as a product of $S_{2}$ for different time intervals, but in contrast to Suzuki’s approach there is not an explicit analytic form for the higher-order formulae. Instead one needs to derive and solve a complicated set of simultaneous nonlinear polynomial equations.

Yoshida gives what appears to be all 6th order solutions but only some 8th order integrators, and did not consider any orders beyond that. More recent work [21, 22, 23, 24] has pushed the search to higher orders. In Ref. [23] a summary of product formulae in the literature is given, including results from orders 4 to 10. The terminology used in that work is stages for the number of $S_{2}$ in the product formula. In Ref. [21] the authors provide 8th order product formulae with 15 and 17 stages which improve over those of Suzuki and Yoshida. Other product formulae based on the ansatz of Yoshida are also presented in Ref. [22]. For 8th order, Ref. [22] replicates the findings of [21] for product formulae with 15 and 17 stages, and finds new solutions with 19 and 21 stages. For 10th order, Ref. [22] provides solutions with 31 and 33 stages.

Another approach is that of processed product formulae [25, 26, 27, 28, 29, 30, 24, 31]. Instead of the product formula being symmetric, it is of the form $P\Sigma P^{-1}$ for kernel $\Sigma$ and processor $P$ . The $P$ and $P^{-1}$ cancel when using the product of many of these for evolution over a long time, so the cost is dominated by the number of stages in the kernel. Since the kernel has fewer restrictions on it, it can have fewer stages than a normal symmetric product formula.

The objective of this work is to find improved product formulae and compare their performance in a consistent manner. We show that the performance of product formulae is better quantified by the error in the eigenvalues, rather than the spectral-norm error as usually considered in prior work. This is because it is the eigenvalue error that dominates the error for evolution at longer times. We also derive a method to consistently compare the performance of product formulae with different numbers of exponentials.

We use these methods to compare between the performance of our product formulae, as well as to compare to prior product formulae in the literature. By our numerical search, we found hundreds of thousands of solutions at 8th order, including some that outperform those previously reported in the literature, with our best being a processed product formula. For 10th order we have found thousands of new solutions including some with extremely low error, though one given in Ref. [22] provides even better performance. A detailed list of the product formulae considered in this work can be found in Section VI.

When comparing product formulae of different orders, it is better to use higher-order formulae for larger values of $T/\epsilon$ , where $T$ is the total evolution time and $\epsilon$ is the required precision. For smaller $T/\epsilon$ it is better to use lower-order formulae, but as it is increased there are threshold values beyond which it is optimal to use higher-order formulae. We derive a methods for determining these thresholds, and show that our 8th order product formula is best for $T/\epsilon$ from about $10^{6}$ to $10^{14}$ , which includes the range of typical values for quantum chemistry applications. This means that the best 8th order product formula we have found in this work will be best suited to real applications. It is possible to reduce the lower threshold by an order of magnitude by adjusting the product formula to provide low error for larger time steps.

The organisation of this work is as follows. In Section II we introduce the background necessary for the rest of the paper. First, we define product formulae and introduce the Suzuki method for generating higher-order product formulae in Section II.1, then in Section II.2 we give a summary of Yoshida’s method, then we explain processors for product formulae in Section II.3. In Section III we present the methods used for solution, including Subsection III.2 where we introduce the Taylor series method to find new product formulae, which is distinct from Yoshida’s method. In Sections IV.1 and IV.2 we present the results regarding new 8th and 10th order product formulae, and then in Section VI we give a comprehensive comparison of the product formulae.

II Background

In this section, we give a summary of the background for our work. We begin by defining product formulae and the Baker-Campbell-Hausdorff formula, then we introduce the Suzuki method and Yoshida method to obtain higher-order formulae, and describe the processed product formulae.

II.1 Product formulae

It is well known that, for any non-commuting operators $X$ and $Y$ ,

\exp((X+Y)t)=\exp(Xt)\exp(Yt)+\order{t^{2}}.

(1)

where we have absorbed the $-i$ factor used in quantum simulation into $X$ and $Y$ . The above equation demonstrates that the exponential of a sum of two operators is, to first order, equal to the product of the exponential of those operators. The above equation is often referred as a ‘first-order product formula’. Higher-order terms can be computed via the Baker-Campbell-Haussdorff (BCH) formula.

Theorem 1 (Baker-Campbell-Haussdorff formula [32]).

Let $X$ and $Y$ be any operators such that $\norm{X}+\norm{Y}<\ln{2}$ . We have for an operator $Z$ that $\exp(X)\exp(Y)=\exp(Z)$ , with

Z=\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}\sum_{\begin{subarray}{c}r_{1}+s_{1}>0\\ \vdots\\ r_{n}+s_{n}>0\end{subarray}}\frac{[X^{r_{1}},Y^{s_{1}},\cdots X^{r_{n}},Y^{s_{n}}]}{\left(\sum_{j=1}^{n}r_{i}+s_{i}\right)\prod_{i=1}^{n}r_{i}!s_{i}!},

(2)

where

[X^{r_{1}},Y^{s_{1}},\dotsm X^{r_{n}},Y^{s_{n}}]=[\underbrace{X,[X,\dotsm[X}_{r_{1}},[\underbrace{Y,[Y,\dotsm[Y}_{s_{1}},\,\dotsm\,[\underbrace{X,[X,\dotsm[X}_{r_{n}},[\underbrace{Y,[Y,\dotsm Y}_{s_{n}}]]\dotsm]].

The standard second-order symmetric product formula is as given in the definition below.

Definition 2 (Symmetric product formula of order two).

Let $X$ and $Y$ be non-commuting operators and let $t$ be a real variable. Define

S_{2}(t):=\exp(\frac{1}{2}Xt)\exp(Yt)\exp(\frac{1}{2}Xt).

(3)

The operators $X$ and $Y$ used in the definition of $S_{2}$ should always be clear from context. More generally, when there is a sum of $X_{j}$ , the product formula is

S_{2}(t):=\left[\prod_{j=1}^{J}\exp(\frac{1}{2}X_{j}t)\right]\left[\prod_{j=J}^{1}\exp(\frac{1}{2}X_{j}t)\right].

(4)

Products are ordered with the starting index on the right and the final one on the left, so for $J=2$ terms the expression in the definition is obtained. Note that $S_{2}$ , and all product formulae considered in this work, have an order independent of $J$ , so there may be any number of terms in the sum. The corollary below describes the form of the error terms in the symmetric product formula, and also implies that it is second order.

Corollary 3 (Symmetric BCH formula [20]).

Let $X$ and $Y$ be any operators such that $\norm{X}+\norm{Y}<\ln{2}$ and let $t$ be a real variable. Define $Z$ such that $S_{2}(t)=\exp(Z)$ . Then there is a sequence $\alpha_{\ell}$ consisting of linear combinations of $\ell$ -term commutators in $X$ and $Y$ such that

Z=\sum_{\ell=1}^{\infty}\alpha_{\ell}t^{\ell}.

(5)

Moreover, $\alpha_{\ell}\equiv 0$ whenever $\ell$ is even.

Reference [20] also shows that even terms are zero for more general symmetric product formulae. The first three non-zero $\alpha_{\ell}$ terms from above are

$\displaystyle\alpha_{1}$	$\displaystyle=X+Y,$	(6)
$\displaystyle\alpha_{3}$	$\displaystyle=\frac{1}{12}[Y,[Y,X]]-\frac{1}{24}[X,[X,Y]],$	(7)
$\displaystyle\alpha_{5}$	$\displaystyle=\frac{7}{5760}[X,X,X,X,Y]-\frac{1}{720}[Y,Y,Y,Y,X]+\frac{1}{360}[X,Y,Y,Y,X]+\frac{1}{360}[Y,X,X,X,Y]-\frac{1}{480}[X,X,Y,Y,X]$
	$\displaystyle\quad+\frac{1}{120}[Y,Y,X,X,Y].$	(8)

Here the square brackets are used to indicate multicommutator expressions similar to the notation in Theorem 1, for example

[Y,Y,X,X,Y]\equiv[Y,[Y,[X,[X,Y]]]].

(9)

Expressions for each $\alpha_{\ell}$ can be derived from two applications of the usual BCH formula.

Definition 4 (Product formula).

Let $X$ and $Y$ be any non-commuting operators. Given a natural number $n$ , a product formula of order $n$ is a pair $(\vec{c},\vec{d})$ with $\vec{c},\vec{d}\in\mathbb{R}^{\ell}$ and $\ell$ a natural number such that for all $t\in\mathbb{R}$

\exp((X+Y)t)=\prod_{j=1}^{\ell}\exp(c_{j}Xt)\exp(d_{j}Yt)+\order{t^{n+1}}.

(10)

We refer to the number of non-zero coefficients in $\left(\vec{c},\vec{d}\right)$ as the length of the product formula.

Hence $S_{2}$ is a length- $3$ product formula. The number of exponentials used in the product formula is a crucial measure of its efficiency, and in quantum simulation it is proportional to the required number of gates.

Suzuki method for generating higher-order product formulae.

Here we describe Suzuki’s fractal methods from [1, 2] to obtain higher-order product formulae. Starting from the symmetrised product formula in Eq. (3), the fractal method produces product formulae at all even orders. Suzuki’s first fractal method to generate a product formula of order $k=2\kappa$ is [1]

S_{2\kappa}(t)=S_{2\kappa-2}(s_{\kappa}t)S_{2\kappa-2}((1-2s_{\kappa})t)S_{2\kappa-2}(s_{\kappa}t),

(11)

where $s_{\kappa}=1/(2-2^{1/(2\kappa-1)})$ . This method can be used to generate even orders starting at $S_{2}$ . A drawback to this method is that both $s_{k}$ and $1-2s_{\kappa}$ are greater than 1, so the coefficients in the higher-order methods are large, causing greater error.

Alternatively to Eq. 11, an order $2\kappa$ product formula can be obtained via Suzuki’s second fractal method [2]

\tilde{S}_{2\kappa}(t)=\tilde{S}_{2\kappa-2}(u_{\kappa}t)^{2}\tilde{S}_{2\kappa-2}((1-4u_{\kappa})t)\tilde{S}_{2\kappa-2}(u_{\kappa}t)^{2},

(12)

where $u_{\kappa}=1/(4-4^{1/(2\kappa-1)})$ . This method has the advantage that both $u_{\kappa}$ and $1-4u_{\kappa}$ are less than 1, so the coefficients of higher-order formulae are small resulting in small error. The drawback is that far more exponentials are required. Each iteration uses 5 copies of the lower-order formula, whereas the previous one uses 3 copies. The virtue of these fractal methods is that they allow one to generate arbitrarily high-order product formulae easily, albeit at the expense of a large number of exponentials.

Exponential length scaling of the Suzuki method.

For Suzuki’s first method, the total number of exponentials for a given order $2\kappa={4,6,8,...}$ in the product formula $S_{2\kappa}$ is given by

2(J-1)3^{\kappa-1}+1\,.

(13)

For example $J=2$ for $X+Y$ and $\kappa=1$ just corresponds to $S_{2}$ , and the expression gives 3 as expected. For Suzuki’s second method $\tilde{S}_{2\kappa}$ the number of exponentials is

2(J-1)5^{\kappa-1}+1.

(14)

The number of exponential operators in both cases of the Suzuki method grows very rapidly. Thus one may be interested in alternative method to obtain product formulae with a lower count, such as the method of Yoshida in the next section.

II.2 Yoshida’s method for deriving 6th order product formulae

Approach.

Rather than using Eqs. 11 and 12, Yoshida uses the general procedure

\displaystyle S^{(m)}(t)=\bigg{(}\prod_{j=1}^{m}S_{2}(w_{m-j+1}t)\bigg{)}S_{2}(w_{0}t)\bigg{(}\prod_{j=1}^{m}S_{2}(w_{j}t)\bigg{)},

(15)

where $w_{j}\in\mathbb{R}$ for $j=0,1,\cdots,m$ are parameters to be determined. Note the number of exponentials in this product is given by $(4m+2)(J-1)+1$ . Given this ansatz, the task becomes to find $m$ and $w_{i}$ such that $S^{(m)}$ is an order $k$ product formula. We will illustrate Yoshida’s method by deriving the result for 6th order.

Expand Yoshida product using Baker-Campbell-Haussdorf formula.

The method is to expand Eq. 15 using the BCH formula from Theorem 1 recursively as follows. First, note that from Corollary 3, $S_{2}(t)=e^{\frac{t}{2}X}e^{tY}e^{\frac{t}{2}X}=e^{t\alpha_{1}+t^{3}\alpha_{3}\cdots}$ , where $\alpha_{\ell}$ is a commutator of $\ell$ operators as described below Corollary 3. We are for the moment considering sums of two terms $X+Y$ . Define $C=\sum_{i=1}^{\infty}t^{2i-1}w_{1}^{2i-1}\alpha_{2i-1}$ and $D=\sum_{i=1}^{\infty}t^{2i-1}w_{0}^{2i-1}\alpha_{2i-1}$ . Then,

		$\displaystyle S_{2}(w_{1}t)S_{2}(w_{0}t)S_{2}(w_{1}t)$
		$\displaystyle=e^{C}e^{D}e^{C}$
		$\displaystyle=\exp\{tw_{1}\alpha_{1}+t^{3}w_{1}^{3}\alpha_{3}+t^{5}w_{1}^{5}\alpha_{5}+\order{t^{7}}\bigg{\}}\exp\{tw_{0}\alpha_{1}+t^{3}w_{0}^{3}\alpha_{3}+tw_{0}^{5}\alpha_{5}+\order{t^{7}}\bigg{\}}$
		$\displaystyle\quad\times\exp\{tw_{1}\alpha_{1}+t^{3}w_{1}^{3}\alpha_{3}+t^{5}w_{1}^{5}\alpha_{5}+\order{t^{7}}\bigg{\}}$
		$\displaystyle=\exp\{t(2w_{1}+w_{0})\alpha_{1}+t^{3}(2w_{1}^{3}+w_{0}^{3})\alpha_{3}+t^{5}(2w_{1}^{5}+w_{0}^{5})\alpha_{5}+\frac{1}{6}([D,D,C]-[C,C,D])+\order{t^{7}}\bigg{\}}.$		(16)

A simple computation shows

\displaystyle[D,D,C]-[C,C,D]

\displaystyle=t^{5}(w_{0}^{2}w_{1}^{3}-w_{1}^{2}w_{0}^{3}+w_{1}^{4}w_{0}-w_{0}^{4}w_{1})[\alpha_{1},\alpha_{1},\alpha_{3}]+\order{t^{7}}.

(17)

Define $\beta_{5}=[\alpha_{1},\alpha_{1},\alpha_{3}]$ so

	$\displaystyle S_{2}(w_{1}t)S_{2}(w_{0}\tau)S_{2}(w_{1}t)$	$\displaystyle=\exp\{t(2w_{1}+w_{0})\alpha_{1}+t^{3}(2w_{1}^{3}+w_{0}^{3})\alpha_{3}+t^{5}(2w_{1}^{5}+w_{0}^{5})\alpha_{5}$
		$\displaystyle\quad+t^{5}\frac{1}{6}(w_{0}^{2}w_{1}^{3}-w_{1}^{2}w_{0}^{3}+w_{0}w_{1}^{4}-w_{0}^{4}w_{1})\beta_{5}+\order{t^{7}}\bigg{\}}.$		(18)

By an induction argument Yoshida shows that

\displaystyle S^{(m)}(t)=\exp\{tA_{1,m}\alpha_{1}+t^{3}A_{3,m}\alpha_{3}+t^{5}(A_{5,m}\alpha_{5}+B_{5,m}\beta_{5})+\order{t^{7}}\bigg{\}},

(19)

where $A_{j,m}$ and $B_{5,m}$ are polynomials on the variables $w_{0},\dots,w_{m}$ .

The case $m=0$ is just the symmetric BCH formula, so it is clear that Eq. (19) holds with

$\displaystyle A_{1,0}$	$\displaystyle=w_{0},$
$\displaystyle A_{3,0}$	$\displaystyle=w_{0}^{3},$
$\displaystyle A_{5,0}$	$\displaystyle=w_{0}^{5},$
$\displaystyle B_{5,0}$	$\displaystyle=0.$	(20)

To prove Eq. (19) for $m>0$ , one needs to show that the exponential is of the form with operator $\alpha_{1}$ for first order in $t$ , operator $\alpha_{3}$ for third order in $t$ , and operators $\alpha_{5}$ and $\beta_{5}$ for fifth order in $t$ . This result may be shown using

$\displaystyle S^{(m+1)}(t)$	$\displaystyle=S_{2}(w_{m+1}t)S^{(m)}(t)S_{2}(w_{m+1}t)$
	$\displaystyle=\exp\{tw_{m+1}\alpha_{1}+t^{3}w_{m+1}^{3}\alpha_{3}+t^{5}w_{m+1}^{5}\alpha_{5}+\order{t^{7}}\bigg{\}}$
	$\displaystyle\quad\times\exp\{tA_{1,m}\alpha_{1}+t^{3}A_{3,m}\alpha_{3}+t^{5}(A_{5,m}\alpha_{3}+B_{5,m}\beta_{5})+\order{t^{7}}\bigg{\}}$
	$\displaystyle\quad\times\exp\{tw_{m+1}\alpha_{1}+t^{3}w_{m+1}^{3}\alpha_{3}+t^{5}w_{m+1}^{5}\alpha_{5}+\order{t^{7}}\bigg{\}}$
	$\displaystyle=\exp\{2tw_{m+1}\alpha_{1}+tA_{1,m}\alpha_{1}+2t^{3}w_{m+1}^{3}\alpha_{3}+t^{3}A_{3,m}\alpha_{3}+2t^{5}w_{m+1}^{5}\alpha_{5}+t^{5}A_{5,m}\alpha_{5}+t^{5}B_{5,m}\beta_{5}$
	$\displaystyle\quad+\frac{1}{6}t^{5}(A_{1,m}^{2}w_{m+1}^{3}-A_{1,m}A_{3,m}w_{m+1}-w_{m+1}^{2}A_{3,m}+w_{m+1}^{4}A_{1,m})\beta_{5}+\order{t^{7}}\bigg{\}}.$	(21)

Hence, if the product formula can be expressed as in the form (19) for $S^{(m)}(t)$ , it can again be expressed in this form for $S^{(m+1)}(t)$ , establishing it for all $m\geq 0$ by induction. This expression also shows that the scalar coefficients can be determined from the formulae

$\displaystyle A_{1,m+1}$	$\displaystyle=2w_{m+1}+A_{1,m},$
$\displaystyle A_{3,m+1}$	$\displaystyle=2w_{m+1}^{3}+A_{3,m},$
$\displaystyle A_{5,m+1}$	$\displaystyle=2w_{m+1}^{5}+A_{5,m},$
$\displaystyle B_{5,m+1}$	$\displaystyle=B_{5,m}+\frac{1}{6}(A_{1,m}^{2}w_{m+1}^{3}-A_{1,m}A_{3,m}w_{m+1}-w_{m+1}^{2}A_{3,m}+w_{m+1}^{4}A_{1,m}).$	(22)

See Appendix A for an explanation of how to extend this approach to 10th order.

Constraints to satisfy in order to derive 6th order formula.

To derive a 6th order formula, the lower-order terms in the exponential in Eq. 19 must be zero (see also Eq. (5.16) of [20]), which gives the four conditions

A_{1,m}=1,\quad A_{3,m}=0,\quad A_{5,m}=0,\quad B_{5,m}=0.

(23)

For $m=3$ there are four unknowns $w_{0}$ to $w_{3}$ , and it can be expected there are solutions because there are the same number of equations as unknowns. In practice $A_{1,m}=1$ is satisfied by taking $w_{0}=1-2\sum_{j}w_{j}$ , so there are then three equations for three unknowns $w_{1},w_{2},w_{3}$ . Whereas it is possible to solve the equations using the Newton-Raphson method, the expression for the appropriate Jacobian matrix is complicated, so Yoshida instead uses the Brent method. Yoshida produced three $m=3$ solutions in this way, and states “It seems that there is no other solution.” We have performed an extensive search and also found no more solutions.

The product formulae obtained through the Yoshida method also work for exponentials of sums of more operators. The $S_{2}$ product formula can again be used as the building block for the product formula, and we can write

S_{2}(t)=\exp\left(\sum_{\ell=0}^{\infty}\tilde{\alpha}_{\ell}t^{\ell}\right),

(24)

where $\tilde{\alpha}_{\ell}$ are now order- $\ell$ multicommutator expressions on the $J$ terms. The reasoning for finding the product formula is entirely based on the construction with ${\alpha}_{\ell}$ , but without taking advantage of its particular form, so exactly the same reasoning holds for $\tilde{\alpha}_{\ell}$ . This immediately implies that the higher-order product formulae work in general. This is an advantage of constructing product formulae as products of $S_{2}$ , because deriving coefficients separately for exponentials of $X$ and $Y$ would not generalise.

II.3 Processed product formulae

Another technique to obtain higher-order product formulae is that of processing [24, 25, 26, 27, 28, 29, 30]. In this technique a product formula $S_{k}$ of order $k$ is generated by the composition of a kernel $\Sigma$ conjugated by a processor operator $P$ as

\displaystyle S_{k}=P\Sigma P^{-1}.

(25)

The advantage of this method is that usually $\Sigma$ requires fewer stages than other methods, and due to the construction, we have that $S_{k}^{n}=P\Sigma^{n}P^{-1}$ . This means that the cost of using the product formula is effectively that of the kernel in the usual application where a long time evolution is partitioned into many repetitions of the product formula for short times.

Typically, one chooses a kernel $\Sigma$ as a product formula that has a smaller order than $k$ , but conjugating by processors gives an order $k$ integrator. A table of available kernels that can be combined with processors is given in Ref. [24], and recently new kernels were reported in [31]. Another advantage of this method is that the number of nonlinear equations required to be solved to find new product formulae is reduced, as one only needs to find new parameters for the kernel $\Sigma$ , and then conjugate it by known processors in the literature to give an order $k$ integrator.

To be more specific, Ref. [24] gives a set of conditions for the kernels and processors in Table 4 of that work. Reference [24] uses the notation $Y_{j}$ rather than $\alpha_{j}$ , so $S_{2}(t)=e^{tY_{1}+t^{3}Y_{3}\cdots}$ . It then uses the notation $E_{j,i}$ for multicommutator expressions, with $E_{j,1}=Y_{j}\equiv\alpha_{j}$ and for example

	$\displaystyle E_{5,2}$	$\displaystyle=[Y_{1},Y_{1},Y_{3}]\equiv\beta_{5}\,,$		(26)
	$\displaystyle E_{7,2}$	$\displaystyle=[Y_{3},Y_{1},Y_{3}]\equiv-\gamma_{7}\,,$		(27)

where $\gamma_{7}$ is given in Appendix A. In Ref. [24] the quantities $f_{j,i}$ are the coefficients of $E_{j,i}$ , so the equivalent of Eq. (19) in that notation is

\displaystyle S^{(m)}(t)=\exp\!\bigg{\{}tf_{1,1}E_{1,1}+t^{3}f_{3,1}E_{3,1}+t^{5}(f_{5,1}E_{5,1}+f_{5,2}E_{5,2})+\order{t^{7}}\bigg{\}}\,.

(28)

To obtain the conditions for a 6th order kernel from Ref. [24], one should use the conditions for $f_{j,i}$ up to $q=5$ in Table 4, which give $f_{1,1}=1$ , $f_{3,1}=0$ , $f_{5,1}=0$ , which in the above notation are equivalent to

A_{1,m}=1,\quad A_{3,m}=0,\quad A_{5,m}=0.

(29)

The conditions for the kernel are the same as for the complete product formula, but missing $B_{5,m}=0$ . That is a general feature of the kernel order conditions; they are a subset of the conditions for the complete product formula, enabling shorter kernels, or more flexibility to choose kernels of the same length but lower error. For 4th order there are the same number of conditions, so the kernel is a valid product formula of that order.

For the processors, Ref. [24] uses the notation $p_{j,i}$ for the coefficients of $E_{j,i}$ in the expansion of $P$ . These need to be determined from the BCH formula rather than symmetric BCH, and instead of only odd-order terms there need to be only even-order terms. Some symmetries mean that low-order odd terms automatically cancel, but higher-order ones need to be made zero by the choice of $P$ . Table 4 in Ref. [24] gives the conditions on the even-order terms. The first processor condition depending on $f_{j,i}$ is for $q=5$ , which means that for 6th order and higher the processors are dependent on the kernel. The 4th order the kernel is already 4th order, and the processors will yield another 4th order product formula, but can be chosen to reduce the error (which arises from higher-order terms). See Appendix B for more detail on the method of determining processors.

III Method of solution

III.1 Numerical methods

The size of the product formula we consider as in Eq. (15) is giverned by the value of $m$ , so there is a product of $M=2m+1$ integrators. The number of free parameters is $m$ , and it needs to be a least a large as the number of independent equations in order for there to be a solution. Yoshida [20] finds that $m=7$ is minimal for 8th order, and similarly to prior work [22] we find that $m=15$ is minimal for 10th order; see Appendix A. Choosing the minimal $m$ typically yields product formulae with large error, and for that reason we also consider larger values of $m$ .

For the numerical solution of the simultaneous nonlinear equations, we found that Matlab’s fsolve (using the Levenberg-Marquardt algorithm) provided rapid results using the vector of errors. The method was to repeatedly solve from random starting vectors $\vec{w}$ generated according to the normal distribution. The best solutions were those with smaller values for the coefficients; for much of the calculations we selected standard deviations of $0.6$ for 8th order and $0.9$ for 10th order. We filtered the large number of solutions by numerically checking the error for example Hamiltonians, and selecting those that provided the best performance. We then perform tests for larger numbers of samples, to better select the best performing product formulae. We also further refined the solutions to give the results to extended precision, which enables testing of the error with smaller values of $t$ . We were also able to further reduce the error by using an alternative solution method using a Taylor series.

III.2 Solution using Taylor expansion

Another solution method we use is based on computing the Taylor expansion of both the exact exponential and its product formula approximation with given $w_{j}$ . This Taylor expansion is performed on both sides up to the desired order of approximation for the integrator. By imposing equality on terms of the same order we obtain a series of equations for $w_{j}$ which can be solved. For higher orders, a large number of simultaneous equations are obtained, so we need an automated way of generating them.

To make precise the problem we are solving, denote as $\mathcal{T}_{k}[\cdot]$ the map giving the Taylor expansion in $t$ around $0$ , truncated at order $k$ , so

	$\displaystyle\mathcal{T}_{k}[e^{t(X+Y)}]$	$\displaystyle=\sum_{p=0}^{k}\frac{t^{p}}{p!}(X+Y)^{p}$
		$\displaystyle=\sum_{p=0}^{k}\frac{t^{p}}{p!}\sum_{r_{1},\cdots,r_{p}=0}^{1}X^{r_{1}}Y^{1-r_{1}}\cdots X^{r_{p}}Y^{1-r_{p}}.$		(30)

We consider a sum of two operators $X+Y$ , but this approach for solving for $w_{j}$ will also be sufficient to provide product formulae for sums of arbitrary numbers of terms. That is because the solutions must also be solutions of the equations derived using Yoshida’s method, and as explained above those equations will be the same when considering sums of arbitrary numbers of terms. Now consider the ansatz operator of Yoshida from Eq. 15

$\displaystyle S^{(m)}(t,w_{1},\cdots,w_{m})$	$\displaystyle=e^{tw_{m}X/2}e^{tw_{m}Y}e^{t(w_{m}+w_{m-1})X/2}e^{tw_{m-1}Y}e^{t(w_{m-1}+w_{m-2})X/2}$
	$\displaystyle\quad\cdots e^{tw_{0}Y}e^{t(w_{1}+w_{0})X/2}e^{tw_{1}Y}\cdots e^{tw_{m}X/2}$
	$\displaystyle=e^{tc_{1}X}e^{tc_{2}Y}\cdots e^{tc_{4m+3}X},$	(31)

where in the last line we have defined the constants $c_{1}=w_{m}/2$ , $c_{2}=w_{m}$ , $c_{3}=(w_{m}+w_{m-1})/2,\cdots$ $c_{4m+3}=w_{m}/2$ . We Taylor expand this ansatz up order $k$ , noting that the total number of exponentials in Yoshida’s ansatz is $4m+3$

\displaystyle\mathcal{T}_{k}[S^{(m)}(t,w_{1},\dots,w_{m})]

\displaystyle=\sum_{\begin{subarray}{c}r_{1},\cdots,r_{4m+3}=0\\ r_{1}+\dots+r_{4m+3}\leq k\end{subarray}}^{k}\frac{t^{r_{1}+\dots+r_{4m+3}}}{r_{1}!\cdots r_{4m+3}!}c_{1}^{r_{1}}\dots c_{4m+3}^{r_{4m+3}}X^{r_{1}}Y^{r_{2}}\dots X^{r_{k-1}}Y^{r_{4m+3}}.

(32)

We require that the product formula obtained from our solution procedure be independent of the choice of $X$ and $Y$ , so need to match the coefficients for each sequence of products of $X$ and $Y$ . Because $X$ and $Y$ are non-commuting, we need to record coefficients for each ordered sequence. To do this in an automated way we construct a data structure to store the coefficients.

Given operators of the form $e^{cX}=I+cX+\frac{c^{2}}{2!}X^{2}+\frac{c^{3}}{3!}X^{3}+\cdots$ and $e^{dY}=I+dY+\frac{d^{2}}{2!}Y^{2}+\frac{d^{3}}{3!}Y^{3}+\cdots$ with $c,d$ scalar coefficients and $X,Y$ general operators, we describe this Taylor expansion up to order $k$ using an array encoding. First, write monomials composed of $X$ and $Y$ operators in lexicographical order and note that these operators can be mapped to binary numbers:

\begin{matrix}I&X&Y&XX&XY&YX&YY&XXX&XXY&\cdots\\ 1&10&11&100&101&110&111&1000&1001&\cdots\\ 1&2&3&4&5&6&7&8&9&\cdots\end{matrix}

(33)

To construct a bit string, we map each $X$ to 0 and each $Y$ to 1, then place a 1 on the left to flag the length of the string, as shown in the second line of Eq. 33. Then, to obtain the operator products, simply remove the leading 1 and then map $0$ to $X$ and $1$ to $Y$ . The empty string corresponds to the identity $I$ . Now, to store the coefficients in a sum of products of $X$ and $Y$ , convert each product to a binary integer as above, then store the coefficient in the corresponding location in a vector.

By way of illustration, consider the polynomial $10I+3X+2Y+2X^{2}+YX$ . This operator would be stored in an array as $[10,3,2,2,0,1,0,\cdots]$ . In this way, any polynomial of $X$ and $Y$ can be efficiently stored in a vector. We denote this vector storing the coefficients of operators of order up to $k$ as $\mathrm{vec}_{k}[p(X,Y)]$ , where $k$ denotes that the vector will only store the coefficients of the corresponding operators up to order $k$ (so a vector of length $2^{k+1}-1$ ) and $p(X,Y)$ is the polynomial in terms of $X$ and $Y$ .

This approach can be used to solve for product formulae, but is also very effective at further improving the performance of solutions that have already been found. If a solution is of 8th order (for example), then we can consider the error using the Taylor expansion up to and including the 9th order. This will give nonzero error for the 8th order product formula, so minimising the error gives a new product formula. That can be used as a starting point to again solve for an 8th order formula, and often it is found that the product formula adjusted in that way has reduced error.

IV Product formula search

IV.1 Improved 8th order

Following the method outlined in Section III, we have searched for new product formulae and found thousands of product formulae of 8th order. The case $m=7$ yields product formulae of minimal length, because the number of equations is equal to the number of unknowns. We have found over 600 solutions, and the search now finds almost only repeated solutions and very few new solutions. This indicates that we have found nearly all solutions, but it is also possible that there are many more solutions with large values of $\vec{w}$ . Indeed, the most recent new solutions we found have significantly larger values of $\vec{w}$ . We find that large values of $\vec{w}$ correspond to worse product formulae with larger error. Therefore, even if there are many more solutions with large $\vec{w}$ , they likely will not give improved performance over those we have already found.

Once we have obtained the potential solutions, we generate random Hamiltonians and compute the product formula errors as a function of time. For 8th order product formulae we know that the product formula error is $\order{t^{9}}$ . We check the error scaling by picking two times $t_{1}$ and $t_{2}$ and computing errors $\delta_{1}$ and $\delta_{2}$ at these times, then we compute $\log(\delta_{1}/\delta_{2})/\log(t_{1}/t_{2})$ and check that it is close to $9$ . As the error for our product formulae is given by $\delta(t)=\chi t^{9}$ where $\chi$ is a constant factor, we can compute $\chi$ for each of them by from $\chi=\delta/t^{9}$ . For each product formula, we compute an average constant factor error; this average corresponds to the geometric mean of the constant factors computed for each random Hamiltonian. This method of comparing the performance of product formulae through the estimation of the constant factor in the error has been used before (see for example [33]) and is considered a good approximation to the performance of the product formula in practice.

The best performing 8th order product formula of minimal length $m=7$ is shown in Table 1. This has average constant factor $\chi=5.8\times 10^{-6}$ , where we are measuring error using the spectral norm. We evaluated all the 8th order solutions of Yoshida, and found Solution D was best, with a constant $\chi=9.7\times 10^{-4}$ . See Section VI for more extensive comparisons other product formulae. The best solution we found for $m=7$ corresponds to that reported as s15odr8 in [21], where the relation between the coefficients given is $w_{j}=\delta_{8-j}$ (so they list the coefficients in reverse order to ours). The extensive nature of our search indicates that the product formula reported in Ref. [21] is optimal for minimal-length 8th order product formulae.

	Best 8th order with $m=7$	Best 8th order with $m=8$
$w_{1}$	$0.315293092396766596632056663811$	$0.29137384767986663096528500968049$
$w_{2}$	$0.33462491824529818378495797988218$	$0.26020394234904150277316667709864$
$w_{3}$	$0.2990641813036559238444635406886$	$0.18669648149540687549831902999911$
$w_{4}$	$-0.57386247111608226665638772663554$	$-0.40049110428180105319963667975074$
$w_{5}$	$0.19075471029623837995387625645037$	$0.15982762208609923217390166127256$
$w_{6}$	$-0.40910082580003159399730009589356$	$-0.38400573301491401473462588779099$
$w_{7}$	$0.74167036435061295344822780178381$	$0.56148845266356446893590729572808$
$w_{8}$	NA	$0.12783360986284110837857554950443$

Table 1: Our best-performing 8th order solutions with

m=7

and 8.

	Best 8th order for spectral-norm error	Best 8th order for eigenvalue error
$w_{1}$	$0.59358060400850625863514059265224$	$0.10467636532245895252340732579853$
$w_{2}$	$-0.46916012347004197296293264921328$	$-0.57896999331780988041471955125778$
$w_{3}$	$0.2743566425898467907228242878146$	$0.57503350160061785946141563279891$
$w_{4}$	$0.17193879484656773059919074965377$	$0.12231011868707029786561397542663$
$w_{5}$	$0.23439874482541384415430578747541$	$0.27793149999039524816733903301747$
$w_{6}$	$-0.48616424480326193899617759997914$	$-0.37349605088056728482635987352576$
$w_{7}$	$0.49617367388114660354871757044906$	$0.11575566589480463220616543972403$
$w_{8}$	$-0.32660218948439130114501815323814$	$0.1464645610975800618712569230326$
$w_{9}$	$0.23271679349369857679445410270557$	$-0.39443578322284085764474498594073$
$w_{10}$	$0.098249557414708533273471906180643$	$0.44370228726021218923197141183196$

Table 2: Our best-performing 8th order solutions when setting

m=10

We have also conducted a search for 8th order solutions with $m=8$ , finding over 600 solutions. Using $m=8$ results in an underdetermined system of equations with continuous sets of solutions, and nearly half of these solutions were close (with a norm less than 0.01). The underdetermined nature of the solutions also gives the flexibility to adjust the solution to reduce the error by further using the method described in Section III. The best solution found is given in Table 1, with an average constant factor of $\chi=5.7\times 10^{-7}$ , which is an order of magnitude improvement with only a slight increase in the number of exponentials. This solution is close to solutions s17odr8a and s17odr8b given in Ref. [21], suggesting that the difference is due to the system of equations being underdetermined.

Extending our search further to $m=10$ we found over 100,000 solutions, including solutions improving over those in prior work. Our best solution is given in Table 2 and has a constant of $\chi=4.9\times 10^{-8}$ . This new 8th order product formula for $m=10$ improves 40,000 times over the ones provided by Yoshida. An alternative measure of the performance of the product formulae is the accuracy of the eigenvalues. When considering this error, the product formula has a somewhat smaller constant factor $\zeta=1.1\times 10^{-8}$ . It is also possible to select product formulae to provide the minimum error in the eigenvaules. When we do this, we obtain the second solution given in Table 2, which has the significantly smaller constant factor in the error $\zeta=1.6\times 10^{-9}$ . That is better than any prior product formulae we have tested when accounting for the number of exponentials.

A further improvement in product formulae can be obtained by using processing [24]. We have performed a search for high-accuracy 8th order product formulae using this method with $m=8$ , and found one particular example that provides excellent performance. In our testing this solution is found repeatedly, and provides significantly lower error than any other solutions. Because the equations are underdetermined, the solution may be adjusted to provide best performance.

	Best 8th order kernel	Processor for kernel
$w_{1}$	$0.21784176681731006074681969186513$	$-0.44324901019570126590495430949294$	$\gamma_{1}$
$w_{2}$	$0.1947017706053903224022456342907$	$0.25459857192003772850622377066944$	$\gamma_{2}$
$w_{3}$	$0.18372413281145589944261642180363$	$-0.73862036266779261573694538099739$	$\gamma_{3}$
$w_{4}$	$-0.37307499512657736825709230652023$	$-0.00024139614958652134370419495289618$	$\gamma_{4}$
$w_{5}$	$0.15757644257569146373033662060461$	$0.73873460354125365739379753874964$	$\gamma_{5}$
$w_{6}$	$-0.33342207567391682979227850551172$	$-0.20285971152536085519251666906017$	$\gamma_{6}$
$w_{7}$	$0.51788649682987924281787142226803$	$0.44989521689676869571827637424046$	$\gamma_{7}$
$w_{8}$	$0.21456475499897766986381219621761$	$0.29538398007876871184026747505657$	$\gamma_{8}$
		$-0.3364996155865700091428329802017$	$\gamma_{9}$

Table 3: Our best-performing 8th order kernel and processor when setting

m=8

. The value of

\gamma_{10}

for the processor is chosen such that the sum is zero.

Our best solution is given in Table 3. For the kernel there are many solutions for the processor, but we list the one with the best performance of the solutions we have found. This one was found by solving for many (hundreds of thousands) of processors, picking the one with the best performance, and further optimising. This kernel provides a constant of $\zeta=2.2\times 10^{-9}$ , only slightly larger than the best length $m=10$ product formula, while providing better performance due to the shorter length. Moreover, using the processor the constant for the spectral-norm error is $\chi=5.4\times 10^{-8}$ . That is only about 10% more than the best length $m=10$ formula, while again providing better performance due to the shorter length.

IV.2 Finding 10th order

We have also used our solution procedure to find new $10$ th order product formulae. We generalised Yoshida’s method, and find $15$ independent equations to be solved (see Appendix A). We performed searches for solutions with $m=15$ (the minimal number) to $m=18$ . Again the larger values of $m$ give the flexibility to adjust the solution to reduce the error. We report our 10th order product formulae for $m=15$ , $16$ and $18$ in Table 4 that are best for spectral-norm error. We don’t report values for $m=17$ as the performance of the product formulae found were worse than for $m=16$ . In Table 5 we report product formulae for $m=17$ and $m=18$ with our best results for eigenvalue error.

As in Section IV.1, we compare the performance of product formulae of $10$ th order by computing the constant factors $\chi$ and $\zeta$ for random Hamiltonians. For the best solution with $m=16$ we have a constant factor of $\chi=1.9\times 10^{-8}$ , and the best solution with $m=15$ has $\chi=4.5\times 10^{-7}$ , which is about a factor of 24 times worse. The far better constant factor for $m=16$ is far more significant than the slightly larger number of exponentials in the product formula.

Further increasing $m$ to 18 gives a product formula with $\chi=3.1\times 10^{-9}$ ; a significant further improvement. Again we can consider the error in terms of the eigenvalues; in this case the constant for eigenvalue error is $\zeta=2.2\times 10^{-9}$ , only slightly less. Selecting the best performing product formula for eigenvalue error gives the significantly smaller constant $\zeta=4.2\times 10^{-10}$ . However, out of our solutions the best eigenvalue performance is given by a $m=17$ solution, which yields $\zeta=1.8\times 10^{-10}$ . These solutions are listed in Table 5.

In the search for 10th order product formulae, unlike in the case of 8th order, we find that almost all new solutions found are different from those found before. That indicates that there is an extremely large number of solutions, and we have only found a very small proportion of them. Indeed, the solutions are distinct from those in prior work.

	Best 10th order solution with $m=15$	Best 10th order solution with $m=16$	Best 10th order solution with $m=18$
$w_{1}$	$0.14552859955499429739088135596618$	$-0.4945013179955571856347147977644$	0.019042478645106035261914181501875
$w_{2}$	$-0.48773512068133537309419933740564$	$0.2904317222970121479878414292093$	$-0.48337326409346903272186302946692$
$w_{3}$	$0.12762011242429535909727342301656$	$0.34781541068705330937913890281003$	$0.035060961741879192451298102625219$
$w_{4}$	$0.70225450019485751220143080587959$	$-0.98828132118546184603769781410676$	$0.20690475331505992081884048319725$
$w_{5}$	$-0.62035679146761710925756521405042$	$0.98855187532756405235733957305613$	$0.039554342269800383312212959879587$
$w_{6}$	$0.39099152412786178133688869373114$	$-0.34622976933123177430694714630668$	$0.062010837356401048997119918637392$
$w_{7}$	$0.17860253604355465807791041367045$	$0.20218952619073117554714280367018$	$-0.46961231983086041266381539270133$
$w_{8}$	$-0.80455783177921776295588528272593$	$0.13064273069786247787208895471461$	$-0.15137223243888068391593992998235$
$w_{9}$	$0.053087216442758242118687385646283$	$-0.26441199183146805554735845490359$	$0.13186222745709395576675594763784$
$w_{10}$	$0.86836307910275556258687030904753$	$0.060999140559210408869096992291531$	$0.44628663303136375145122785014895$
$w_{11}$	$-0.85326297197907834671536254437991$	$-0.6855442489606141359108973267028$	$-0.31721379667717916478350053562451$
$w_{12}$	$-0.11732457198874083224967699358383$	$-0.15843692473786584550599206557006$	$0.44313588649776693705154231063871$
$w_{13}$	$0.03827345494186056632406947772047$	$0.15414691779958299150286452215575$	$0.16887007584153091511395119434171$
$w_{14}$	$0.74843529029532498233997793305357$	$0.66715205827214320371061839297055$	$-0.22652658662557993653900899346103$
$w_{15}$	$0.30208715621975773712410948025906$	$0.20411874474696598289603677693511$	$0.13053736297137232483181427384048$
$w_{16}$	NA	$0.081207318210272593225087711441684$	$0.11337301050285651053819309187802$
$w_{17}$	NA	NA	$0.056199557660148108798028960238124$
$w_{18}$	NA	NA	$0.038918323115794012069868989863952$

Table 4: Our best performing 10th order solutions for spectral-norm error with

m=15

m=16

and

m=18

	Best 10th order solution with $m=17$	Best 10th order solution with $m=18$
$w_{1}$	$-0.28371232689144296279654621726493$	$0.025722554623006480493726308396586$
$w_{2}$	$0.046779504778147381605331000278223$	$0.024673923089392154535100643510344$
$w_{3}$	$0.36845892382797770619657504217539$	$-0.40545153312882551694596948883526$
$w_{4}$	$0.19186204094674514739760408197461$	$0.086870323364257282181073061915168$
$w_{5}$	$-0.53123134392680669702873064192428$	$0.12368899347772019656137276541942$
$w_{6}$	$-0.0081253242720827266680816105600661$	$0.34599591069083361101791099618656$
$w_{7}$	$-0.16389450414378567860032917538393$	$0.046765678517740550705548061486811$
$w_{8}$	$0.18514766119291405032528647881$	$-0.27103335145245847800657868572535$
$w_{9}$	$0.5383584694754681989174668806505$	$0.13398594471200943261255065567866$
$w_{10}$	$-0.30583981835573485697292316732177$	$-0.45010365706956744617357917877887$
$w_{11}$	$0.43199935609523301289295473774488$	$0.33699858113023399397587906362881$
$w_{12}$	$0.1510502301631786853020124612813$	$0.14286479024077276505929263927029$
$w_{13}$	$-0.35051099204829676098801520498121$	$-0.30679647776174213774450994020067$
$w_{14}$	$0.1032971125844291674511513007661$	$0.048785861198921384322572380948858$
$w_{15}$	$0.15043936943817152697371946806229$	$0.035258483631052620304882207189439$
$w_{16}$	$0.12118469498650736511410491586846$	$-0.22380268023236595677874655821875$
$w_{17}$	$0.10437742779547826358296681557444$	$0.42346449759412505872094526232433$
$w_{18}$	NA	$0.14888705463805455702454629353763$

Table 5: Our best performing 10th order solutions for eigenvalue error with

m=17

and

m=18

V Methods for comparison of product formulae

V.1 Same order comparison

We make a fair comparison between product formulae of different length in the following way. An order $k$ integrator for time $t$ will have an error $\delta=\chi t^{k+1}$ where $\chi$ is a real constant. Let $T$ be the total evolution time for an integrator of order $k$ , and $\varepsilon$ be the maximum allowable error. Subdivide the evolution time $T$ into $r$ subintervals, so $t={T}/{r}$ is the length of each time subinterval. We thus have $\chi\left({T}/{r}\right)^{k+1}\approx{\varepsilon}/{r}$ , which gives

\displaystyle r\approx\left(\frac{\chi T}{\epsilon}\right)^{{1}/{k}}T.

(34)

As explained above, the number of exponentials in the product is $(4m+2)(J-1)+1$ . When applying products of these product formulae, two exponentials can be combined, so the effective number for each is $(4m+2)(J-1)$ . As a result, the total number of exponentials can be given as proportional to

\displaystyle M\left(\frac{\chi T}{\epsilon}\right)^{{1}/{k}}T\,,

(35)

where we have ignored a common factor of $2(J-1)$ , and $M=2m+1$ . If we wish to compare product formulae of the same order, then we need only compare the values of $M\chi^{1/k}$ , and the one with the smaller value is the more efficient product formula. Similarly, if we consider eigenvalue error then we should compare $M\zeta^{1/k}$ between the formulae.

V.2 Thresholds for different order

If we wish to compare product formulae of different order, then we need to take account of the values of $T$ and $\varepsilon$ . Assume we have two integrators of order $k_{1}$ and $k_{2}$ , with corresponding constants $\chi_{1}$ , $\chi_{2}$ . Given $T$ and $\varepsilon$ , when the two integrators use the same number of exponentials we have $M_{1}r_{1}=M_{2}r_{2}$ , thus

	$\displaystyle M_{1}\left(\frac{\chi_{1}T}{\varepsilon}\right)^{{1}/{k_{1}}}T$	$\displaystyle=M_{2}\left(\frac{\chi_{2}T}{\varepsilon}\right)^{{1}/{k_{2}}}T$		(36)
	$\displaystyle\implies\quad\frac{T}{\varepsilon}$	$\displaystyle=\left(\frac{M_{2}\chi_{2}^{1/k_{2}}}{M_{1}\chi_{1}^{1/k_{1}}}\right)^{\frac{1}{\frac{1}{k_{1}}-\frac{1}{k_{2}}}}.$		(37)

This gives the threshold beyond which the higher-order product formula should be used for improved performance. Again, when considering eigenvalue error, $\chi$ would be replaced with $\zeta$ .

A limitation of this analysis it is assumed that the time step $(\epsilon/(\chi T))^{1/k}$ is small enough for the scaling law for the error to hold. To adjust the threshold for cases where the time step is large, we can consider a more general functional dependence of the error on the time interval as $f(t)$ . Then we require for time $T$ and error $\epsilon$ that

f(t)\times T/t=\epsilon\,.

(38)

The cost is then proportional to $M\times r=MT/t$ . For the threshold between two product formulae, we require

M_{1}T/t_{1}=M_{2}T/t_{2}\,.

(39)

That implies $t_{2}=t_{1}M_{2}/M_{1}$ . For the total error to be equal to $\epsilon$ in each case, we then require

f_{1}(t_{1})\times T/t_{1}=f_{2}(t_{2})\times T/t_{2}\,.

(40)

That can be rearranged to give

f_{1}(t_{1})=f_{2}(t_{1}M_{2}/M_{1})\times M_{1}/M_{2}\,.

(41)

This is an expression we can solve for $t_{1}$ , which may then be used to determine the $T/\epsilon$ threshold from

T/\epsilon=t_{1}/f_{1}(t_{1})\,.

(42)

This means that the threshold is still for the ratio $T/\epsilon$ , rather than having separate dependence on $T$ and $\epsilon$ .

V.3 Motivation for considering eigenvalue error

A further question is the measure of the error to be used. If the goal is to estimate eigenvalues of the Hamiltonian (as is often the case in quantum chemistry), then the measure of the error should be that in the eigenvalues. That is, how close are the eigenvalues of the product formula to those of the exact Hamiltonian evolution. In the case that the goal is to accurately reproduce the final state after the evolution, then an appropriate measure of the error should be the spectral norm of the difference of the unitary operators. This error accounts for both error in the eigenvalues and basis.

That error will upper bound the 2-norm error in the generated quantum state, but using the triangle inequality to bound the error for long evolution times overestimates the error. This is because the error beyond that in the eigenvalues cancels when using a product of many short time intervals. Let us denote the exact evolution operator for short time $U$ , and the approximate evolution operator provided by the product formula as $\widetilde{U}$ . In the basis of the Hamiltonian’s eigenstates, $U$ is diagonal, and we can diagonalise $\widetilde{U}$ in this basis as $\widetilde{U}=VDV^{\dagger}$ . Then the difference between $D$ and $U$ describes the eigenvalue error, since these are diagonal matrices with the eigenvalues on the diagonal. The matrix $V$ describes the basis error.

It is convenient to write $V=e^{i\tilde{H}}$ for some Hermitian matrix $\tilde{H}$ . Then the error can be given as

$\displaystyle\\|\widetilde{U}-U\\|$	$\displaystyle=\\|e^{i\tilde{H}}De^{-i\tilde{H}}-U\\|$
	$\displaystyle=\\|(I+i\tilde{H}+\mathcal{O}(\tilde{H}^{2}))D(I-i\tilde{H}+\mathcal{O}(\tilde{H}^{2}))-U\\|$
	$\displaystyle=\\|i[\tilde{H},D]+D-U+\mathcal{O}(\tilde{H}^{2})\\|$
	$\displaystyle\leq\\|[\tilde{H},D-I]\\|+\\|D-U\\|+\mathcal{O}(\tilde{H}^{2})\,.$	(43)

That is, the error can be split up between a part $[\tilde{H},D-I]$ corresponding to basis error, and a part $D-U$ corresponding to eigenvalue error. Note that $D-I$ is proportional to $t$ , so the error in the operator due to the basis error is one order higher than the error in the basis.

Then, for the case of a large number of time steps $r$ , we can bound the overall error as

$\displaystyle\\|\widetilde{U}^{r}-U^{r}\\|$	$\displaystyle=\\|VD^{r}V^{\dagger}-U^{r}\\|$
	$\displaystyle=\\|(VD^{r}V^{\dagger}-D^{r}V^{\dagger})+(D^{r}V^{\dagger}-U^{r}V^{\dagger})+(U^{r}V^{\dagger}-U^{r})\\|$
	$\displaystyle\leq\\|V-I\\|+\\|D^{r}-U^{r}\\|+\\|V^{\dagger}-I\\|$
	$\displaystyle\leq 2\\|V-I\\|+r\\|D-U\\|\,.$	(44)

That is, the error due to the eigenvalues is multiplied by a factor of $r$ , but the error due to the basis is not. The error in the basis can be larger than for a single time step, but only by a factor of $\sim 1/t$ . That factor arises because for a single time step the basis error is commuted with $D-I$ which is proportional to $t$ , so the expression here without $D-I$ can be a factor $\propto 1/t$ larger.

In practice we find that in many cases the error in the basis scales as $\mathcal{O}(t^{k})$ , resulting in the expected contribution to the spectral-norm error $\mathcal{O}(t^{k+1})$ . However, in other cases we find scaling of the error in the basis as $\mathcal{O}(t^{k+1})$ , resulting in a higher-order contribution to the spectral-norm error. Let us give the error in the eigenvalues as $\zeta t^{k+1}$ , and the error in the basis as $\mu t^{k+\nu}$ where $\nu\in[0,1]$ to account for the various scalings found numerically.

Then the choice of $r$ to make the error in the eigenvalues smaller than $\epsilon$ is

r\approx\left(\frac{\zeta T}{\epsilon}\right)^{1/k}T.

(45)

The contribution to the total error from the basis error is

2\mu(T/r)^{k+\nu}\approx 2\mu\left(\frac{\epsilon}{\zeta T}\right)^{1+\nu/k}\,.

(46)

That gives the ratio of the basis error to the eigenvalue error as approximately

\frac{2\mu}{\epsilon}\left(\frac{\epsilon}{\zeta T}\right)^{1+\nu/k}=\frac{2\mu}{\zeta T}\left(\frac{\epsilon}{\zeta T}\right)^{\nu/k}\,.

(47)

Hence, for large $T$ the eigenvalue error should be dominant. This means that the eigenvalue error for a single step is the relevant measure to estimate the spectral-norm error for long-time evolution.

VI Numerical comparison of product formulae

VI.1 Comparison for general matrices

In what follows, we report on the comparison between product formulae of 4th, 6th, 8th and 10th order based on the threshold provided above. We analyse those found in this work, and those from prior work including processed product formulae. When giving the number of stages of a processed product formula, we refer to the number of stages in the kernel (not $P$ ), because those are the ones that would be repeated for a simulation over an extended time.

We list 4th order product formulae in Table 6, 6th order in Table 7, 8th order in Table 8, and 10th order in Table 9, and give the constant factors in the error scaling $\chi,\zeta$ (spectral norm or eigenvalue error). These product formulae include some constructed with even $M$ rather than $M=2m+1$ as for our product formulae. We generated 10,000 random Hamiltonians for each product formula, and then estimated the constants $\chi$ and $\zeta$ by taking the geometric mean. Each random Hamiltonian is generated as a sum of two random Hamiltonians $H=A+B$ , where $A$ and $B$ are random Hermitian matrices of dimension $6\times 6$ and norm 1.

label	$M$	processing	reference	$\chi$	$M\chi^{1/k}$	$\zeta$	$M\zeta^{1/k}$
S4m1	$3$	no	first Suzuki product [1]	$4.5\times 10^{-2}$	1.38	$3.0\times 10^{-2}$	1.25
S4m2	$5$	no	second Suzuki product [2]	$2.6\times 10^{-3}$	1.13	$4.2\times 10^{-4}$	0.72
O4M5	$5$	no	Ostmeyer Eq. (40) of [34]	$2.9\times 10^{-4}$	0.65	$1.2\times 10^{-4}$	0.52
BM4M6	$6$	no	$S_{6}$ Table 2 of [35]	$1.5\times 10^{-4}$	0.67	$3.6\times 10^{-5}$	0.47
BCE4m6b	$12$	no	BM ${}_{6}^{[4]}$ Table 1 of [31]	$4.5\times 10^{-5}$	0.99	$2.8\times 10^{-5}$	0.87
PPBCM4m6	$12$	yes	$P_{6}4$ Table 5 of [24]	$5.0\times 10^{-5}$	1.01	$1.4\times 10^{-5}$	0.74
BCE4m3	$6$	yes	$\psi_{3}^{[4]}$ Table 6 of [31]	$1.5\times 10^{-2}$	2.09	$1.8\times 10^{-3}$	1.23
BCE4m4	$8$	yes	$\psi_{4}^{[4]}$ Table 6 of [31]	$5.0\times 10^{-4}$	1.19	$1.3\times 10^{-4}$	0.85
BCE4m5	$10$	yes	$\psi_{5}^{[4]}$ Table 6 of [31]	$5.9\times 10^{-5}$	0.88	$2.7\times 10^{-5}$	0.72
BCE4m6	$12$	yes	$\psi_{6}^{[4]}$ Table 6 of [31]	$2.3\times 10^{-5}$	0.83	$8.9\times 10^{-6}$	0.66
BCE4m7	$14$	yes	$\psi_{7}^{[4]}$ Table 6 of [31]	$1.2\times 10^{-5}$	0.82	$3.9\times 10^{-6}$	0.62
BCE4m8	$16$	yes	$\psi_{8}^{[4]}$ Table 6 of [31]	$6.8\times 10^{-6}$	0.82	$2.0\times 10^{-6}$	0.60
BCE4m9	$18$	yes	$\psi_{9}^{[4]}$ Table 6 of [31]	$4.5\times 10^{-6}$	0.83	$1.1\times 10^{-6}$	0.58

Table 6: The list of 4th order product formulae. The column labelled “processing” indicates whether the formula uses processor. The label is the name we will refer to the formula by. The constant factor in the error is denoted

\chi

for spectral-norm error and

\zeta

for eigenvalue error, and the corresponding quantities

M\chi^{1/k}

and

M\zeta^{1/k}

are given. The best results are highlighted in bold. Results for product formulae that are not products of

S_{2}

are given in italics.

label	$M$	processing	reference	$\chi$	$M\chi^{1/k}$	$\zeta$	$M\zeta^{1/k}$
S6m1	$9$	no	first Suzuki product [1]	$4.0\times 10^{-2}$	5.26	$3.2\times 10^{-2}$	5.08
S6m2	$25$	no	second Suzuki product [2]	$1.0\times 10^{-5}$	3.68	$2.6\times 10^{-7}$	1.99
Y6m3a	$7$	no	Solution A in Table 1 of [20]	$1.7\times 10^{-3}$	2.42	$1.3\times 10^{-3}$	2.32
KL6s9a	$9$	no	s9odr6a in Appendix A of [21]	$2.5\times 10^{-4}$	2.26	$2.0\times 10^{-4}$	2.18
KL6s9b	$9$	no	s9odr6b in Appendix A of [21]	$2.5\times 10^{-4}$	2.25	$2.0\times 10^{-4}$	2.18
SS6s11	$11$	no	Section 4.2 of [22]	$3.4\times 10^{-5}$	1.98	$1.7\times 10^{-5}$	1.77
SS6s13	$13$	no	Section 4.2 of [22]	$1.7\times 10^{-5}$	2.08	$4.0\times 10^{-6}$	1.64
BM6M10	$10$	no	$S_{10}$ Table 2 of [35]	$5.5\times 10^{-6}$	1.33	$2.0\times 10^{-6}$	1.12
BCE6m10b	$20$	no	BM ${}_{10}^{[6]}$ Table 1 of [31]	$5.4\times 10^{-6}$	2.65	$7.8\times 10^{-7}$	1.92
PPBCM6m9	$9$	yes	$P_{9}6$ Table 5 of [24]	$2.0\times 10^{-6}$	2.02	$1.2\times 10^{-7}$	1.26
PPBCM6m5	$11$	yes	$P_{11}6$ in Table 6 of [24]	$1.6\times 10^{-6}$	1.19	$9.6\times 10^{-7}$	1.09
PPBCM6m6	$13$	yes	$P_{13}6$ in Table 6 of [24]	$4.2\times 10^{-7}$	1.13	$2.6\times 10^{-7}$	1.04
BCE6m5	$10$	yes	$\psi_{5}^{[6]}$ Table 8 of [31]	$5.1\times 10^{-3}$	4.15	$4.0\times 10^{-3}$	3.98
BCE6m6	$12$	yes	$\psi_{6}^{[6]}$ Table 8 of [31]	$6.5\times 10^{-5}$	2.41	$2.6\times 10^{-5}$	2.07
BCE6m7	$14$	yes	$\psi_{7}^{[6]}$ Table 8 of [31]	$1.3\times 10^{-5}$	2.14	$3.0\times 10^{-6}$	1.69
BCE6m8	$16$	yes	$\psi_{8}^{[6]}$ Table 8 of [31]	$3.6\times 10^{-6}$	1.98	$4.3\times 10^{-7}$	1.39
BCE6m9	$18$	yes	$\psi_{9}^{[6]}$ Table 8 of [31]	$5.7\times 10^{-6}$	2.41	$1.7\times 10^{-7}$	1.33
BCE6m10	$20$	yes	$\psi_{10}^{[6]}$ Table 8 of [31]	$2.4\times 10^{-6}$	2.32	$9.9\times 10^{-9}$	0.93
BCE6m11	$22$	yes	$\psi_{11}^{[6]}$ Table 8 of [31]	$6.1\times 10^{-6}$	2.97	$1.5\times 10^{-8}$	1.10

Table 7: The list of 6th order product formulae with their average errors. The best results are highlighted in bold. Results for product formulae that are not products of

S_{2}

are given in italics.

label	$M$	processing	reference	$\chi$	$M\chi^{1/k}$	$\zeta$	$M\zeta^{1/k}$
S8m1	$27$	no	first Suzuki product [1]	$4.8\times 10^{-2}$	18.5	$2.3\times 10^{-2}$	16.9
S8m2	$125$	no	second Suzuki product [2]	$4.8\times 10^{-9}$	11.4	$5.0\times 10^{-13}$	3.62
Y8m7d	$15$	no	Solution D in Table 2 of [20]	$9.7\times 10^{-4}$	6.30	$1.9\times 10^{-4}$	5.15
KL8s15	$15$	no	s15odr8 in [21] (our Table 1)	$5.9\times 10^{-6}$	3.33	$2.7\times 10^{-6}$	3.02
KL8s17a	$17$	no	s17odr8a in [21]	$5.9\times 10^{-7}$	2.83	$2.3\times 10^{-7}$	2.52
KL8s17b	$17$	no	s17odr8b in [21]	$5.8\times 10^{-7}$	2.83	$2.1\times 10^{-7}$	2.49
SS8s19	$19$	no	Section 4.3 of [22]	$1.6\times 10^{-7}$	2.68	$7.2\times 10^{-8}$	2.43
SS8s21	$21$	no	Section 4.3 of [22]	$2.6\times 10^{-7}$	3.16	$8.1\times 10^{-8}$	2.73
PP8s13	$13$	yes	$P_{13}8$ in Table 6 of [24]	$1.1\times 10^{-6}$	2.33	$8.6\times 10^{-7}$	2.27
PP8s19	$19$	yes	$P_{19}8$ in Table 6 of [24]	NA	NA	$2.6\times 10^{-7}$	2.85
Y8m10	$21$	no	Table 2 (our new result)	$4.9\times 10^{-8}$	2.56	$1.1\times 10^{-8}$	2.13
Y8m10b	$21$	no	Table 2 (our new result)	$5.4\times 10^{-7}$	3.45	$1.6\times 10^{-9}$	1.67
YP8m8	$17$	yes	Table 3 (our new result)	$5.4\times 10^{-8}$	2.10	$2.2\times 10^{-9}$	1.41

Table 8: The list of 8th order product formulae with their average errors. The best results are highlighted in bold.

label	$M$	processing	reference	$\chi$	$M\chi^{1/k}$	$\zeta$	$M\zeta^{1/k}$
S10m1	$81$	no	first Suzuki product [1]	$7.5\times 10^{-2}$	62.5	$8.1\times 10^{-3}$	50.0
S10m2	$625$	no	second Suzuki product [2]	$2.6\times 10^{-13}$	34.5	$5.9\times 10^{-19}$	9.39
KL10s31a	$31$	no	s31odr10a in Appendix A of [21]	$6.1\times 10^{-6}$	9.33	$5.4\times 10^{-6}$	9.21
KL10s31b	$31$	no	s31odr10b in Appendix A of [21]	$6.4\times 10^{-5}$	11.88	$4.2\times 10^{-5}$	11.3
SS10s31	$31$	no	Section 4.4 of [22]	$3.4\times 10^{-8}$	5.55	$2.7\times 10^{-8}$	5.43
SS10s33	$33$	no	Section 4.4 of [22]	$1.0\times 10^{-8}$	5.25	$8.1\times 10^{-9}$	5.12
SS10s35	$35$	no	Section 4.4 of [22]	$8.0\times 10^{-10}$	4.31	$4.3\times 10^{-11}$	3.22
Alberdi31	$31$	no	Appendix A of [36]	$1.1\times 10^{-7}$	6.23	$1.0\times 10^{-7}$	6.18
Alberdi33	$33$	no	Appendix A of [36]	$5.9\times 10^{-8}$	6.25	$5.2\times 10^{-8}$	6.17
Alberdi35	$35$	no	Appendix A of [36]	$1.1\times 10^{-8}$	5.62	$9.3\times 10^{-9}$	5.51
PP10s19	$19$	yes	$P_{19}10$ in Table 6 of [24]	NA	NA	$5.7\times 10^{-6}$	5.68
PP10s23	$23$	yes	$P_{23}10$ in Table 6 of [24]	$2.7\times 10^{-5}$	8.04	$3.2\times 10^{-8}$	4.09
Y10m15	$31$	no	Table 4	$4.5\times 10^{-7}$	7.19	$4.1\times 10^{-7}$	7.12
Y10m16	$33$	no	Table 4	$1.9\times 10^{-8}$	5.57	$7.5\times 10^{-9}$	5.08
Y10m17	$35$	no	Table 5	$1.4\times 10^{-8}$	5.75	$1.8\times 10^{-10}$	3.71
Y10m18	$37$	no	Table 4	$3.1\times 10^{-9}$	5.22	$2.2\times 10^{-9}$	5.05
Y10m18b	$37$	no	Table 5	$2.6\times 10^{-8}$	6.46	$4.2\times 10^{-10}$	4.27

Table 9: The list of 10th order product formulae with their average errors. The best results are highlighted in bold, and the second-best are underlined.

To compare among product formulae of the same order, we compare $M\chi^{1/k}$ for spectral norm error, or $M\zeta^{1/k}$ for eigenvalue error, and the best results in each table are indicated by the values in bold. First we consider the comparison of $\chi$ for spectral-norm error. For 4th order, the lowest value for $M\chi^{1/k}$ is that of the processed formulae BCE4m7 and BCE4m8, which yield very similar values. That is out of product formulae that are products of $S_{2}$ , we will discuss the others below. In the 6th order case we find that the best is PPBCM6m6. Among the 8th order, the best performing is our processed formula YP8m8. The best 10th order is given by SS10s35, but the second-best is given by our Y10m18.

For eigenvalue error, we compare the quantity $M\zeta^{1/k}$ to establish the performance of the product formulae when they have the same order. The best performing 4th order formula for eigenvalue error is BCE4m9 for product formulae based on $S_{2}$ . For 6th order, the best one is BCE6m10. Among the 8th order product formulae, our product formula YP8m8 is still the best performing. The best 10th order is still SS10s35, though the second-best is our Y10m17. Therefore, in both cases the best 10th order is SS10s35 from Ref. [22], though we have found the second-best product formulae. More importantly, our new 8th order product formula provides the best performance both in terms of spectral norm and eigenvalue error.

In Table 6 and Table 7 we also list results for some product formulae that are specific to Hamiltonians that are a sum of only two terms, and are not constructed as a product of $S_{2}$ . In the case of fourth order, these give better performance than any of the other product formulae listed. In particular, for eigenvalue error the best performance is given by the formula BM4M6 of Blanes and Moan [35], whereas for spectral-norm error, a recently proposed formula O4M5 of Ostmeyer [34] gives the best performance. In the case of 6th order, the formula BM6M10 of Blanes and Moan gives very good performance, but still not better than some of the processed product formulae.

The Suzuki product formulae give poor performance in all cases as compared to the best product formulae, and Suzuki’s second formula always outperforms the first. In the case of 6th order, Suzuki’s second formula provides better performance than the best 6th order product formula of Yoshida in terms of eigenvalue error, but not spectral norm error. Suzuki’s second product formula also outperforms the best 8th order product formula of Yoshida. A surprising result is that Suzuki’s second product formula provides much smaller eigenvalue error than spectral-norm error, particularly for higher orders. For 10th order it is smaller by nearly 6 orders of magnitude. Nevertheless, despite that extremely small error, the large number of factors in the product mean that it has significantly worse performance than most other product formulae in the list.

To consider the threshold to use a higher-order product formula over a lower-order one, we first consider the asymptotic expression using $M\zeta^{1/k}$ for eigenvalue error. For 4th order the value of this quantity is 0.58 (BCE4m9) and for 6th order is 0.93 (BCE6m10), which gives a threshold of 290. The value for the best 8th order is 1.41 (YP8m8), which gives a threshold of around 1,200 when compared with 4th order and 22,000 when compared with 6th order. The thresholds from 4th order are larger using BM4M6, which is specific to Hamiltonians that are a sum of two parts. Then the threshold from 4th to 6th is 3800, and from 4th to 8th is 7000. At 10th order the best value is 3.22 (SS10s35), which when compared with 8th order gives a threshold of $2.2\times 10^{14}$ . This is far greater than can be expected for realistic simulations.

Most of these thresholds based on asymptotic expressions are inaccurate, because the thresholds correspond to large time steps where the asymptotic expressions break down. An exception is the threshold for 8th to 10th order. That corresponds to a time interval around $0.2$ that is sufficiently small that the threshold estimate is accurate. In contrast, the calculated threshold between 6th and 8th order would correspond to a time step larger than 3, which is far too large for the scaling of the error to be accurate.

In this case we find numerically that the threshold between these 6th and 8th order formulae (BCE6m10 and YP8m8) is about $7000$ , which is less than the threshold calculated above using the asymptotic scaling. On the other hand, the threshold for the 6th order formula PPBCM6m6 to our 8th order formula YP8m8 is 940,000. That corresponds to about 590,000 time steps, with a time step size of about $1.6$ (for the 8th order formula).

It is possible to optimise our kernel for the larger time step size in order to provide better performance. This optimisation provides a kernel (see Table 10) such that the threshold from PPBCM6m6 is about 82,000. This now corresponds to 36,000 time steps with a time step size of about $2.3$ . Thus the optimisation provides another order of magnitude where 8th order product formulae are optimal, albeit with a product formula tailored to the time step size.

	8th order kernel for large time steps
$w_{1}$	$0.17292977711543507576156846186751$
$w_{2}$	$0.27170302438610082629780995303455$
$w_{3}$	$0.21909548236117584671732473611764$
$w_{4}$	$-0.37248751509173994928726577188503$
$w_{5}$	$0.12371215242829562284552662906028$
$w_{6}$	$-0.38248795584080401246916516638535$
$w_{7}$	$0.54458116939694551608378519711877$
$w_{8}$	$0.21703219494512028923053251824215$

Table 10: The 8th order kernel tailored for large time step size to improve the threshold between 6th and 8th order.

One must also be careful in considering thresholds based on eigenvalue error, because there needs to be a large number of time steps in order for the spectral-norm error in the complete evolution time to be dominated by the eigenvalue error. There is a factor of about 24 between the spectral-norm error and the eigenvalue error for the 8th order product formula YP8m8. For the threshold between 6th order and 8th order the number of time steps is orders of magnitude larger than this, which indicates that the eigenvalue error is an appropriate measure.

These threshold calculations are for Hamiltonians that are sums of two terms, each of which is normalised. In the case where the Hamiltonian is not normalised, the threshold should be scaled by the norm. In the simple case where the Hamiltonian is a sum of two terms with equal norms $\|A\|=\|B\|$ , then the threshold will be for $\|A\|T/\epsilon$ rather than $T/\epsilon$ . More generally the Hamiltonian may be a sum of any number of terms with different norms, which would change the threshold and make it unclear what quantity the threshold should be considered in terms of. There is also the possibility that the threshold can be changed for Hamiltonians that have structure to them instead of being random, for example those for quantum chemistry.

VI.2 Comparison for fermionic Hamiltonians

Fermionic Hamiltonians encountered in quantum chemistry often have the form

\displaystyle\sum_{p,q=1}^{d}\tau_{pq}a^{{\dagger}}_{p}a_{q}+\sum_{p,q=1}^{d}\nu_{pq}a^{{\dagger}}_{p}a_{p}a^{{\dagger}}_{q}a_{q}

(48)

where $a^{\dagger}_{p}$ and $a_{p}$ are the fermionic creation and annihilation operators acting on orbital $p$ , and there are a total of $d$ orbitals. Each entry $\tau_{pq}$ , $\nu_{pq}$ is real and there is symmetry in exchanging indices. The behaviour of the error as the size of system is changed can be predicted based on the result in Theorem 4 of [37].

Theorem 5 (Theorem 4 in [37]).

Let $H=T+V=\sum_{p,q}^{d}\tau_{pq}a^{{\dagger}}_{p}a_{q}+\sum_{p,q}^{d}\nu_{pq}a^{{\dagger}}_{p}a_{p}a^{{\dagger}}_{q}a_{q}$ be an interacting-electronic Hamiltonian, and $S_{k}(t)$ be a $k$ th order product formula splitting the evolutions under $T$ and $V$ . Then

\norm{S_{k}(t)-e^{-itH}}_{W_{\eta}}=\order{(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})^{k-1}\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta t^{k+1}}\,,

(49)

where $\norm{\cdot}_{W_{\eta}}$ corresponds to the operator norm on the operator acting in the $\eta$ -electron subspace, $\norm{\tau}_{1}=\max_{p}\sum_{q}\absolutevalue{\tau_{pq}}$ and $\norm{\nu}_{1,[\eta]}=\max_{p}\max_{q_{1}<q_{1}<...<q_{\eta}}(\absolutevalue{\nu_{pq_{1}}}+\cdots+\absolutevalue{\nu_{pq_{\eta}}})$ .

To fairly compare product formulae as applied to quantum chemistry, we can define $\xi,\omega$ so that

	$\displaystyle\chi$	$\displaystyle=\xi\,(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})^{k-1}\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta\,,$		(50)
	$\displaystyle\zeta$	$\displaystyle=\omega\,(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})^{k-1}\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta\,.$		(51)

That is, $\chi,\xi$ are for the spectral-norm error, and $\zeta,\omega$ are for the eigenvalue error, with $\xi,\omega$ being the quantities defined for fermionic Hamiltonians. Note that the expression in Ref. [37] was derived for error in the spectral norm, but that also provides an upper bound on the error in the eigenvalues, so it is reasonable to consider the constant defined for eigenvalue error. Then the formula for the threshold $T/\varepsilon$ becomes, for $\omega$ (eigenvalue error)

$\displaystyle\frac{T}{\varepsilon}$	$\displaystyle=\left(\frac{M_{2}[\omega_{2}(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})^{k_{2}-1}\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta]^{1/k_{2}}}{M_{1}[\omega_{1}(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})^{k_{1}-1}\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta]^{1/k_{1}}}\right)^{\frac{1}{\frac{1}{k_{1}}-\frac{1}{k_{2}}}}$
	$\displaystyle=\left(\frac{M_{2}[\omega_{2}\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta/(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})]^{1/k_{2}}}{M_{1}[\omega_{1}\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta/(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})]^{1/k_{1}}}\right)^{\frac{1}{\frac{1}{k_{1}}-\frac{1}{k_{2}}}}$
$\displaystyle\frac{\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta}{\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}}\frac{T}{\varepsilon}$	$\displaystyle=\left(\frac{M_{2}\omega_{2}^{1/k_{2}}}{M_{1}\omega_{1}^{1/k_{1}}}\right)^{\frac{1}{\frac{1}{k_{1}}-\frac{1}{k_{2}}}}.$	(52)

Thus we see that the ratio ${\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta}/(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]})$ governs the threshold where a 10th-order product formula will improve over an 8th-order product formula (and similarly for other orders). This is somewhat different from the factor of the norm of the Hamiltonian that might otherwise be expected for unstructured Hamiltonians.

To analyse the non-asymptotic regime, we rewrite the error to be proportional to

\frac{\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta}{\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}}t\times\left[{t}\left(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\right)\right]^{k}.

(53)

Following the derivation of Ref. [37], it is easily seen that the contributions to the error from each of the orders higher than $k+1$ will be of a similar form. These higher-order contributions to the error can be used to determine the error in the non-asymptotic regime where the order $k+1$ term is not dominant.

To see why that is the case, note that Lemma 2 of Ref. [37] describes the error at order $k+1$ ( $p+1$ in the notation of that work) in terms of a sum over norms of multicommutators of $k+1$ terms in the Hamiltonian. It is easily seen that the error at each of the higher orders will also be given by similar multicommutator forms. The method of proof of Theorem 4 in Ref. [37] is based on simplifying the multicommutator expressions from Lemma 2, and proceeds in exactly the same way for the higher orders.

Therefore, the sum over errors of all orders would give a total error of the form

\frac{\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta}{\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}}t\times\sum_{j=k+1}^{\infty}\xi_{j}\left[t\left(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\right)\right]^{j}.

(54)

This sum would then give the non-asymptotic form for the error

\frac{\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta}{\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}}t\times g\!\left({t}\left(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\right)\right)\,,

(55)

for some function $g$ .

Then following the same procedure to determine the threshold as before, we will have a total error given by

\frac{\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta}{\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}}T\times g\!\left({t}\left(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\right)\right)\,.

(56)

That would imply that comparing two product formulae, the threshold will be obtained for

g_{1}\!\left({t_{1}}\left(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\right)\right)=g_{2}\!\left({t_{2}}\left(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\right)\right)\,.

(57)

Using $t_{2}=t_{1}M_{2}/M_{1}$ and solving for $t_{1}$ as before, the $T/\epsilon$ threshold would then be given by substituting the solution for $t_{1}$ into

\frac{\norm{\tau}_{1}\norm{\nu}_{1,[\eta]}\eta}{\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}}\frac{T}{\epsilon}=1/g_{1}\!\left({t_{1}}\left(\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\right)\right)\,.

(58)

The threshold can be calculated for Hamiltonians with a given value of $\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}$ , then for some new Hamiltonian with different values of the norms, the threshold should still be for the same value of ${t_{1}}\big{(}\norm{\tau}_{1}+\norm{\nu}_{1,[\eta]}\big{)}$ . This implies that the right-hand-side of Eq. (58) should be the same for the threshold with the new Hamiltonian. As a result, we expect that the threshold should be in terms of the left-hand-side of Eq. (58), which is the same expression as found in the asymptotic regime.

In our numerical testing, the coefficients $\tau_{pq}$ , $\nu_{pq}$ were chosen uniformly from the interval $[-1,1]$ . We computed $\zeta$ for a selection of the product formulae with the best performance for random Hamiltonians. With the errors computed numerically, we can compute $\omega$ . In Table 11, Table 12 and Table 13 we give the computed result for $d=6$ orbitals, assuming half-filling of the orbitals. Our numerics indicate that the error is roughly proportional to the bound in Eq. 49, independent of the number of orbitals, though the constant factors are small. In Appendix C we give the constant factors for the case $d=4$ , showing that the computed constant does not change much with a different $d$ .

label	$\omega$	$M\omega^{1/k}$
BCE6m10	$3.4\times 10^{-11}$	$0.36$
PPBCM6m6	$1.2\times 10^{-9}$	$0.42$

Table 11: Comparison of constant factors

\omega

for the best product formulae for 6th order. We generate

1000

random Hamiltonians with

d=6

orbitals as in Eq. 48 and compute the average

\omega

label	$\omega$	$M\omega^{1/k}$
SS8s19	$3.5\times 10^{-11}$	$0.94$
Y8m10	$8.5\times 10^{-12}$	$0.87$
Y8m10b	$1.3\times 10^{-12}$	$0.68$
PP8s13	$4.4\times 10^{-10}$	$0.88$
YP8m8	$1.7\times 10^{-12}$	$0.57$

Table 12: Comparison of constant factors

\omega

for the best product formulae for 8th order. We generate

1000

random Hamiltonians with

d=6

orbitals as in Eq. 48 and compute the average

\omega

label	$\omega$	$M\omega^{1/k}$
SS10s35	$3.0\times 10^{-15}$	$1.24$
Y10m17	$2.6\times 10^{-14}$	1.53
Y10m18b	$5.1\times 10^{-14}$	1.74
PP10s23	$2.3\times 10^{-12}$	1.58

Table 13: Comparison of constant factors

\omega

for the best product formulae for 10th order. We generate

1000

random Hamiltonians with

d=6

orbitals as in Eq. 48 and compute the average

\omega

For calculating the thresholds based on the asymptotic formula in Eq. 52, we should compare our best 8th order formula YP8m8 to SS10s35 in the 10th order case, and BCE6m10 in the 6th order case. Then we obtain a threshold of $3\times 10^{13}$ for 8th to 10th order, and $6\times 10^{4}$ for 6th to 8th order, indicating about 9 orders of magnitude for $T/\epsilon$ where our 8th order product formula is optimal.

Similar to the case for general matrices, we expect that the threshold for 8th to 10th order is sufficiently large that the time steps are small and the asymptotic formulae are accurate, but the threshold between 6th and 8th order may be unreliable. In that case we solve Eq. 57 for $t_{1}$ with $t_{2}=t_{1}M_{2}/M_{2}$ , and find the threshold using Eq. 58 is increased to $5\times 10^{5}$ . Similar to the case with general matrices, we expect that the 6th order product formula PPBCM6m6 can give better performance for larger time steps, and indeed we find that the threshold is further increased to $8\times 10^{6}$ . By using our 8th order kernel optimized for large time step in Table 10 we can improve the threshold between PPBCM6m6 and our 8th order formula in the non-asymptotic case. In this case we obtain a threshold of around $2\times 10^{6}$ , which is an improvement by a factor of 4 times; not as significant as the improvement for general random matrices.

Next we compare the performance of the product formulae for specific systems. From Ref. [37] the norms can be expected to scale as

\norm{\nu}_{1,[\eta]}=\mathcal{O}\left(\frac{\eta^{2/3}N^{1/3}}{\Omega^{1/3}}\right),\qquad\norm{\tau}_{1}=\mathcal{O}\left(\frac{N^{2/3}}{\Omega^{2/3}}\right),

(59)

where $N$ is the number of orbitals and $\Omega$ is the volume (denoted $n$ and $\omega$ in [37]). The constants for the scaling of these norms are derived in Ref. [38] as

\norm{\nu}_{1,[\eta]}\approx\pi^{1/3}(3/4)^{2/3}\frac{\eta^{2/3}N^{1/3}}{\Omega^{1/3}},\qquad\norm{\tau}_{1}\approx\frac{3\pi^{2}N^{2/3}}{2\Omega^{2/3}}.

(60)

For an initial order of magnitude estimate, the chemical accuracy required for the phase estimation is about $0.001$ Hartree, which implies $T$ of about $1000\pi$ , and $\varepsilon$ can be taken to be of order 1. With high estimates $\eta\approx 100$ and $N/\Omega\approx 10^{9}$ , the left-hand-side (LHS) of Eq. (52) would be on the order of $10^{10}$ , which means our 8th order should be preferred to 10th order. Therefore, in this case either 8th order or 10th order could be chosen depending on the details of the system

For comparison, we also consider the Alpha + Hydrogen system from Ref. [38]. The parameters $\eta$ , $N$ and $\Omega$ are given in Table II of that work. Some of the most challenging values of $T,\epsilon$ are given in Figure 8 of that work, with for example $T=40$ and $\epsilon=10^{-4}$ . With these parameters the LHS of Eq. (52) is approximately $8.6\times 10^{9}$ , which is below the threshold (between 8th and 10th order) by a few orders of magnitude. The smallest values of $T/\epsilon$ from Ref. [38] are obtained with $T=10$ and $\epsilon=0.01$ (see Table IV). That results in a value for the LHS of Eq. (52) of $2.2\times 10^{7}$ , which is above the threshold between 6th and 8th order, so our 8th order product formula would still be optimal. This indicates that the new 8th order product formula found here should be optimal for the full range of parameter values used in Ref. [38]. Nonetheless, for smaller values of $T/\epsilon$ it may be advantageous to further customise the product formula to provide better performance at larger time steps, as described above to improve the threshold between 6th order and 8th order.

Note also that the fact that the threshold here is different than for random Hamiltonians means that the threshold will depend on the class of Hamiltonians. Selecting a different type of Hamiltonians can change the threshold. Similarly, it will be expected that choosing a particular class of fermionic Hamiltonians (such as those arising from plane waves) rather than random fermionic Hamiltonians can give a different threshold. Nevertheless, changing the class of Hamiltonians still resulted in an extremely large threshold for 10th order product formulae to be optimal. This suggests that the 8th order formula will be optimal for realistic simulations regardless of the class of Hamiltonians.

VII Conclusion

In this work we have developed a method of fairly comparing product formulae with different numbers of factors and different orders, as well as searched for improved higher-order product formulae. We find that it is better to compare the eigenvalue error, rather than the spectral-norm error as in prior work, because it is the eigenvalue error in a single step which dominates the error over longer time evolutions. The optimal order of integrator to use depends primarily on the ratio $T/\epsilon$ , the ratio of the total evolution time to the required error.

This is because, as $T/\epsilon$ becomes larger, the error allowed for each individual time step becomes smaller, and the time steps must become shorter to make the error sufficiently small. That results in a larger overall number of time steps, and so a larger complexity. The higher-order integrators better reduce the error, with shorter time steps, and so for large $T/\epsilon$ that improvement more than compensates for the larger number of exponentials for the higher-order product formulae.

In our search for higher-order product formulae we found very large numbers of solutions at both 8th and 10th order. Our solution at 8th order improves over all prior 8th order product formulae that we found in prior work. The best 8th order product formula we have found is a processed product formula, which means that in a long time simulation the kernel would be repeated many times, and the processor would only be used at the beginning and end of the simulation. For applications in eigenvalue estimation the processor would not be needed at all. We found highly accurate solutions at 10th order, though they are outperformed by one result found in prior work that is extremely accurate. These product formulae greatly outperform the fractal product formulae of Suzuki, which are commonly considered in quantum simulation but require a large number of exponentials.

We show via numerical testing for random Hamiltonians that there is a range of about eight orders of magnitude for $T/\epsilon$ where our new 8th-order product formula outperforms all other known product formulae. Moreover, our 8th order product formula can be further adjusted for smaller values of $T/\epsilon$ to provide another order of magnitude where it is optimal. This range includes reasonable combinations of parameters that one would use in applications. A lower-order formula would only be optimal for small simulations that are likely to be classically tractable. A 10th order formula would only be optimal for simulations well beyond the scale considered for quantum algorithms.

A particularly important class of Hamiltonians is those corresponding to fermions (rather than just random), because those are relevant to simulations in quantum chemistry. Using the formula from Ref. [37], those give a threshold for $T/\epsilon$ combined with norms of the matrices describing the fermionic system. As an example we consider parameters for the simulations from Ref. [38], and estimate that the threshold for 10th order product formulae to outperform our 8th order product formula is about three orders of magnitude larger than the largest simulations considered there. Those simulations are already expected to require more than $10^{13}$ Toffolis, so even larger simulations where 10th order product formulae are optimal would be well beyond the scale that could be realistically implemented on foreseeable quantum computers.

A potential topic for future work is further customising product formulae for large time steps. We have customised the 8th order product formula to improve the threshold, but it would also be possible to customise the 6th and 4th order product formulae. It would also be possible to customise the product formulae for larger time steps with the fermionic Hamiltonians. So far we have only customised the 8th order product formula for general matrices.

VIII Acknowledgements

MESM was supported by the ARC Centre of Excellence for Quantum Computation and Communication Technology (CQC2T), project number CE170100012 and a scholarship top-up and extension from the Sydney Quantum Academy. MESM and YRS were supported by the Defense Advanced Research Projects Agency under Contract No. HR001122C0074. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the Defense Advanced Research Projects Agency. DKB acknowledges funding by the Australian Research Council (project numbers FT190100106, DP210101367, CE170100009). DWB worked on this project under a sponsored research agreement with Google Quantum AI. DWB is also supported by Australian Research Council Discovery Projects DP210101367 and DP220101602.

References

Suzuki [1990] M. Suzuki, Physics Letters A 146, 319 (1990).
Suzuki [1991] M. Suzuki, Journal of Mathematical Physics 32, 400 (1991).
Berry et al. [2007] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Communications in Mathematical Physics 270, 359 (2007).
Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
Aharonov and Ta-Shma [2003] D. Aharonov and A. Ta-Shma, in Proceedings of the Thirty-Fifth Annual ACM Symposium on Theory of Computing, STOC ’03 (Association for Computing Machinery, New York, NY, USA, 2003) p. 20–29.
Berry and Childs [2012] D. W. Berry and A. M. Childs, Quantum Information and Computation 12, 29–62 (2012).
Berry et al. [2014] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, in Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing, STOC ’14 (Association for Computing Machinery, New York, NY, USA, 2014) p. 283–292.
Berry et al. [2015a] D. W. Berry, A. M. Childs, and R. Kothari, in 2015 IEEE 56th Annual Symposium on Foundations of Computer Science (2015) pp. 792–809.
Low [2019] G. H. Low, in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) p. 491–502.
Childs et al. [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Physical Review X 11, 011020 (2021).
Childs and Wiebe [2012] A. M. Childs and N. Wiebe, Quantum Information and Computation 12, 901–924 (2012).
Berry et al. [2015b] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Physical Review Letters 114, 90502 (2015b).
Low and Chuang [2017] G. H. Low and I. L. Chuang, Physical Review Letters 118, 010501 (2017).
Babbush et al. [2015] R. Babbush, J. McClean, D. Wecker, A. Aspuru-Guzik, and N. Wiebe, Physical Review A 91, 022311 (2015).
Childs et al. [2018] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences 115, 9456 (2018).
Su et al. [2021] Y. Su, H.-Y. Huang, and E. T. Campbell, Quantum 5, 495 (2021).
Zhang [2012] C. Zhang, in Monte Carlo and Quasi-Monte Carlo Methods 2010, edited by L. Plaskota and H. Woźniakowski (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012) pp. 709–719.
Campbell [2019] E. Campbell, Physical Review Letters 123, 070503 (2019).
Childs et al. [2019] A. M. Childs, A. Ostrander, and Y. Su, Quantum 3, 182 (2019).
Yoshida [1990] H. Yoshida, Physics Letters A 150, 262 (1990).
Kahan and Li [1997] W. Kahan and R.-C. Li, Mathematics of Computation 66, 1089–1099 (1997).
Sofroniou and Spaletta [2005] M. Sofroniou and G. Spaletta, Optimization Methods and Software 20, 597 (2005).
Blanes et al. [2008] S. Blanes, F. Casas, and A. Murua, arXiv:0812.0377 (2008).
Blanes et al. [2006] S. Blanes, F. Casas, and A. Murua, SIAM Journal on Scientific Computing 27, 1817 (2006).
Blanes [2001] S. Blanes, Applied Numerical Nathematics 37, 289 (2001).
Blanes et al. [1999] S. Blanes, F. Casas, and J. Ros, SIAM Journal on Scientific Computing 21, 711 (1999).
Butcher [1969] J. C. Butcher, in Conference on the Numerical Solution of Differential Equations, edited by J. L. Morris (Springer Berlin Heidelberg, Berlin, Heidelberg, 1969) pp. 133–139.
Butcher and Sanz-Serna [1996] J. Butcher and J. Sanz-Serna, Applied Numerical Mathematics 22, 103 (1996).
McLachlan [1996] R. I. McLachlan, Integration Algorithms and Classical Mechanics 10, 141 (1996).
Wisdom et al. [1996] J. Wisdom, M. Holman, and J. Touma, Fields Institute Communications 10, 217 (1996).
Blanes et al. [2024] S. Blanes, F. Casas, and A. Escorihuela-Tomàs, arXiv: 2404.04340 (2024).
Blanes and Casas [2004] S. Blanes and F. Casas, Linear Algebra and its Applications 378, 135 (2004).
Blanes et al. [2013] S. Blanes, F. Casas, P. Chartier, and A. Murua, Mathematics of Computation 82, 1559 (2013).
Ostmeyer [2023] J. Ostmeyer, Journal of Physics A: Mathematical and Theoretical 56, 285303 (2023).
Blanes and Moan [2002] S. Blanes and P. Moan, Journal of Computational and Applied Mathematics 142, 313 (2002).
Alberdi et al. [2019] E. Alberdi, M. Antoñana, J. Makazaga, and A. Murua, arXiv: 1909.07263 (2019).
Low et al. [2023] G. H. Low, Y. Su, Y. Tong, and M. C. Tran, PRX Quantum 4, 020323 (2023).
Rubin et al. [2024] N. C. Rubin, D. W. Berry, A. Kononov, F. D. Malone, T. Khattar, A. White, J. Lee, H. Neven, R. Babbush, and A. D. Baczewski, Proceedings of the National Academy of Sciences 121, e2317772121 (2024).
Van-Brunt and Visser [2016] A. Van-Brunt and M. Visser, Journal of Mathematical Physics 57, 023507 (2016).
Duleba and Karcz-Duleba [2020] I. Duleba and I. Karcz-Duleba, in Computer Aided Systems Theory – EUROCAST 2019, edited by R. Moreno-Díaz, F. Pichler, and A. Quesada-Arencibia (Springer International Publishing, Cham, 2020) pp. 465–473.

Appendix A Extending Yoshida’s method to 10th order

Here we explain how to extend the method of Yoshida to obtain the equations for a 10th order integrator. See Section II.2 for an introductory explanation of how the method is used for 6th order. The general principle used there was to provide an expression for $S^{(m)}(\tau)$ in Eq. (19) with the expression in the exponential given up to 6th order. Then the coefficients of the multicommutators of operators needed to be made equal to zero in order to obtain a 6th order approximation. Here we apply the same principle, except now we need to derive the terms up to 10th order.

In order to do this, we start by expressing $S_{2}(t)$ up to 10th order as

S_{2}(t)=e^{t\alpha_{1}+t^{3}\alpha_{3}+t^{5}\alpha_{5}+t^{7}\alpha_{7}+t^{9}\alpha_{9}+\mathcal{O}(t^{11})}

(61)

where we follow the notation from Corollary 3 with $\alpha_{j}$ defined as commutators of operators. We will then iteratively apply the symmetric BCH expansion for $Z$ in $e^{Z}=e^{C}e^{D}e^{C}$ in order to obtain an expression for $S^{(m)}(\tau)$ up to 10th order. It will be helpful to consider the following notation for commutators. For any operators $X_{1},X_{2},\cdots,X_{L}$ we define

[X_{1}^{n_{1}},X_{2}^{n_{1}},\cdots,X_{L}^{n_{L}}]:=[\underbrace{X_{1},\cdots,X_{1}}_{\text{{\hbox{n_1}} times }},\underbrace{X_{2},\cdots,X_{2}}_{\text{{\hbox{n_2}} times }},\cdots,\underbrace{X_{L},\cdots,X_{L}}_{\text{{\hbox{n_L}} times }}],

(62)

where the commutator on the right hand should be understood as in Eq. 9. To express the results it will be useful to define the following commutators

$\displaystyle\beta_{9}$	$\displaystyle=[\alpha_{1},\alpha_{1},\alpha_{7}],$	(63)
$\displaystyle\gamma^{(1)}_{9}$	$\displaystyle=[\alpha_{1},\alpha_{3},\alpha_{5}],$	(64)
$\displaystyle\gamma^{(2)}_{9}$	$\displaystyle=[\alpha_{3},\alpha_{1},\alpha_{5}],$	(65)
$\displaystyle\gamma^{(3)}_{9}$	$\displaystyle=[\alpha_{5},\alpha_{1},\alpha_{3}],$	(66)
$\displaystyle\delta^{(1)}_{9}$	$\displaystyle=[\alpha_{1}^{4},\alpha_{5}],$	(67)
$\displaystyle\delta^{(2)}_{9}$	$\displaystyle=[\alpha_{3},\alpha_{1}^{3},\alpha_{3}],$	(68)
$\displaystyle\delta^{(3)}_{9}$	$\displaystyle=[\alpha_{1},\alpha_{3},\alpha_{1}^{2},\alpha_{3}],$	(69)
$\displaystyle\epsilon_{9}$	$\displaystyle=[\alpha_{1}^{6},\alpha_{3}].$	(70)

We also use the notation of Yoshida [20] for the commutators

$\displaystyle\beta_{5}$	$\displaystyle=[\alpha_{1},\alpha_{1},\alpha_{3}]\equiv E_{5,2},$	(71)
$\displaystyle\beta_{7}$	$\displaystyle=[\alpha_{1},\alpha_{1},\alpha_{5}]\equiv E_{7,3},$	(72)
$\displaystyle\gamma_{7}$	$\displaystyle=[\alpha_{3},\alpha_{3},\alpha_{1}]\equiv-E_{7,2},$	(73)
$\displaystyle\delta_{7}$	$\displaystyle=[\alpha_{1},\alpha_{1},\alpha_{1},\alpha_{1},\alpha_{3}]\equiv E_{7,4},$	(74)

where we have indicated the equivalence to the multicommutators in the notation of Ref. [24].

Note that here we have only defined commutators of up to 7 of the $\alpha_{j}$ operators. In the following we will be expanding expressions up to 9th order in $t$ , because the symmetry of the formulae means that 10th order terms (and all even-order terms) must be zero. The only way to obtain 9th order in $t$ with commutators of $t\alpha_{1}$ , $t^{3}\alpha_{3}$ , etc., is to have commutators of all $\alpha_{1}$ as $[\alpha_{1},\alpha_{1},\cdots,\alpha_{1}]$ . But that expression must be zero, because $[\alpha_{1},\alpha_{1}]=0$ . Then when we express the expansion of $Z$ in $e^{Z}=e^{C}e^{D}e^{C}$ , the first-order terms in $C$ and $D$ will be proportional to $\alpha_{1}$ . That means the only 9th order terms coming from commutators of 9 of $C,D$ would correspond to $[\alpha_{1},\alpha_{1},\cdots,\alpha_{1}]$ and therefore be zero. This means we will only need commutators of up to 7 operators when expressing the expansion of $Z$ in $e^{Z}=e^{C}e^{D}e^{C}$ as well.

To obtain the coefficients multiplying the commutators with up to 7 operators in the symmetric BCH expansion for $e^{Z}=e^{C}e^{D}e^{C}$ , we first use the algorithm defined in Section V of Ref. [39]. The algorithm in that work generates the scalar coefficients multiplying products of operators rather than their commutators, as we need here. In order to derive the corresponding expressions in terms of commutators, we express the symmetric BCH expansion in the Ph. Hall basis, which is a basis for writing Lie monomials consisting of commutators of the generators of the Lie algebra. For example, the elements up to 4th order in $C,D$ are the operators $C,D$ themselves, as well as

[C,D],\quad[C,C,D],\quad[D,C,D],\quad[C,C,C,D],\quad[D,C,C,D],\quad[D,D,C,D].

(75)

For a list of operators in this basis up to 7th order, see Table 1 in [40].

We obtain the coefficients for the Ph. Hall basis by solving the corresponding linear problem of changing from one basis to another. As an example, consider the term with $3$ operators in the symmetric BCH expansion from Corollary 3, given as $\alpha_{3}=\frac{1}{12}[Y,[Y,X]]-\frac{1}{24}[X,[X,Y]]$ . We can also express the commutators as products by expanding out the commutators, which gives $\tilde{\alpha}_{3}=\frac{1}{24}\left(2Y^{2}X-4YXY-2XY^{2}-X^{2}Y+2XYX+YX^{2}\right)$ . The algorithm in [39] outputs expressions with the commutators expanded out as in $\tilde{\alpha}_{3}$ . In order to obtain the original expression $\alpha_{3}$ , we write $\tilde{\alpha}_{3}=a[Y,[Y,X]]+b[X,[X,Y]]$ with $a,b\in\mathbb{R}$ and expand the commutators $[Y,[Y,X]]$ and $[X,[X,Y]]$ . This gives several linear equations that can be written in terms of a matrix. By inverting this matrix, we obtain the coefficients $a$ and $b$ .

By using this method, we obtain the symmetric BCH expansion for $e^{Z}=e^{C}e^{D}e^{C}$ up to 7th order as

$\displaystyle Z$	$\displaystyle=2C+D+\frac{1}{6}([D,D,C]-[C,C,D])$
	$\displaystyle\quad+\frac{7}{360}[C,C,C,C,D]-\frac{1}{360}[D,D,D,D,C]$
	$\displaystyle\quad+\frac{1}{90}[C,D,D,D,C]+\frac{1}{45}[D,C,C,C,D]$
	$\displaystyle\quad-\frac{1}{60}[C,C,D,D,C]+\frac{1}{30}[D,D,C,C,D]$
	$\displaystyle\quad-\frac{31}{15120}[C,C,C,C,C,C,D]-\frac{31}{5040}[D,C,C,C,C,C,D]$
	$\displaystyle\quad-\frac{13}{1890}[D,D,C,C,C,C,D]-\frac{53}{15120}[D,D,D,C,C,C,D]$
	$\displaystyle\quad-\frac{1}{1260}[D,D,D,D,C,C,D]-\frac{1}{15120}[D,D,D,D,D,C,D]+\mathcal{R}_{(9\leq)}.$	(76)

Here $\mathcal{R}_{(9\leq)}$ is an infinite sum with commutators of an odd number of operators equal to or higher than $9$ . Now we prove the following Lemma which gives us the expansion up to 10th order of $S^{(m)}(\tau)$ . This will allow us to derive the equations for 10th order product formulae.

Lemma 6.

Using the definition for $S^{(m)}(\tau)$ as given in Eq. 15, we have that for all $m\in\mathbb{N}$

$\displaystyle S^{(m)}(\tau)$	$\displaystyle=\exp\{\tau A_{1,m}\alpha_{1}+\tau^{3}A_{3,m}\alpha_{3}+\tau^{5}(A_{5,m}\alpha_{5}+B_{5,m}\beta_{5})$
	$\displaystyle\quad+\tau^{7}(A_{7,m}\alpha_{7}+B_{7,m}\beta_{7}+C_{7,m}\gamma_{7}+D_{7,m}\delta_{7})$
	$\displaystyle\quad+\tau^{9}(A_{9,m}\alpha_{9}+B_{9,m}\beta_{9}+C^{(1)}_{9,m}\gamma^{(1)}_{9}+C^{(2)}_{9,m}\gamma^{(2)}_{9}+C^{(3)}_{9,m}\gamma^{(3)}_{9}$
	$\displaystyle\quad+D^{(1)}_{9,m}\delta^{(1)}_{9}+D^{(2)}_{9,m}\delta^{(2)}_{9}+D^{(3)}_{9,m}\delta^{(3)}_{9}+E_{9,m}\epsilon_{9})+\order{\tau^{11}}\bigg{\}},$	(77)

where the variables in upper case denote polynomials in the variables $(w_{1},\cdots,w_{m})$ .

Proof.

We proceed by induction. First, note that the statement is true for the case $m=0$ ,

	$\displaystyle S^{(m=0)}(\tau)$	$\displaystyle=S_{2}(w_{0}\tau),$		(78)
		$\displaystyle=\exp{tw_{0}\alpha_{1}+t^{3}w_{0}^{3}\alpha_{3}+t^{5}w_{0}^{5}\alpha_{5}+t^{7}w_{0}^{7}\alpha_{7}+t^{9}w_{0}^{9}\alpha_{9}+\mathcal{O}(t^{11})}.$		(79)

This clearly has the form of Lemma 6 by taking $A_{j,m=0}=w_{0}^{j}$ and all other scalar variables as $0$ .

Assume now that Lemma 6 is correct, we want to derive an expression for $S^{(m+1)}$ . We then have $S^{(m+1)}(\tau)=S_{2}(w_{m+1}\tau)S^{(m)}(\tau)S_{2}(w_{m+1}\tau)$ and thus

$\displaystyle S_{2}(w_{m+1}\tau)S^{(m)}(\tau)S_{2}(w_{m+1}\tau)$	$\displaystyle=\exp\{\tau w_{m+1}\alpha_{1}+\tau^{3}w_{m+1}^{3}\alpha_{3}+\tau^{5}w_{m+1}^{5}\alpha_{5}+\tau^{7}w_{m+1}^{7}\alpha_{7}+\tau^{9}w_{m+1}^{9}\alpha_{9}+\mathcal{O}(\tau^{11})\bigg{\}}$
	$\displaystyle\quad\times\exp\{\tau A_{1,m}\alpha_{1}+\tau^{3}A_{3,m}\alpha_{3}+\tau^{5}(A_{5,m}\alpha_{5}+B_{5,m}\beta_{5})$
	$\displaystyle\quad+\tau^{7}(A_{7,m}\alpha_{7}+B_{7,m}\beta_{7}+C_{7,m}\gamma_{7}+D_{7,m}\delta_{7})$
	$\displaystyle\quad+\tau^{9}(A_{9,m}\alpha_{9}+B_{9,m}\beta_{9}+C^{(1)}_{9,m}\gamma^{(1)}_{9}+C^{(2)}_{9,m}\gamma^{(2)}_{9}+C^{(3)}_{9,m}\gamma^{(3)}_{9}$
	$\displaystyle\quad+D^{(1)}_{9,m}\delta^{(1)}_{9}+D^{(2)}_{9,m}\delta^{(2)}_{9}+D^{(3)}_{9,m}\delta^{(3)}_{9}+E_{9,m}\epsilon_{9})+\order{\tau^{11}}\bigg{\}}$
	$\displaystyle\quad\times\exp\{\tau w_{m+1}\alpha_{1}+\tau^{3}w_{m+1}^{3}\alpha_{3}+\tau^{5}w_{m+1}^{5}\alpha_{5}+\tau^{7}w_{m+1}^{7}\alpha_{7}+\tau^{9}w_{m+1}^{9}\alpha_{9}+\order{\tau^{11}}\bigg{\}}.$	(80)

We compute the right-hand-side (RHS) of Appendix A applying the symmetric BCH formula from Corollary 3. Writing the RHS as $e^{C}e^{D}e^{C}$ , we have that

$\displaystyle C$	$\displaystyle=\tau w_{m+1}\alpha_{1}+\tau^{3}w_{m+1}^{3}\alpha_{3}+\tau^{5}w_{m+1}^{5}\alpha_{5}+\tau^{7}w_{m+1}^{7}\alpha_{7}+\tau^{9}w_{m+1}^{9}\alpha_{9}+\mathcal{O}(\tau^{11})$	(81)
$\displaystyle D$	$\displaystyle=\tau A_{1,m}\alpha_{1}+\tau^{3}A_{3,m}\alpha_{3}+\tau^{5}(A_{5,m}\alpha_{5}+B_{5,m}\beta_{5})+\tau^{7}(A_{7,m}\alpha_{7}+B_{7,m}\beta_{7}+C_{7,m}\gamma_{7}+D_{7,m}\delta_{7})$
	$\displaystyle\quad+\tau^{9}(A_{9,m}\alpha_{9}+B_{9,m}\beta_{9}+C^{(1)}_{9,m}\gamma^{(1)}_{9}+C^{(2)}_{9,m}\gamma^{(2)}_{9}+C^{(3)}_{9,m}\gamma^{(3)}_{9}+D^{(1)}_{9,m}\delta^{(1)}_{9}+D^{(2)}_{9,m}\delta^{(2)}_{9}+D^{(3)}_{9,m}\delta^{(3)}_{9}+E_{9,m}\epsilon_{9})$
	$\displaystyle\quad+\order{\tau^{11}}.$	(82)

We then compute the commutators of $C$ and $D$ that appear in the symmetric BCH formula, here we give the resulting 9th order operators after applying the commutators. When we write $[C,D,\cdots,C]_{9}$ , the subscript indicates that we are only keeping the $9$ th order terms when expanding the commutator. We will explain in detail how to compute the commutator $\mathcal{C}=[D,D,C]_{9}$ , the other commutators are computed in a similar way. Since we only need to consider terms of 9th order when computing $\mathcal{C}$ , each term will have contributions from each operator inside $\mathcal{C}$ (in this case two operators $D$ and one $C$ ) which is comprised of odd numbers that sum up to $9$ such that the commutator is non-zero. We then have that

$\displaystyle[D,D,C]_{9}$	$\displaystyle=\tau^{9}(A^{2}_{1,m}w_{m+1}^{7}-A_{1,m}A_{7,m}w_{m+1})[\alpha_{1},\alpha_{1},\alpha_{7}]$
	$\displaystyle\quad+\tau^{9}A_{1,m}B_{7,m}w_{m+1}[\alpha_{1},\beta_{7},\alpha_{1}]$
	$\displaystyle\quad+\tau^{9}A_{1,m}C_{7,m}w_{m+1}[\alpha_{1},\gamma_{7},\alpha_{1}]$
	$\displaystyle\quad+\tau^{9}A_{1,m}D_{7,m}w_{m+1}[\alpha_{1},\delta_{7},\alpha_{1}]$
	$\displaystyle\quad+\tau^{9}(A_{1,m}A_{3,m}w_{m+1}-A_{1,m}A_{5,m}w_{m+1}^{3})[\alpha_{1},\alpha_{3},\alpha_{5}]$
	$\displaystyle\quad+\tau^{9}A_{1,m}B_{5,m}w_{m+1}^{3}[\alpha_{1},\beta_{5},\alpha_{3}]$
	$\displaystyle\quad+\tau^{9}(A_{3,m}A_{1,m}w_{m+1}^{5}-A_{3,m}A_{5,m}w_{m+1})[\alpha_{3},\alpha_{1},\alpha_{5}]$
	$\displaystyle\quad+\tau^{9}A_{1,m}B_{5,m}w_{m+1}[\alpha_{3},\beta_{5},\alpha_{1}]$
	$\displaystyle\quad+\tau^{9}(A_{5,m}A_{1,m}w_{m+1}^{3}-A_{5,m}A_{3,m}w_{m+1})[\alpha_{5},\alpha_{1},\alpha_{3}]$
	$\displaystyle\quad+\tau^{9}B_{5,m}A_{3,m}w_{m+1}[\beta_{5},\alpha_{3},\alpha_{1}].$	(83)

Given how we have defined the commutator, we have then

$\displaystyle[D,D,C]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}^{2}w_{m+1}^{7}-A_{1,m}A_{7,m}w_{m+1})\beta_{9}-\tau^{9}A_{1,m}B_{7,m}w_{m+1}\delta^{(1)}_{9}$
	$\displaystyle\quad+\tau^{9}A_{1,m}C_{7,m}w_{m+1}\delta^{(3)}_{9}-\tau^{9}A_{1,m}D_{7,m}w_{m+1}\epsilon_{9}$
	$\displaystyle\quad+\tau^{9}(A_{1,m}A_{3,m}w_{m+1}^{5}-A_{1,m}A_{5,m}w_{1}^{3})\gamma^{(1)}_{9}-\tau^{9}A_{1,m}B_{5,m}w_{m+1}^{3}\delta^{(3)}_{9}$
	$\displaystyle\quad+\tau^{9}(A_{3,m}A_{1,m}w_{m+1}^{5}-A_{3,m}A_{5,m}w_{m+1})\gamma^{(2)}_{9}-\tau^{9}A_{3,m}B_{5,m}w_{m+1}\delta^{(2)}_{9}$
	$\displaystyle\quad+\tau^{9}(A_{5,m}A_{1,m}w_{m+1}^{3}-A_{5,m}A_{3,m}w_{m+1}^{2})\gamma^{(3)}_{9}$
	$\displaystyle\quad+\tau^{9}(B_{5,m}A_{1,m}w_{m+1}^{3}-B_{5,m}A_{3,m}w_{m+1})(\delta^{(2)}_{9}-\delta^{(3)}_{9})$	(84)
$\displaystyle[C,C,D]_{9}$	$\displaystyle=\tau^{9}(w_{m+1}^{2}A_{7,m}-w_{m+1}^{8}A_{1,m})\beta_{9}+\tau^{9}w_{m+1}^{2}B_{7,m}\delta^{(1)}_{9}$
	$\displaystyle\quad-\tau^{9}w_{m+1}^{2}C_{7,m}\delta^{(3)}_{9}+\tau^{9}w_{m+1}^{2}D_{7,m}\epsilon_{9}$
	$\displaystyle\quad+\tau^{9}(w_{m+1}^{4}A_{5,m}-w_{m+1}^{6}A_{3,m})\gamma^{(1)}_{9}+\tau^{9}w_{m+1}^{4}B_{5,m}\delta^{(3)}_{9}$
	$\displaystyle\quad+\tau^{9}(w_{m+1}^{4}A_{5,m}-w_{m+1}^{8}A_{1,m})\gamma^{(2)}_{9}+\tau^{9}w_{m+1}^{4}B_{5,m}\delta^{(2)}_{9}$
	$\displaystyle\quad+\tau^{9}(w_{m+1}^{6}A_{3,m}-w_{m+1}^{8}A_{1,m})\gamma^{(3)}_{9}$	(85)
$\displaystyle[C,C,C,C,D]_{9}$	$\displaystyle=\tau^{9}(w_{m+1}^{4}A_{5,m}-w_{m+1}^{8}A_{1,m})\delta^{(1)}_{9}+\tau^{9}A_{1,m}^{3}B_{5,m}w_{m+1}\epsilon_{9}$
	$\displaystyle\quad+\tau^{9}(w_{m+1}^{6}A_{3,m}-w_{m+1}^{8}A_{1,m})\delta^{(2)}_{9}$
	$\displaystyle\quad+\tau^{9}2(w_{m+1}^{6}A_{3,m}-w_{m+1}^{8}A_{1,m})\delta^{(3)}_{9}$	(86)
$\displaystyle[D,D,D,D,C]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}^{4}w_{m+1}^{5}-A_{1,m}^{3}A_{5,m}w_{m+1})\delta^{(1)}_{9}-\tau^{9}A_{1,m}^{3}B_{5,m}w_{m+1}\epsilon_{9}$
	$\displaystyle\quad+\tau^{9}(A_{3,m}A_{1,m}^{3}w_{m+1}^{3}-A_{3,m}^{2}A_{1,m}^{2}w_{m+1})\delta^{(2)}_{9}$
	$\displaystyle\quad+\tau^{9}2(A_{1,m}^{3}A_{3,m}w_{m+1}^{3}-A_{1,m}^{2}A_{3,m}^{2}w_{m+1})\delta^{(3)}_{9}$	(87)
$\displaystyle[C,D,D,D,C]_{9}$	$\displaystyle=\tau^{9}(w_{m+1}^{6}A_{1,m}^{3}-w_{m+1}^{2}A_{1,m}^{2}A_{5,m})\delta^{(1)}_{9}-\tau^{9}A_{1,m}^{2}B_{5,m}w_{m+1}^{2}\epsilon_{9}$
	$\displaystyle\quad+\tau^{9}(w_{m+1}^{6}A_{1,m}^{3}-w_{m+1}^{4}A_{1,m}^{2}A_{3,m})\delta^{(2)}_{9}$
	$\displaystyle\quad+\tau^{9}2(w_{m+1}^{4}A_{3,m}A_{1,m}^{2}-w_{m+1}^{2}A_{3,m}^{2}A_{1,m})\delta^{(3)}_{9}$	(88)
$\displaystyle[D,C,C,C,D]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}A_{5,m}w_{m+1}^{3}-A_{1,m}^{2}w_{m+1}^{7})\delta^{(1)}_{9}+\tau^{9}w_{m+1}^{3}A_{1,m}B_{5,m}\epsilon_{9}$
	$\displaystyle\quad+\tau^{9}(A_{3,m}^{2}w_{m+1}^{3}-A_{1,m}A_{3,m}w_{m+1}^{5})\delta^{(2)}$
	$\displaystyle\quad+\tau^{9}2(A_{1,m}A_{3,m}w_{m+1}^{5}-A_{1,m}^{2}w_{m+1}^{7})\delta^{(3)}_{9}$	(89)
$\displaystyle[C,C,D,D,C]_{9}$	$\displaystyle=\tau^{9}(w_{m+1}^{7}A_{1,m}^{2}-w_{m+1}^{3}A_{1,m}A_{5,m})\delta^{(1)}_{9}-w_{m+1}^{3}A_{1,m}B_{5,m}\epsilon_{9}$
	$\displaystyle\quad+\tau^{9}(w_{m+1}^{7}A_{1,m}^{2}-w_{m+1}^{5}A_{1,m}A_{3,m})\delta^{(2)}_{9}$
	$\displaystyle\quad+\tau^{9}(w_{m+1}^{5}A_{1,m}A_{3,m}-w_{m+1}^{3}A_{3,m}^{2})\delta^{(3)}_{9}$	(90)
$\displaystyle[D,D,C,C,D]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}^{2}A_{5,m}w_{m+1}^{2}-A_{1,m}^{3}w_{m+1}^{6})\delta^{(1)}_{9}+A_{1,m}^{2}B_{5,m}w_{m+1}^{2}\epsilon_{9}$
	$\displaystyle\quad+(A_{3,m}^{2}A_{1,m}w_{m+1}^{2}-A_{3,m}A_{1,m}^{2}w_{m+1}^{4})\delta^{(2)}_{9}$
	$\displaystyle\quad+\tau^{9}(A_{3,m}^{2}A_{1,m}w_{m+1}^{2}-A_{3,m}A_{1,m}^{2}w_{m+1}^{4})\delta^{(3)}_{9}+\tau^{9}(A_{1,m}^{2}A_{3,m}w_{m+1}^{4}-A_{1,m}^{3}w_{m+1}^{6})\delta^{(3)}_{9}$	(91)
$\displaystyle[C,C,C,C,C,C,D]_{9}$	$\displaystyle=\tau^{9}(w_{m+1}^{6}A_{3,m}-w_{m+1}^{8}A_{1,m})\epsilon_{9}$
$\displaystyle[D,C,C,C,C,C,D]_{9}$	$\displaystyle=\tau^{9}(w_{m+1}^{5}A_{1,m}A_{3,m}-w_{m+1}^{7}A_{1,m}^{2})\epsilon_{9}$	(92)
$\displaystyle[D,D,C,C,C,C,D]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}^{2}A_{3,m}w_{m+1}^{4}-A_{1,m}^{3}w_{m+1}^{6})\epsilon_{9}$	(93)
$\displaystyle[D,D,D,C,C,C,D]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}^{3}A_{3,m}w_{m+1}^{3}-A_{1,m}^{4}w_{m+1}^{5})\epsilon_{9}$	(94)
$\displaystyle[D,D,D,D,C,C,D]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}^{4}A_{3,m}w_{m+1}^{2}-A_{1,m}^{5}w_{m+1}^{4})\epsilon_{9}$	(95)
$\displaystyle[D,D,D,D,D,C,D]_{9}$	$\displaystyle=\tau^{9}(A_{1,m}^{5}A_{3,m}w_{m+1}-A_{1,m}^{6}w_{m+1}^{3})\epsilon_{9}$	(96)

Note that all the terms previously computed have terms that can be written as in Lemma 6, thus proving that $S^{(m+1)}$ can also be written in this way. ∎

Having proved Lemma 6, we can now compute the polynomials in Lemma 6. The polynomials are obtained from the recursion in Appendix A, the left hand side corresponds to $S^{(m+1)}$ and can can be written as a single exponential, the same is true of the right side which is written as a single exponential. We have then the following polynomials:

$\displaystyle A_{9,m+1}$	$\displaystyle=A_{9,m}+2w_{m+1}^{9}$	(97)
$\displaystyle B_{9,m+1}$	$\displaystyle=B_{9,m}+\frac{1}{6}(A_{1,m}^{2}w_{m+1}^{7}-A_{1,m}A_{7,m}w_{m+1})$
	$\displaystyle\quad-\frac{1}{6}(A_{7,m}w_{m+1}^{2}-A_{1,m}w_{m+1}^{8})$	(98)
$\displaystyle C_{9,m+1}^{(1)}$	$\displaystyle=C_{9,m}^{(1)}+\frac{1}{6}(A_{3,m}^{2}A_{1,m}w_{m+1}^{5}-A_{1,m}A_{5,m}w_{m+1}^{3})$
	$\displaystyle\quad-\frac{1}{6}(A_{5,m}w_{m+1}^{4}-A_{3,m}w_{m+1}^{6})$	(99)
$\displaystyle C_{9,m+1}^{(2)}$	$\displaystyle=C_{9,m}^{(2)}+\frac{1}{6}(A_{3,m}^{2}A_{1,m}w_{m+1}^{5}-A_{3,m}A_{5,m}w_{m+1})$
	$\displaystyle\quad-\frac{1}{6}(A_{5,m}w_{m+1}^{4}-A_{1,m}w_{m+1}^{8})$	(100)
$\displaystyle C_{9,m+1}^{(3)}$	$\displaystyle=C_{9,m}^{(3)}+\frac{1}{6}(A_{5,m}A_{1,m}w_{m+1}^{3}-A_{3,m}A_{5,m}w_{m+1})$
	$\displaystyle\quad-\frac{1}{6}(A_{3,m}w_{m+1}^{6}-A_{1,m}w_{m+1}^{8})$	(101)
$\displaystyle D_{9,m+1}^{(1)}$	$\displaystyle=D_{9,m}^{(1)}-\frac{1}{6}(A_{1,m}B_{7,m}w_{m+1}+w_{m+1}^{2}B_{7,m})$
	$\displaystyle\quad+\frac{7}{360}(A_{5,m}w_{m+1}^{4}-w_{m+1}^{8}A_{1,m})$
	$\displaystyle\quad-\frac{1}{360}(A_{1,m}^{4}w_{m+1}^{5}-A_{1,m}^{3}A_{5,m}w_{m+1})$
	$\displaystyle\quad+\frac{1}{90}(A_{1,m}^{3}w_{m+1}^{6}-A_{1,m}^{2}A_{5,m}w_{m+1}^{2})$
	$\displaystyle\quad+\frac{1}{45}(A_{1,m}A_{5,m}w_{m+1}^{3}-A_{1,m}^{2}w_{m+1}^{7})$
	$\displaystyle\quad-\frac{1}{60}(A_{1,m}^{2}w_{m+1}^{7}-A_{1,m}A_{5,m}w_{m+1}^{3})$
	$\displaystyle\quad+\frac{1}{30}(A_{1,m}^{2}A_{5,m}w_{m+1}^{2}-A_{1,m}^{3}w_{m+1}^{6})$	(102)
$\displaystyle D_{9,m+1}^{(2)}$	$\displaystyle=D_{9,m}^{(2)}-\frac{1}{6}(A_{3,m}B_{5,m}w_{m+1}+w_{m+1}^{4}B_{5,m})$
	$\displaystyle\quad+\frac{7}{360}(A_{3,m}w_{m+1}^{6}-w_{m+1}^{8}A_{1,m})$
	$\displaystyle\quad-\frac{1}{360}(A_{1,m}^{3}A_{3,m}w_{m+1}^{3}-A_{1,m}^{2}A_{3,m}^{2}w_{m+1})$
	$\displaystyle\quad+\frac{1}{90}(A_{1,m}^{3}w_{m+1}^{6}-A_{1,m}^{2}A_{3,m}w_{m+1}^{4})$
	$\displaystyle\quad+\frac{1}{45}(A_{3,m}^{2}w_{m+1}^{3}-A_{1,m}A_{3,m}w_{m+1}^{5})$
	$\displaystyle\quad-\frac{1}{60}(A_{1,m}^{2}w_{m+1}^{7}-A_{1,m}A_{3,m}w_{m+1}^{5})$
	$\displaystyle\quad+\frac{1}{30}(A_{3,m}^{2}A_{1,m}w_{m+1}^{2}-A_{1,m}^{2}A_{3,m}w_{m+1}^{4})$	(103)
$\displaystyle D_{9,m+1}^{(3)}$	$\displaystyle=D_{9,m}^{(3)}-\frac{1}{6}(A_{1,m}B_{5,m}w_{m+1}^{3}+w_{m+1}^{4}B_{5,m})$
	$\displaystyle\quad+\frac{14}{360}(A_{3,m}w_{m+1}^{6}-w_{m+1}^{8}A_{1,m})$
	$\displaystyle\quad-\frac{2}{360}(A_{1,m}^{3}A_{3,m}w_{m+1}^{3}-A_{1,m}^{2}A_{3,m}^{2}w_{m+1})$
	$\displaystyle\quad+\frac{2}{90}(A_{1,m}^{2}A_{3,m}w_{m+1}^{4}-A_{1,m}A_{3,m}^{2}w_{m+1}^{2})$
	$\displaystyle\quad+\frac{2}{45}(A_{3,m}A_{1,m}w_{m+1}^{5}-A_{1,m}^{2}w_{m+1}^{7})$
	$\displaystyle\quad-\frac{1}{60}(A_{1,m}^{2}w_{m+1}^{7}-A_{1,m}A_{3,m}w_{m+1}^{5})$
	$\displaystyle\quad+\frac{1}{30}(A_{3,m}^{2}A_{1,m}w_{m+1}^{2}-A_{1,m}^{2}A_{3,m}w_{m+1}^{4})$
	$\displaystyle\quad+\frac{1}{6}(A_{1,m}C_{7,m}w_{m+1}+w_{m+1}^{2}c_{7,m})$
	$\displaystyle\quad-\frac{1}{60}(w_{m+1}^{5}A_{1,m}A_{3,m}-w_{m+1}^{3}A_{3,m}^{2})$
	$\displaystyle\quad+\frac{1}{30}(A_{1,m}^{2}A_{3,m}w_{m+1}^{4}-A_{1,m}^{3}w_{m+1}^{6})$
	$\displaystyle\quad-\frac{1}{6}(B_{5,m}A_{1,m}w_{m+1}^{3}-B_{5,m}A_{3,m}w_{m+1})$	(104)
$\displaystyle E_{9,m+1}$	$\displaystyle=E_{9,m}-\frac{1}{6}(A_{1,m}D_{7,m}-w_{m+1}^{2}D_{7,m})$
	$\displaystyle\quad+\frac{7}{360}w_{m+1}^{4}B_{5,m}$
	$\displaystyle\quad+\frac{1}{360}A_{1,m}^{3}B_{5,m}w_{m+1}$
	$\displaystyle\quad-\frac{1}{90}A_{1,m}^{2}B_{5,m}w_{m+1}^{2}$
	$\displaystyle\quad+\frac{1}{45}A_{1,m}B_{5,m}w_{m+1}^{3}$
	$\displaystyle\quad+\frac{1}{60}A_{1,m}B_{5,m}w_{m+1}^{3}$
	$\displaystyle\quad+\frac{1}{30}A_{1,m}^{2}B_{5,m}w_{m+1}^{2}$
	$\displaystyle\quad-\frac{31}{15120}(w_{m+1}^{6}A_{3,m}-w_{m+1}^{8}A_{1,m})$
	$\displaystyle\quad-\frac{31}{5040}(w_{m+1}^{5}A_{1,m}A_{3,m}-w_{m+1}^{7}A_{1,m}^{2})$
	$\displaystyle\quad-\frac{13}{1890}(A_{1,m}^{2}A_{3,m}w_{m+1}^{4}-A_{1,m}^{3}w_{m+1}^{6})$
	$\displaystyle\quad-\frac{53}{15120}(A_{1,m}^{3}A_{3,m}w_{m+1}^{3}-A_{1,m}^{4}w_{m+1}^{5})$
	$\displaystyle\quad-\frac{1}{1260}(A_{1,m}^{4}A_{3,m}w_{m+1}^{2}-A_{1,m}^{5}w_{m+1}^{4})$
	$\displaystyle\quad-\frac{1}{15120}(A_{1,m}^{5}A_{3,m}w_{m+1}-A_{1,m}^{6}w_{m+1}^{3}).$	(105)

We obtain the polynomial equations for the tenth order product formula by imposing that $A_{1,m}=1$ and all other terms are equal to zero. Because $C^{(2)}_{9,m}=C^{(1)}_{9,m}+C^{(3)}_{9,m}$ , one equation is eliminated, there are 15 equations to solve.

Appendix B Method for determining processors

A processed formula is composed of two elements: a kernel, $\Sigma(t)$ , and a processor, $P(t)$ . The effective order $k$ captures up to which order in $t$ the full product formula, including the processor, reproduces the the target dynamics, $P(t)\Sigma(t)P(t)^{-1}=e^{(X+Y)t}+\order{t^{k+1}}$ . In this work, we use processors that are constructed with the same procedure as Ref. [24]. This type of processors are products of $S_{2}(\tau w_{i})$ arranged as $Q^{(m)}(\tau)Q^{(m)}(-\tau)$ , where $Q^{(m)}(\tau)=S_{2}(\tau w_{m})S_{2}(\tau w_{m-1})\cdots S_{2}(\tau w_{0})$ . In the language of Ref. [24], $Q^{(m)}(\pm\tau)$ is an element of the group of integrators $\mathcal{G}_{3}$ . The same reference gives a basis for the generating algebra, which we used to perform the following calculations. The basis can be found in Table 2 of Ref. [24], and we use the same naming scheme for the algebra elements.

We can obtain recursive formulas for $Q^{(m)}(\tau)=S_{2}(\tau w_{m})Q^{(m-1)}(\tau)$ via successive applications of the BCH formula (as opposed to the symmetric BCH formula in the case of the kernel). After taking the product $Q^{(m)}(\tau)Q^{(m)}(-\tau)$ , the processor simplifies compared to $Q^{(m)}$ . In particular, the logarithm of $P(\tau)$ is zero up to terms of order $t^{3}$ . We illustrate the procedure that gives an iterative expression for the processor up to order 8, omitting the calculations. First, we write the product $S_{2}(\tau w_{1})S_{2}(\tau w_{0})$ in the $E_{i,j}$ -basis. We have, for the logarithm of $Q^{(1)}(\tau)=S_{2}(w_{1}\tau)S_{2}(w_{0}\tau)$ ,

\displaystyle\log(Q^{(1)}(\tau))

\displaystyle=p_{1,1}^{(1)}Y_{1}t+p_{3,1}^{(1)}Y_{3}t^{3}+p_{4,1}^{(1)}E_{4,1}t^{4}+\left(p_{5,1}^{(1)}E_{5,1}t^{5}+p_{5,2}^{(1)}E_{5,2}\right)+\left(p_{6,1}^{(1)}E_{6,1}+p_{6,2}^{(1)}E_{6,2}\right)t^{6}+\order{t^{7}}\,,

(106)

where each coefficient is a function of $w_{0},w_{1}$ . Using an inductive argument analogous to the derivation in the case of the symmetric BCH formula, this expression will give a recursive expression for $Q^{(m)}(\tau)=S_{2}(\tau w_{m})Q^{(m-1)}(\tau)$ , whose logarithm is an expansion in terms of $p^{(m)}_{i,j}(w_{0},\dots,w_{m})$ . The resulting expression is identical to the one above, with the superscripts $(1)$ replaced by $(m)$ , and each coefficient is a polynomial of $w_{0},\dots,w_{m}$ . In fact, it is not necessary to compute all the $p$ -coefficients to find $P^{(m)}(\tau)$ . For example, only three of them (out of seven) are necessary to find the processor up to order $6$ , since

\displaystyle\log(P^{(m)}(\tau))=2p^{(m)}_{4,1}E_{4,1}t^{4}+p^{(m)}_{1,1}p^{(m)}_{4,1}E_{5,3}t^{5}+\frac{1}{3}\left(\left(p^{(m)}_{1,1}\right)^{2}p^{(m)}_{4,1}E_{6,1}+6p^{(m)}_{6,2}E_{6,2}\right)t^{6}+\order{t^{7}}\,.

(107)

Similar cancellations occur at any order. In our simulations, terms up to eighth order are required. The seventh-order term in the expansion is (omitting the dependency on $m$ )

\displaystyle\frac{1}{2}p_{3,1}p_{4,1}E_{7,2}+p_{1,1}p_{6,1}E_{7,3}+\left(\frac{1}{12}p_{1,1}^{3}p_{4,1}+p_{1,1}p_{6,2}\right)E_{7,4}\,,

(108)

while the eighth-order term is

	$\displaystyle\frac{1}{3}E_{8,4}p_{1,1}^{2}p_{6,1}+2E_{8,1}p_{8,1}+2E_{8,2}p_{8,2}+\frac{1}{2}E_{8,3}\left(p_{1,1}p_{3,1}p_{4,1}-p_{3,1}p_{5,2}+4p_{8,3}\right)$
	$\displaystyle+2E_{8,4}p_{8,4}+E_{8,5}\left(\frac{1}{60}P_{1,1}^{4}p_{4,1}+\frac{1}{3}p_{1,1}^{2}p_{6,2}+2p_{8,5}\right)\,.$		(109)

Note that not all the coefficients that appear in the expansion of $Q^{(m)}(\tau)$ are necessary.

Appendix C Error constants for fermionic Hamiltonians with 4 orbitals

In this appendix we provide the average error constant $\omega$ for the eigenvalue error for the evolution of fermionic Hamiltonians in the case that the number of orbitals is $d=4$ . We find that the $\omega$ obtained are close to those we obtain when $d=6$ .

label	$\omega$	$M\omega^{1/k}$
BCE6m10	$3.3\times 10^{-11}$	$0.36$
PPBCM6m6	$1.4\times 10^{-9}$	$0.43$

Table 14: Comparison of constant factors

\omega

for the best product formulae for 6th order. We generate

1000

random Hamiltonians with

d=4

orbitals as in Eq. 48 and compute the average

\omega

label	$\omega$	$M\omega^{1/k}$
SS8s19	$3.4\times 10^{-11}$	$0.93$
Y8m10	$8.7\times 10^{-12}$	$0.87$
Y8m10b	$1.5\times 10^{-12}$	$0.68$
PP8s13	$4.2\times 10^{-10}$	$0.87$
YP8m8	$2.3\times 10^{-12}$	$0.59$

Table 15: Comparison of constant factors

\omega

for the best product formulae for 8th order. We generate

1000

random Hamiltonians with

d=4

orbitals as in Eq. 48 and compute the average

\omega

label	$\omega$	$M\omega^{1/k}$
SS10s35	$3.0\times 10^{-15}$	$1.24$
Y10m17	$2.5\times 10^{-14}$	1.53
Y10m18b	$3.8\times 10^{-14}$	1.69
PP10s23	$1.5\times 10^{-12}$	1.51

Table 16: Comparison of constant factors

\omega

for the best product formulae for 10th order. We generate

1000

random Hamiltonians with

d=4

orbitals as in Eq. 48 and compute the average

\omega

$\displaystyle\\|\widetilde{U}-U\\|$	$\displaystyle=\\|e^{i\tilde{H}}De^{-i\tilde{H}}-U\\|$
	$\displaystyle=\\|(I+i\tilde{H}+\mathcal{O}(\tilde{H}^{2}))D(I-i\tilde{H}+\mathcal{O}(\tilde{H}^{2}))-U\\|$
	$\displaystyle=\\|i[\tilde{H},D]+D-U+\mathcal{O}(\tilde{H}^{2})\\|$
	$\displaystyle\leq\\|[\tilde{H},D-I]\\|+\\|D-U\\|+\mathcal{O}(\tilde{H}^{2})\,.$	(43)

$\displaystyle\\|\widetilde{U}^{r}-U^{r}\\|$	$\displaystyle=\\|VD^{r}V^{\dagger}-U^{r}\\|$
	$\displaystyle=\\|(VD^{r}V^{\dagger}-D^{r}V^{\dagger})+(D^{r}V^{\dagger}-U^{r}V^{\dagger})+(U^{r}V^{\dagger}-U^{r})\\|$
	$\displaystyle\leq\\|V-I\\|+\\|D^{r}-U^{r}\\|+\\|V^{\dagger}-I\\|$
	$\displaystyle\leq 2\\|V-I\\|+r\\|D-U\\|\,.$	(44)