Novel Trotter formulas for digital quantum simulation

TL;DR

This paper introduces new Trotter formulas for digital quantum simulation that optimize fidelity over different time scales, enhancing the efficiency and robustness of quantum simulations on noisy devices.

Contribution

It presents a geometric framework for optimizing Trotter expansion order, leading to two novel second-order formulas with improved fidelity for quantum simulation.

Findings

Two new second-order Trotter expansions with optimized fidelity.

Framework for operation order optimization based on geometric principles.

Enhanced robustness and efficiency in digital quantum simulation.

Abstract

Quantum simulation promises to address many challenges in fields ranging from quantum chemistry to material science, and high-energy physics, and could be implemented in noisy intermediate-scale quantum devices. A challenge in building good digital quantum simulators is the fidelity of the engineered dynamics given a finite set of elementary operations. Here we present a framework for optimizing the order of operations based on a geometric picture, thus abstracting from the operation details and achieving computational efficiency. Based on this geometric framework, we provide two alternative second-order Trotter expansions, one with optimal fidelity at a short time scale, and the second robust at a long time scale. Thanks to the improved fidelity at different time scale, the two expansions we introduce can form the basis for experimental-constrained digital quantum simulation.

Figures

Equations

$$e^{p{a}_{1}A/n}{e}^{q{b}_{1}B/n}\dots e^{p{a}_{M}A/n}{e}^{q{b}_{M}B/n}\approx e^{pA+qB},$$

$$F_U=e^{-(pA+qB)}(e^{p{a}_{1}A/n}{e}^{q{b}_{1}B/n}\dots e^{p{a}_{M}A/n}{e}^{q{b}_{M}B/n})=$$

$$\text{openone}+{\mathcal{E}}_2[A,B]+{\mathcal{E}}_{3,A}[A,[A,B]]+{\mathcal{E}}_{3,B}[B,[A,B]]+\mathcal{O}(t^4)$$

$$(e^{pA/(2n)}{e}^{qB/n}{e}^{pA/(2n)}{)}^n,\text{with}$$

$$F_U=\text{openone}-\frac{1}{24n^2}pq(p[A,[A,B]]+2q[B,[A,B])+\mathcal{O}(t^4)$$

$$\frac{1}{2}[A,B](\int x_{1}d{x}_2-x_{2}d{x}_1)=[A,B]\iint d{x}_{1}d{x}_2,$$

$$\frac{1}{3}[A,[A,B]]\int x_1(x_{1}d{x}_2-x_{2}d{x}_1)$$

$$+\frac{1}{3}[B,[A,B]]\int x_2(x_{1}d{x}_2-x_{2}d{x}_1)$$

$$=[A,[A,B]]\iint x_{1}d{x}_{1}d{x}_2+[B,[A,B]]\iint x_{2}d{x}_{1}d{x}_2.$$

$$F$$

$$\omega =Ad{x}_1+Bd{x}_2,$$

$$\mathcal{P}\exp ({\int }_l\omega ),$$

$$e^{-pA-qB}{\left(e^{p{a}_{1}A/n}{e}^{q{b}_{1}B/n}\dots e^{p{a}_{M}A/n}{e}^{q{b}_{M}B/n}\right)}^n$$

$$=\mathcal{P}\exp ({\oint }_l\omega ),$$

$$\mathcal{P}\exp ({\oint }_l\omega )\approx \exp (\iint \Omega ),$$

$$\frac{1}{2}{\Omega }_{ij}{\int }_l{x}_{i}d{x}_j+\frac{1}{3}({\partial }_k{\Omega }_{ij}+[{\omega }_k,{\Omega }_{ij}]){\int }_l{x}_i{x}_{k}d{x}_j+O(l^4).$$

$$\Omega =d\omega +\omega \wedge \omega =[A,B]d{x}_{1}d{x}_2,$$

$$\frac{1}{2}[A,B](\int x_{1}d{x}_2-x_{2}d{x}_1)+\frac{1}{3}[A,[A,B]]\int x_1(x_{1}d{x}_2-x_{2}d{x}_1)+\frac{1}{3}[B,[A,B]]\int x_2(x_{1}d{x}_2-x_{2}d{x}_1)$$

$$=$$

$$e^{\frac{1}{2}xpA}{e}^{xqB}{e}^{\frac{1}{2}xpA}$$

$$e^{\frac{1}{2}xqB}{e}^{xpA}{e}^{\frac{1}{2}xqB}$$

$$\mathcal{E}=-\frac{pq}{24}(pC+2qD),$$

$${\mathcal{E}}^{\prime }=-\frac{pq}{24}(2pC+qD),$$

$$e^{\frac{a_{M+1}}{2}xpA}{e}^{\frac{b_M}{2}xqB}\cdots e^{\frac{a_2}{2}xpA}(e^{\frac{b_1}{2}xqB}{e}^{a_{1}xpA}{e}^{\frac{b_1}{2}xqB})$$

$$e^{\frac{a_2}{2}xpA}\cdots e^{\frac{b_M}{2}xqB}{e}^{\frac{a_{M+1}}{2}xpA},$$

$$e^{\frac{1}{2}{b}_{1}xqB}{e}^{a_{1}xpA}{e}^{\frac{1}{2}{b}_{1}xqB}$$

$${\mathcal{E}}_1$$

$$e^{\frac{1}{2}{a}_{2}xpA}(e^{a_{1}xpA+b_{1}xqB}+x^3{\mathcal{E}}_1)e^{\frac{1}{2}{a}_{2}xpA}$$

$$=e^{(a_1+a_2)xpA+b_{1}xqB}+x^3({\mathcal{E}}_2+{\mathcal{E}}_1)+O(x^4),$$

$${\mathcal{E}}_2$$

$${\mathcal{E}}_{2k-1}$$

$${\mathcal{E}}_{2k}$$

$$\mathcal{C}=$$

$$\mathcal{D}=$$

$${\mathcal{E}}_{\text{2D}}=\frac{pq}{24}(p\mathcal{C}C+q\mathcal{D}D)$$

Full text

Simulation has been at the core of quantum information processing right from its inception, starting from Feynman’s vision of simulating physics using a quantum system Feynman (1982). Quantum simulators are posed to be one of the first quantum devices to show task-specific quantum supremacy Preskill (2012). Quantum simulation has great potential impact on quantum chemistry O’Malley et al. (2016a); Aspuru-Guzik et al. (2005); Hempel et al. (2018), material science Babbush et al. (2018), condensed matter Lanyon et al. (2011); Salathé et al. (2015); Langford et al. (2017); Wei et al. (2018), and high-energy physics Martinez et al. (2016); Zohar et al. (2016). The most flexible strategy to achieve quantum simulation is via digital quantum simulation Lloyd (1996), where a target time-evolution operator is represented by a sequence of elementary quantum gates, usually involving one or two qubits. The strategy of approximating the continuous evolution with discrete gates is also known as Trotter expansion Trotter (1959). First- and second-order Trotter expansions with simple alternating patterns are most commonly used, although there has been a recent interest also in randomized Trotter expansions Andrew M. Childs (2018); Campbell (2019).

Finding a sequence of elementary operations that improves the fidelity to a desired simulated Hamiltonian is crucial to practical implementations of digital quantum simulation. Experimental platforms include superconducting qubits Lamata et al. (2018); Langford et al. (2017); Salathé et al. (2015); O’Malley et al. (2016b), trapped ions Arrazola et al. (2016); Lanyon et al. (2011), atomic systems Weimer et al. (2010), and spin systems Cai and Barthel (2013); Wei et al. (2018); Ajoy and Cappellaro (2013).

While the fidelity of digital quantum simulation can be increased by increasing the number of Trotter steps, on all experimental platforms, the smallest Trotter step will eventually hit some practical limitations. In many cases, it is desirable to engineer the simulated evolution from a finite set of elementary operations (gates). For example, in fault-tolerant quantum computation, the circuit is built up from a set of universal gates that can be implemented fault tolerantly Nielsen and Chuang (2000); Jones et al. (2012). In Hamiltonian simulation, the rotation angle under a specific Hamiltonian can not be made arbitrarily small due to experimental constraints. For example, in a recent implementation of digital quantum simulation with trapped ions Lanyon et al. (2011), the smallest flip angle of the unitary evolution block was $Jt=\pi/16$, and similar numbers can be obtained in other experimental platforms Barends et al. (2015). Then, even the second-order Trotter sequence might not yield good enough fidelity, while higher-order expansions are usually hard to implement experimentally due to the complicated, sometimes negative Suzuki (1991) and even imaginary Janke and Sauer (1992), coefficients.

In this paper, we focus on finding second-order Trotter expansions that yield a better fidelity than the conventional second-order Trotter (2T) expansion Suzuki (1986, 1990), given the set of implementable gates. Finding the optimal sequence that yields the best fidelity requires integer optimization over the product of large operators describing the exact form of each gate, which can be computationally expensive. Here we present a geometric framework to optimize the sequence without calculating the quantum mechanical propagators. Using this picture, we present two expansions: one minimizes the third-order Trotter error, and we call it 2-Optimal (2O); the second one has better performance at a longer time scale, and we name it 2-Diagonal (2D). Both methods have better fidelity than the 2T expansion and even 3rd-order Trotter when the time step is larger.

Geometric method to calculate the Trotter error –

We consider the problem of engineering a target operator $e^{-it\sum_{L}w_{L}H_{L}}$, which encapsulates a general quantum simulation task, using implementable gates $e^{-iH_{k}t/n}$. For illustration purpose let’s first consider $L=2$ and the target is $e^{pA+qB}$, where $A=-iH_{1}t$ and $B=-iH_{2}t$. Here we assume $p$ and $q$ to be integers and that the implementable gates are $e^{A/n}$ and $e^{B/n}$ (where $n>0$ is also an integer). Under the constraints mentioned above, we can generally write any expansion as

[TABLE]

where $\sum_{k\!=\!1}^{M}\!a_{k}\!=\!\!\sum_{k=1}^{M}\!b_{k}\!=\!n$ and $pa_{k}$, $qb_{k}$ are integers $\forall k$. The unitary fidelity $F_{U}$ of expansion (1) up to $\mathcal{O}(t^{3})$ is

[TABLE]

where $\mathcal{E}_{2}$ is of order $t^{2}$, and $\mathcal{E}_{3,A},~{}\mathcal{E}_{3,B}$ are of order $t^{3}$. For example, the first-order Trotter expansion $(e^{pA/n}e^{qB/n})^{n}$ gives $F_{U}=\openone+\frac{1}{2n}pq[A,B]+\mathcal{O}(t^{3})$. If $p$ (or $q$) is even, we can engineer the 2T expansion,

[TABLE]

to achieve a better approximation (a similar expression holds for $q$ even).

In general, finding the coefficients $\{pa_{k},qb_{k}\}$ is a hard problem since it requires an integer optimization search performed by calculating $F_{U}$ based on the exact form of the Hamiltonians $H_{1,2}$, which might be prohibitive when considering large quantum systems. Here we present a geometric method to approach this algebraic problem, which enables us to optimize the ordering without considering the exact forms of the Hamiltonians.

For any expansion in Eq. (1), we can accurately calculate the Trotter error for all possible paths up to at least the 3rd-order using a simple grid picture, as shown in Fig. 1(a). The ideal operator $e^{pA+qB}$ can be represented on the grid as a vector $(np,nq)$. One elementary gate is represented by a unit vector: $e^{A/n}$ is the unit vector (1,0) and $e^{B/n}$ is the unit vector (0,1). Starting from the origin $(0,0)$, at each step we can move right (evolving under $e^{A/n}$) or up (evolving under $e^{B/n}$). A given ordering is then represented as a directed path on the grid. Fig. 1(b) shows three different orderings for $p=4,q=3$ and $n=1$. This grid representation can be regarded as a “projection” of the Trotterized process on the unitary group. More rigorously, the Trotterized unitary operators live in a fiber bundle, for which the grid is the base space (see Appendix. A for details). On this grid picture, any path that ends at $(np,nq)$ is accurate to first order. The second-order error of the path is given by (see Appendix. A for derivation)

[TABLE]

where the coefficient of $[A,B]$ can be interpreted as the total area defined by the path and the diagonal: $\iint dx_{1}dx_{2}$, where the areas below (above) the diagonal have density $+1$ ($-1$). Then, to cancel the second-order error, the area should sum up to zero. The third-order error is (see Appendix. A for derivation)

[TABLE]

In the third-order error, the coefficient of $[A,[A,B]]$($[B,[A,B]]$) is the total moment about the $x_{2}$-axis($x_{1}$-axis) defined by the path [see Fig. 1(a)]. If all steps are restricted to be directional, that is, either along $(0,1)$ or $(1,0)$, the third-order error is always non-zero Suzuki (1991). In principle, a third-order expansion, which has zero total moments about both axes, can also be represented by this picture, by allowing walking left (-1,0) ($e^{-A/n}$) and down (0,-1) ($e^{-B/n}$). This indeed leads, e.g., to the known 3rd-order Trotter expansion known as Ruth’s formula Ruth (1983). Given experimental constraints that typically do not allow simple inversion of the Hamiltonian arrow of time, we do not consider this scenario.

Using this geometric framework, we can now discuss two second-order sequences optimizing over different cost functions. The first sequence (2O) has a global minimum third-order error. The other sequence (2D) stays as close to the diagonal as possible and minimizes the distance from each node on the path to the diagonal. 2D has a better performance at a longer time scale, even if it is not optimal in terms of the third-order error. The algorithms used to find the two sequences have different computational complexities.

2-Optimal Sequence –

In order to minimize the third-order error, we start by assigning each edge on the grid a triplet of weights given by Eq. (4) and the first line of Eq. (Geometric method to calculate the Trotter error –). Using dynamic programming, we can find a path that puts the second-order error to zero and minimizes the third-order error. The main idea is that at each step, we keep track of the accumulated weights for each node, and finally select the path corresponding to the smallest absolute value weight. The details of the algorithm can be found in Appendix. B. In general, the optimization should be done over all Trotter steps, thus on a grid $pn\times qn$. When the expansions become deep and optimizing over all Trotter steps is too expensive to calculate, one can optimize within a coprime $(p,q)$ and repeat, or mirror symmetrize the smallest unit 111When both $p$ and $q$ are odd, it is not possible to cancel the second-order error and one needs to perform the optimization on a $(2p,2q)$ unit grid..

2-Diagonal Sequence –

This algorithm is not based on the exact expression of the Trotter error, but rather on the intuition that being as close to the diagonal as possible should give a good approximation. There are different metrics to quantify “being close”; 2D algorithm minimizes the distance from each node to the diagonal, thus minimizing the total distance, which can be shown to minimize the total unsigned area enclosed by the diagonal and the path as well. What is interesting is that the global distance optimal path can be found in a computationally efficient way by a greedy optimization (see Appendix C for the proof): at each step, we choose the move that ends closer to the diagonal. It can be shown that when $p,~{}q$ are mutually prime, and when at least one of them is even, 2D always finds a unique path, and the 2D path guarantees the second-order error to be zero. When both $p$ and $q$ are odd, 2D finds two opposite-order paths with opposite second-order error. We can then symmetrize the sequence to cancel the second-order error. That is, 2D can be made into a second-order expansion by optimizing and symmetrizing on the smallest unit grid ($2(p+q)$ steps), and the whole sequence is generated by repeating the smallest unit.

For some $(p,q)$s, 2D and 2O share the same ordering, but in general, 2D does not have a global minimum third-order error. Still, we can prove that 2D still has smaller 3rd-order error compared with the 2T expansion. The details of the proof can be found in Appendix D.

Both 2O and 2D can be easily generalized to higher dimension (see Fig. 4(b) in Appendix. E), that is, to the scenario where we want to combine a larger number of propagators, $e^{A/n},e^{B/n},e^{C/n}$, etc. Fig. 1(b) shows 2O, 2D, and 2T paths for $p=4,~{}q=3$ and $n=1$ on the grid picture. Fig. 2(a) illustrates paths for $p=12,~{}q=8$. Especially when higher dimensions are involved, the simpler greedy optimization of 2D is beneficial to reduce the computational complexity.

Numerical evaluations –

To demonstrate and compare the advantage of 2D and 2O sequences, we numerically evaluate their fidelity at short and long time scale. Though the dynamic programming time and space complexities are polynomial in the grid dimension, the computation for 2O quickly becomes both time-consuming, and space-consuming as $p$ and $q$ become larger. Here we only compare expansions up to $p+q=20$.

The metric we use to quantify the performance is the average fidelity of a quantum gate Nielsen (2002), which reads

[TABLE]

where $U_{1}$ is the ideal operator and $U_{2}$ is the approximated one. To further highlight differences in fidelity, we plot the log-fidelity $F_{l}=-\log_{10}(1-F)$. We also considered $\|U_{1}-U_{2}\|$, where $\|\cdot\|$ is Frobenius norm, as a metric of error, and the results are qualitatively the same.

The orderings of the three second-order Trotter expansions for $p\!=\!12,~{}q\!=\!8$ are shown in Fig. 2(a). Fig. 2(b) shows numerical results for a two-spin transverse-field Ising Hamiltonian, $H=12H_{1}+8H_{2}$, where $H_{1}=\frac{1}{2}(\sigma_{z}^{1}+\sigma_{z}^{2}$) and $H_{2}=\sigma_{x}^{1}\sigma_{x}^{2}$. Then $A=-iH_{1}t$, $B=-iH_{2}t$. For different $(p,q)$ pairs (here with $p+q=20$) the results look qualitatively the same. We also verified that the same qualitative results apply for different Hamiltonians, such as single-particle operators, larger transverse-field Ising models with $N=2-10$ spins, and high-dimension random Hermitian matrices.

We can first verify that all three expansions are second-order expansions with respect to $t$ by fitting the data to $F_{l}=-a\log_{10}t+b$ and extract the slope $a$ at small $t$ (see the solid lines in Fig. 2(b)). From the Taylor expansion we expect $a=6$ and the fitted values for 2O, 2D, 2T are 6.07, 5.99, and 5.99 respectively. In the short $t$ regime, 2O shows the best fidelity, as expected since it minimizes the third-order error. 2D also provides a higher fidelity than the conventional second-order Trotter (2T) expansion for all $t$. At some $t$, 2D even starts to outperform 2O, the crossover point being around $t\|A\|/n=0.14$. Fig. 2(c) shows a histogram of crossover point positions over 1,000 pairs of random hermitian matrices $H_{1,2}$. The center of the Gaussian depends highly on specific $(p,q)$, but it is always on the order of 0.1. For typical experimental values, such as $Jt=\pi/16$ Lanyon et al. (2011), and $U\Delta t/2=5/16$ Barends et al. (2015), 2D is then the best sequence.

Conversely, these results imply that 2D does not need as small a Trotter step to reach the same fidelity. Then, the number of elementary gates needed to reach a given fidelity is smaller, a practical advantage for digital simulation Campbell (2019), especially when fault-tolerant gates are required. We further evaluated the 2D performance with respect not only to the time step but also to the number of “switching” operations needed, finding that 2D performs well even under this metric (see Appendix. F-G). We also compared 2D and 2T assuming that the time-step resolution in 2T is not limited, and found that the performance of 2D is still comparable to 2T even though this condition is in favor of 2T. We further find that 2D can outperform the 3rd-order Trotter expansion and the simplest sequence achieved by alternating the two operators (see Appendix. H).

Conclusions and outlook –

Finding the exact global optimal Trotter sequence or even just the sequence minimizing the minimal third-order error can be expensive to calculate. In this letter, we introduced two Trotter sequences to achieve higher fidelity in digital quantum simulation, by exploiting a novel geometric framework to estimate the Trotterization error. By optimization the fidelity on a simple grid picture, our 2O and 2D sequences outperform the widely used 2T sequence in complementary Trotter-step regimes. In addition, the 2D solution is based on an intuitive ordering that can be found with an efficient, greedy optimization algorithm. Compared to 2T, less Trotter steps are required to reach the same fidelity with our sequences, thus also providing an advantage when the Trotter steps number is a performance metric.

We note that while we mostly presented numerical results for two qubits, we numerically obtain similar performances for larger systems, and indeed the scalings found are expected to apply to many-qubit systems. As the minimization algorithm is very efficient, since it does not require to evaluate quantum propagators, and it can be easily extended to deal with multiple building block operators, 2O and 2D can become versatile and powerful tools to improve the fidelity of digital quantum simulation for many-qubit systems. It could be interesting to further investigate applying our 2D and 2O sequences to other quantum algorithms beyond quantum simulation, such as quantum phase estimation Abrams and Lloyd (1999); Campbell (2019) or to find a good ansatz for the Quantum Approximate Optimization Algorithm Edward Farhi (2014).

Acknowledgements

We thank Chao Yin gbsn(尹超) for helpful discussions. This work was supported in part by DARPA DRINQS, NSF grant EECS1702716, and PHY1734011.

Appendix A Geometrical framework for Trotter error based on holonomy on a fiber bundle

In this section, we provide a geometrical perspective on our problem, based on the notion of holonomy on a fiber bundle Kobayashi and Nomizu (1996).

Consider the space $G\times\mathbb{R}^{2}$ as a trivial bundle over $\mathbb{R}^{2}$, where $G=U(n)$ is the unitary group whose elements are unitary evolutions on the Hilbert space. We define a connection $\omega$ as follows:

[TABLE]

where $A$ and $B$ are $i\times$ the Hamiltonians in our problem. This connection defines a parallel transport over the path $l$ as:

[TABLE]

where $\mathcal{P}$ denotes path-ordering. Under this geometric setting, the error is nothing but the holonomy along a specific loop:

[TABLE]

where the loop $l$ goes from $(0,0)$ to $(p,q)$ along the path defining the Trotter sequence (as shown in Fig. 1 of the main paper), and then goes back to $(0,0)$ along the diagonal.

Geometrically, since the holonomy is provided by the curvature (the topology is trivial and therefore has no contribution), one expects that:

[TABLE]

where $\Omega$ is the curvature 2-form (with a good choice of gauge) and the double integral is over the area surrounded by the loop $l$. Indeed, it was proved in Ref. Erlend and Pierre (2018) that this expectation is true up to at least 3rd-order under a gauge in which the connection form is parallel in radial direction. One can transform any connection $\omega$ into a radially parallel form with a gauge transformation, resulting in:

[TABLE]

In our case,

[TABLE]

and a simple calculation gives (see Eqs. 4-Geometric method to calculate the Trotter error – of the main text) :

[TABLE]

Some nice features appear here:

• •

This formula is true for arbitrary loops, not necessary piecewise linear as those in Fig. 1 of the main paper.

• •

From the first line of Eq. (13), we see that the error can be expressed as an integral of local contributions along the loop. This will make it possible for a dynamic programming algorithm to find the global optimal path.

• •

From the last line of Eq. (13)), the coefficients of $[A,B]$ can be explained geometrically as the area surrounded by the loop, while the coefficients of $[A,[A,B]]$ and $[B,[A,B]]$ are moments about the $x_{2}$-axis and $x_{1}$-axis. Note that the areas and moments here are oriented: positive for anti-clockwise boundary and negative for clockwise boundary; or equivalently (since the path can only go up or right), positive for areas on the right side of the diagonal and negative otherwise.

Appendix B Algorithm for constructing the 2-Optimal sequence

In this section we provide one possible algorithm to find a path that has the minimal third-order error (2O) on a 2-d $p\times q$ grid using dynamic programming. Using the first line of Eq. (13) we can assign a triplet of weights $\{\mathcal{E}_{2},\mathcal{E}_{3A},\mathcal{E}_{3B}\}$ to each edge, where $\mathcal{E}_{2}$ is the second-order error and $\mathcal{E}_{3A},\mathcal{E}_{3B}$ are coefficients of $[A,[A,B]]$ and $[B,[A,B]]$ respectively in the third-order error. For convenience we associate two edges with node $(i,j)$ on the grid. $E_{r,i,j}$ is the vector weight for the edge connecting node $(i-1,j)$ and node $(i,j)$; $E_{u,i,j}$ is the vector weight for edge connecting node $(i,j-1)$ and node $(i,j)$. The algorithm has three steps:

Generate ErrorTable. $ErrorTable(i,j)$ records the (cumulative) vector weight triplets of all possible paths from the origin –node $(0,0)$– to node $(i,j)$. $ErrorTable$ can be generated recursively, See Algorithm 1 2. 2.

Among the paths such that $ErrorTable(p,q)\{1\}=0$, find $MinErr=\min(|ErrorTable(p,q)\{2\}|+|ErrorTable(p,q)\{3\}|)$. 3. 3.

Find the path that corresponds to $MinErr$ by tracing back from node $(p,q)$ to node $(0,0)$.

Assuming $p\sim q\sim n$, the total number of coefficients for the second-order term is $O(n^{2})$, while the number of coefficients for each third-order term is $O(n^{3})$. Therefore, the size of the error set at any point is at most $O(n^{8})$. At each step in the loop, we need to calculate two additions, which cost $O(n^{8})$ computational time, and a union, which also cost at least $O(n^{8})$ time (with some good data structure like union-find or hash table). The total time needed is then $O(n^{10})$. We note that to find the optimal path, we need to optimize over the whole grid: While in the 2D finding algorithm we first normalize $p,q$ by their greatest common divisor, $g$, and repeat $g$ times the optimal construction found for $(p/g,q/g)$, this simplification is not guaranteed to provide the correct solution in 2O. Indeed, one can verify that 2O does not possess translational symmetry; for example, the 2O ordering for (4,4) is not simply repeating the 2O ordering (2,2) twice. Optimizing over all Trotter steps further increases the complexity.

Appendix C Proof that 2-Diagonal path is also the global minimal distance path

In this section, we use mathematical induction to show that any point $A_{n}$ on the 2D path $A_{0}\rightarrow A_{1}\rightarrow\cdots\rightarrow A_{p+q}$ is the closest one to the diagonal among all possible points after $n$ steps. So as a greedy method, 2D finds the global minimal distance path.

To be more concrete, we want to show that for any $n$ and any $x+y=n$, $dist(A_{n},L_{D})\leq dist((x,y),L_{D})$, where $dist(\cdot,\cdot)$ is the distance between two inputs and and $L_{D}$ is the diagonal.

• •

When $n=0$, the statement is easily seen to be true.

• •

Assume the induction hypothesis that $dist(A_{k},L_{D})\leq dist((x,y),L_{D})$, where $x+y=k$, is true for any $k\geq 0$.

• •

From the induction hypothesis, we know that the diagonal $L_{D}$ intersects with $L_{k}:x+y=k$ between [B,C], where B and C are the centers of their unit grids, see Fig. 3. We also know that $L_{D}$ points towards the upper right. So $L_{D}$ intersects with $L_{k+1}:x+y=k+1$ between (D,E), which guarantees that at step $k+1$, either point F or $A_{k+1}$ is the closest to $L_{D}$. Between [D,E] there are two choices and we know that the hypothesis holds for k+1.

By mathematical induction, we know that the statement holds for any natural number n. Fig. 3 shows the least obvious case where the two possible steps following $A_{k}$ lie on the same side of $L_{D}$.

Appendix D Analytical bounds on 2-Diagonal fidelity

Assume that the target evolution is $e^{x(pA+qB)}$, where $p$ and $q$ are mutually prime integers 222If that is not the case and they have a common divider $m$, we can set $p^{\prime}=p/m$, $q^{\prime}=q/m$ and find the construction for $U^{\prime}=e^{x(p^{\prime}A+q^{\prime}B)}$ for $n$ and then repeat it $m$ times. and $x$ is a small parameters.

Approximating the desired evolution with a 2T expansion gives

[TABLE]

or

[TABLE]

where $\mathcal{E}$ is the error operator. For Eq. (14) we have

[TABLE]

and for Eq. (15) we have

[TABLE]

where we defined $C=[A,[A,B]]$ and $D=[B,[A,B]]$. The symmetrized 2D decomposition has a general form

[TABLE]

where $\sum_{k=1}^{M+1}a_{k}=\sum_{k=1}^{M}b_{i}=1$, and $a_{k}\geq 0$, $b_{k}\geq 0$.

We can evaluate the error by starting from the middle of the sequence and keeping track of the third-order error:

[TABLE]

The error $\mathcal{E}_{1}$ can be written as

[TABLE]

Consider the next layer

[TABLE]

where

[TABLE]

In general, we can find the formulas for the odd/even errors:

[TABLE]

where $\mathcal{A}_{k}=\sum_{i=1}^{k}a_{i}$ and $\mathcal{B}_{k}=\sum_{i=1}^{k}b_{i}$. We want to compare the total error of 2D, $\mathcal{E}_{\textrm{2D}}=\sum_{k=1}^{M}(\mathcal{E}_{2k}+\mathcal{E}_{2k-1})$ with the conventional 2nd Trotter error $\mathcal{E}$ and $\mathcal{E}^{\prime}$. We use the Frobenius norm to compare operators, $\|E\|^{2}_{F}=\textrm{Tr}(E^{\dagger}E)$, and without loss of generality we assume $\|C\|_{F}=\|D\|_{F}$, and to simplify the presentation we also take $\textrm{Tr}(C^{\dagger}D)=0$. If $\|C\|_{F}\neq\|D\|_{F}$ we can rescale $p$ to make $\|C\|_{F}=\|D\|_{F}$. Hereafter we can drop the subscript $F$ without causing confusion. With these assumptions, it is straightforward to show that $\|\mathcal{E}_{\textrm{2D}}\|\leq\|\mathcal{E}\|$ iif $p^{2}\mathcal{C}^{2}+q^{2}\mathcal{D}^{2}\leq\text{min}(p^{2}+4q^{2},q^{2}+4p^{2})$, where

[TABLE]

Note that $\mathcal{B}_{k-1}$ only exists when $k\geq 2$. Indeed, we have

[TABLE]

If we rewrite Eq. (27) in form of integration, Eq. (13) is regained. Define a function $f(a_{1},\cdots,a_{M+1},b_{1},\cdots,b_{M})=p^{2}\,\mathcal{C}^{2}+q^{2}\,\mathcal{D}^{2}$. We can show that any choice of the coefficients $a_{i}$, $b_{i}$ that does not collapses the 2D expansion onto the two 2nd Trotter ones gives a lower error than $\text{max}(\|\mathcal{E}\|,\|\mathcal{E}^{\prime}\|)$. The strategy is to find the upper bound of $f$ by varying $a_{i}$ and $b_{i}$ recursively.

We can start by varying the $a_{M+1}$ parameter. First change $a_{1},a_{2},\cdots,a_{M},a_{M+1}$ to a set of independent parameters $\mathcal{A}_{1},\mathcal{A}_{2},\cdots,\mathcal{A}_{M-1},a_{M}=1-a_{M+1}-\mathcal{A}_{M-1},a_{M+1}$. After this change $a_{M+1}$ is an independent variable. All the terms in $\mathcal{C}$ that depend on $a_{M+1}$ are $3b_{M}(1-a_{M+1})^{2}-\mathcal{B}_{M}+\mathcal{A}_{M-1}(\mathcal{B}_{M}-b_{M})$ therefore

[TABLE]

$\mathcal{C}(a_{M+1})$ is a decreasing function and $\mathcal{C}(a_{M+1}=0)=2$, $\mathcal{C}(a_{M+1}=1)=-1$. Actually $C(a_{M+1})=3b_{M}a_{M+1}^{2}-6b_{M}a_{M+1}+2$.

Similarly we find

[TABLE]

$\mathcal{D}(a_{M+1})$ is a decreasing function and $\mathcal{D}(a_{M+1}=0)=1$, $\mathcal{C}(a_{M+1}=1)=-2$. Actually $\mathcal{D}(a_{M+1})=3b_{M}(b_{M}-2)a_{M+1}+1$.

The derivatives of the function $f(a_{1},\cdots,a_{M+1},b_{i},\cdots,b_{M})=p^{2}\,\mathcal{C}^{2}+q^{2}\,\mathcal{D}^{2}$ are given by

[TABLE]

Plugging Eq. (28),(29),(30),(31) into (33) we find that when $\mathcal{C}>0$, $f$ is a convex function of $a_{M+1}$ and we have $f\leq\text{max}\{f(a_{M+1}\!=\!0),f(a_{M+1}\!=\!1)\}$. When $\mathcal{C}\leq 0$, it is not difficult to verify that when $\mathcal{C}\leq 0$, we also have $\mathcal{D}\leq 0$. Even if $\frac{\partial^{2}f}{\partial a_{M+1}^{2}}\leq 0$, we still have $\frac{\partial f}{\partial a_{M+1}}\geq 0$ and $f\leq\text{max}\{f(a_{M+1}\!=\!0),f(a_{M+1}\!=\!1)\}$ is still true.

When $a_{M+1}=1$ (that is, all other coefficients $a_{1},\cdots,a_{M}$ are zero, meaning that $M=0$) $\mathcal{C}=-1$, $\mathcal{D}=-2$ and $f=p^{2}+4q^{2}$, as expected, since we retrieve the simpler Trotter expansion $e^{\frac{1}{2}xpA}e^{xqB}e^{\frac{1}{2}xpA}$.

When $a_{M+1}=0$, we can repeat the same analysis above, now taking the derivatives with respect to $b_{M}$. We find again that for $b_{M}=1$ the construction collapses to the Trotter expansion $e^{\frac{1}{2}xqB}e^{xpA}e^{\frac{1}{2}xqB}$, and we can analyze the case $b_{M}=0$ as done above. Recursively, we thus find that $f\leq 4p^{2}+q^{2}$, with equality in the case where the construction simplifies to the Trotter expansion $e^{\frac{1}{2}xqB}e^{xpA}e^{\frac{1}{2}xqB}$ or $f\leq p^{2}+4q^{2}$, with equality in the case where the construction simplifies to the Trotter expansion $e^{\frac{1}{2}xpA}e^{xqB}e^{\frac{1}{2}xpA}$.

When $p=q$, this proves that any decomposition of the exponential operators with more and smaller steps than 2T yields a smaller error than 2T.

When $p\neq q$, the construction following the form in Eq. (18) has an error that is bounded by $\text{max}(\|\mathcal{E}\|,\|\mathcal{E}^{\prime}\|)$. We can further show that the error of 2D is actually bounded by $\text{min}(\|\mathcal{E}\|,\|\mathcal{E}^{\prime}\|)$. When $p>q$, $\|\mathcal{E}\|<\|\mathcal{E}^{\prime}\|$ and we want to show that $f\leq p^{2}+4q^{2}$. According to the construction of 2D, when $p>q$, we first take $m$ steps along $A$, where $\frac{mq}{p+q}<\frac{1}{2}$ and $\frac{(m+1)q}{p+q}>\frac{1}{2}$, i.e., $m=\lfloor\frac{p+q}{2q}\rfloor$. The construction that is closest to $e^{\frac{q}{2}xB}e^{xpA}e^{\frac{q}{2}xB}$ allowed by 2D writes

[TABLE]

For the decomposition in Eq. (34), we find that $\mathcal{C}=3(1-\frac{2m}{p})^{2}-1$ and $\mathcal{D}=1-\frac{6m}{p}$. It’s not difficult to verify that $p^{2}\mathcal{C}^{2}+q^{2}\mathcal{D}<p^{2}+4q^{2}$, when $p>q$. If we further decompose the central part $e^{\frac{q}{2}B}e^{(1-\frac{2m}{p})pA}e^{\frac{q}{2}B}$, since $q\geq 1-\frac{2m}{p}$, the worst case will converge to $e^{(\frac{1}{2}-\frac{m}{p})pA}e^{qB}e^{(\frac{1}{2}-\frac{m}{p})pA}$, with which the whole formula becomes $e^{\frac{1}{2}pA}e^{qB}e^{\frac{1}{2}pA}$. In other words, however you decompose the central part in Eq. (34), the error of the whole formula is bounded by $p^{2}+4q^{2}$. In summary, when $p>q$, if we construct the nested structure according to the principle of 2D, the error will be smaller than $p^{2}+4q^{2}$. Similarly when $p<q$, $f\leq q^{2}+4p^{2}$. In summary $f\leq\text{min}(p^{2}+4q^{2},q^{2}+4p^{2})$ and thus $\|\mathcal{E}_{2D}\|\leq\text{min}(\|\mathcal{E}\|,\|\mathcal{E}^{\prime}\|)$.

Appendix E Extension to three operators

Using the geometric picture, 2D can be easily extended to three or more operators. The greedy algorithm always find the step that keeps the trajectory closest to the diagonal. In principle one could also use dynamic programming to find a path that minimizes the third-order error. Nevertheless, as the number of non-commutative operators increases, the complexity of the algorithm increases fast, while 2D remains efficient. Here we tested the performance of 2D on a two-qubit Ising Hamiltonian with both transverse and longitudinal fields. The three implementable Hamiltonians read

[TABLE]

and the target evolution is $e^{-i(pH_{1}+qH_{2}+rH_{3})t}$, where $p,q,r$ are integers with greatest common divisor 1.

We compare in Fig. 4(a) the results for 2D and the first-order Trotter expansion (for simplicity, we did not symmetrize either 2D or first-order Trotter, so that both are now first-order approximations). For case $p=3,~{}q=4,~{}r=5$, first-order Trotter expansion is $(e^{-i3H_{1}\frac{t}{n}}e^{-i4H_{2}\frac{t}{n}}e^{-i5H_{3}\frac{t}{n}})^{n}$ and 2D is

[TABLE]

For case $p=6,~{}q=4,~{}r=2$, first-order Trotter expansion is $(e^{-i3H_{1}\frac{t}{n}}e^{-i2H_{2}\frac{t}{n}}e^{-i1H_{3}\frac{t}{n}})^{2n}$ and 2D expansion is

[TABLE]

The trajectories of the two cases are shown in Fig. 4(b). We see that 2D outperforms the first-order Trotter expansion for all $(p,q,r)$s up to $p+q+r=12$. Similar results can be obtained with different Hamiltonians.

Appendix F Number of gates and number of Trotter steps

While 2D reaches better fidelity than 2T, it generally does so by increasing the number of alternations between the two (or more) building blocks $e^{A}$ and $e^{B}$. When simulating spin dynamics, the total number of switches and gates is not so crucial since the number of Hamiltonians is not large, and the total simulation time remains the same for different Trotter expansions. In many experimental implementations, switching between the two Hamiltonians comes at a cost, either because of additional time overheads, or because the switch itself, obtained, e.g., by applying a short control pulse, might introduce control errors. Still, large number of pulses have routinely been applied, e.g., in dynamical decoupling Abobeih et al. (2018); Bar-Gill et al. (2013) and Hamiltonian engineering Wei et al. (2018), so a large number of switchings is not infeasible. Quantifying the decrease in fidelity due to such errors is complex, because it would be platform dependent. Still, we can evaluate the number of switchings required to achieve a certain fidelity for different constructions. We find that 2D does not always impose an overhead in the number of switches with respect to, e.g., the 2T expansion, although this does depend on the particular $(p,q)$ pair considered (see Fig. 5). In particular, for more “unbalanced” pairs ($p\ll q$ or vice-versa), 2D can outperform 2T.

When the elementary unitary evolutions are composed by a set of discrete gates that cannot be merged Jones et al. (2012), the number of switches no longer quantify the amount of resources needed. The total number of Trotter steps, which is proportional $(p+q)n$, best describes the resources needed. The lower panels in Fig. 5 shows that to reach the same fidelity, 2D requires a smaller number of Trotter steps, thus providing an advantage with respect to this metric.

Appendix G Performance when the time-step resolution is not limited

2D outperforms 2T under the constraints that only integer number of elementary gates are implementable. For example $e^{A/n}$ is implementable but $e^{1.2A/n}$ is not allowed. The constraints could come from experimental limitations on both the shortest time step and the time step resolution. Often in practical implementations, the duration of each Trotter step, even when it is above the smallest time step, cannot be set at will, but it is typically an integer multiple of a minimum time duration. We have assumed that this minimum time duration is the smallest Trotter step. Qubits with the fastest dynamics suffer from these two limitations since typically the timing of each operation is set by a digital clock with a minimum resolution. In addition, often, a single Trotter step is given by a fixed (fault-tolerant) quantum gate (such as CNOT gate), thus imposing an effective minimum “time resolution” constraint. However, for qubits that present a slow dynamics, the time-resolution constraint might be negligible. In this case, assuming that there are no constraints on time resolution, 2T typically performs better than 2D that by construction assumes a finite time-step resolution (equal to the minimum time step). In Fig. 6 we compare 2T and 2D constructions for different $p,~{}q$. Note that for different $p$, $q$ values, there is an optimal 2T expansion, which is the one we selected. Still, we see that 2D has better performance than 2T for certain $p$, $q$ values.

Appendix H Comparison with other Trotter expansions

Higher order Trotter expansions have the general form of $e^{pA+qB}\approx e^{c_{1}pA}e^{c_{2}qB}e^{c_{3}pA}\cdots e^{c_{M}qB}$, where however the coefficients $pc_{1,3,\cdots}$ and $qc_{2,4,\cdots}$ are usually not integers and can even be negative and thus typically challenging (or impossible) to implement experimentally. One example of third order Trotter expansion is Ruth’s formula Ruth (1983),

[TABLE]

where negative coefficients are used. Indeed, some negative coefficients are unavoidable in third- and higher-order Trotter expansions Suzuki (1991). Implementing “time reversal” to construct the unitaries with negative coefficients is not always possible. In Fig. 7 we show that when Trotter step is not very small ($t/n=0.1$), 2nd-Near-optimal has better fidelity than the conventional first-, second-, and third-order Trotter expansions. In the numerical evaluation, we in total $(p+q)n=100$ Trotter steps are used. For 100 steps it is difficult to find 2O path so we only focus on 2D.

While the comparison to well-known Trotter decomposition is instructive, one might think that simple alternative constructions could be found heuristically. For example, the intuition behind 2D is to maximize the alternation (spreading out) of the two operators, $e^{A/n}$ and $e^{B/n}$. For many $(p,q)$ pairs, however, it is not simple to find an exact alternating ordering. Then, a naive alternation construction is still outperformed by 2D, see, e.g., Fig. 8 for several examples.