Simulating the dynamics of a quantum system is one of the most natural applications of a quantum computer. Indeed, the idea of quantum computation, as suggested by Feynman  and others, was primarily motivated by the problem of quantum simulation. Lloyd gave the first explicit quantum algorithm for simulating local Hamiltonians , and subsequent work studied the simulation of general sparse Hamiltonians [2, 8]. Quantum computers can simulate a variety of physical systems, including quantum chemistry [4, 53, 43, 15], quantum field theory [29, 30], and many-body physics , and could ultimately lead to practical applications such as designing new pharmaceuticals, catalysts, and materials [5, 19].
The original quantum simulation algorithm proposed by Lloyd makes use of the Lie-Trotter product formula. Given and a Hamiltonian with , the Lie-Trotter formula
approximates the ideal dynamics
in the sense that . With a finite value of , we obtain the approximation . In other words, is accurate to first order in , so we call it a first-order formula. Suzuki systematically constructed ()th-order formulas
with . We refer to all such formulas as product formulas. The product-formula approach has been studied extensively (see for example [44, 21, 34, 32, 22]), and is a natural choice for experimental simulations [13, 6, 35] due to its simplicity and the fact that it does not require any ancilla qubits.
The major difficulty in making the product-formula algorithm concrete is to choose so that the simulation error is at most some allowed threshold . Specifically, to simulate the Hamiltonian for time , rigorous analysis suggests choosing for the first-order formula and for the ()th-order formula. Within each segment, the number of elementary operations scales linearly in , giving and as the gate complexities for the first- and ()th-order formulas, respectively. These asymptotic complexities are improved by more recent algorithms [11, 9, 10, 38, 39, 37, 12, 14, 31, 40]. However, product formulas appear to be empirically advantageous for simulating various practical systems [17, 4, 46], performing significantly better than the rigorous bounds suggest. This dramatic gap between the provable and the actual behavior of product formulas suggests that it may be possible to significantly tighten their error analysis.
A natural class of Hamiltonians that includes many physically reasonable systems is the class of lattice Hamiltonians [29, 41, 25, 24]. Lattice Hamiltonians arise in many models of condensed matter physics, including systems of spins (e.g., Ising, XY, and Heisenberg models; Kitaev’s toric code and honeycomb models; etc.), fermions (e.g., the Hubbard model and the - model), and bosons (e.g., the Bose-Hubbard model). Note that fermion models can be simulated using local interactions among qubits by using a mapping that preserves locality . Digital simulations of quantum field theory also typically involve approximation by a lattice system .
For simplicity, we focus here on nearest-neighbor lattice systems in one dimension. In this case, qubits are laid out on a one-dimensional lattice and the Hamiltonian only involves nearest-neighbor interactions. A formal definition is given as follows.
Definition 1 (Lattice Hamiltonian).
Let be an integer. We say that a Hamiltonian is a lattice Hamiltonian if it acts on qubits and can be decomposed as , where each is a Hermitian operator that acts nontrivially only on qubits and . We assume that , for otherwise we evolve under the normalized Hamiltonian for time .
Assuming that is fixed, product formulas can simulate such a lattice system with gate complexity in the first-order case and in the ()th-order case. However, it is natural to expect a more efficient simulation. Roughly speaking, a system simulates its own evolution for constant time using only constant circuit depth—and hence an extensive number of gates—so one might expect a true simulation complexity of . Indeed, Jordan, Lee, and Preskill claimed that the product-formula algorithm can simulate an -qubit lattice system with gates , but they did not provide rigorous justification and it is unclear how to formalize their argument. Subsequent work improves the analysis of the product-formula algorithm using information about commutation among terms in the Hamiltonian [48, 17, 43, 53], the distribution of norms of terms , and by randomizing the ordering of terms [18, 56]. However, to the best of our knowledge, none of these improvements achieves the gate complexity for lattice simulation.
Recently, Haah, Hastings, Kothari and Low proposed a new algorithm for lattice simulation . Instead of analyzing the product-formula approach, they develop an approach motivated by the Lieb-Robinson bound [42, 28], which quantifies how fast information can propagate in a system with local interactions. This approach decomposes the entire evolution into blocks, where each block involves forward and backward evolution on a small region. If each block is then simulated using the approach of  or , this gives an algorithm with complexity , nearly matching a lower bound of that Haah et al. also establish. On the other hand, if each block is simulated with a high-order product formula, we obtain a lattice simulation algorithm with total gate complexity that does not require ancilla qubits, also nearly matching the lower bound.
However, this does not address the long-standing claim about product-formula simulation . Indeed, product formulas can harness the commutation of terms without the help of a Lieb-Robinson decomposition. Pinning down this capability for lattice systems is not only a natural challenge in its own right, but could also illuminate the simulation of other systems, where an effective Lieb-Robinson decomposition may be difficult to find. Furthermore, the extensive prior results on the experimental performance of product formulas are not directly applicable to the Lieb-Robinson-based algorithm.
This paper studies the extent to which product formulas can harness the commutation of terms in a lattice Hamiltonian. We consider canonical product formulas with stages, order , and a uniform upper bound on the coefficients (as defined in Section 3), which models well-known constructions including the Lie-Trotter-Suzuki formulas. We prove that a th-order formula can simulate an -qubit lattice Hamiltonian with accuracy using
elementary operations, assuming that , , and are constant. This resolves the open question about product-formula simulation by choosing a sufficiently large and fixing to be a constant. Combining with the previous lower bound of gates , our result implies that product formulas can simulate lattice Hamiltonians with nearly optimal gate complexity, yielding a simpler approach than the previous algorithm .
The key technique that allows us to achieve the gate complexity is an integral representation of the error that we develop based on Descombes and Thalhammer’s “local error analysis” of product formulas . In the local error representation, the integrand is expressed as a linear combination of nested commutators, where the number of summands and the number of layers of nesting are both independent of and . We use this representation to get the correct scaling of the gate complexity as a function of and . In contrast, the conventional approach uses the Baker-Campbell-Hausdorff formula or naive Taylor expansion, which requires the manipulation of infinite series and appears to be technically challenging to analyze .
To illustrate the proof idea, we show how to obtain the gate complexity using the first-order formula, i.e., the Lie-Trotter formula. In particular, we prove that this complexity is achievable when the terms of the lattice Hamiltonian are simulated using the even-odd
However, this ordering is not essential for the first-order formula. We prove that the first-order simulation is robust with respect to the ordering of the summands: no matter how terms from the lattice Hamiltonian are ordered, the asymptotic gate complexity is always . These results on the first-order product formula are presented in Section 2.
Analyzing higher-order formulas is more challenging. To address general product formulas, we represent them in a canonical form, along the lines of Descombes and Thalhammer’s work . Canonical product formulas are easy to manipulate and encompass well-known constructions such as the Lie-Trotter-Suzuki formulas as special cases. As pointed out by [20, 3], the error of a product formula admits an integral representation that obeys certain order conditions (i.e., conditions under which the formula is accurate to a given order). Canonical product formulas and their order conditions are discussed in Section 3.
To best harness the order conditions, we explore the commutator structure of the integrand. Specifically, we express the integrand as a linear combination of nested commutators, where both the number of summands and the number of layers of nesting are independent of and . These technical details are discussed in Section 4. Previous work of Descombes and Thalhammer achieves a similar goal by introducing auxiliary functions  whose recursive structure is hard to unravel. Our simplified analysis avoids this issue. We also identify and correct mistakes in the analysis of .
We give error bounds and analyze the asymptotic gate complexity of higher-order formulas in Section 5. In particular, we show that a th-order canonical formula has gate complexity when the terms in the lattice are ordered in the even-odd pattern, resolving the aforementioned claim about product-formula simulation .
For time-dependent Hamiltonians, there is no closed-form formula for the evolution operator and the analysis is significantly more difficult. Nevertheless, we show in Section 6 that a time-dependent lattice Hamiltonian can be simulated using product formulas with nearly optimal gate complexity. Along the way, we establish order conditions for time-dependent Hamiltonian simulation and develop a corresponding local error representation, which may be of independent interest.
The Hamiltonian introduced in Definition 1 is a one-dimensional lattice system with open boundary conditions. In Section 7, we discuss how our analysis can be adapted to handle periodic boundary conditions, constant-range interactions, and -dimensional systems with .
We have primarily worked with an idealized case where quantum gates can be performed perfectly. However, in an experimental realization, the performance of product formulas may be affected by noise. We analyze the trade-off between the algorithmic error and the gate error in Section 8.
As noted in , it may be difficult to apply the simulation algorithm based on the Lieb-Robinson decomposition to Hamiltonians whose interactions are described by general graphs: in such cases, an effective Lieb-Robinson decomposition may be hard to find. Our improved product-formula algorithm harnesses the commutation of Hamiltonian terms directly without the help of Lieb-Robinson bounds. For future work, it might be fruitful to use product formulas to speed up Hamiltonian simulation for systems that lack spatial locality .
Our local error analysis represents the product-formula error as an integral of a linear combination of nested commutators. Interestingly, similar techniques have been used to establish the Lieb-Robinson bound and to study the computational complexity aspects of local Hamiltonians and ground states [42, 28, 50, 1]. We leave it as an avenue for future research to explore whether these techniques could find more applications in the study of locality of quantum systems.
Finally, we show that the complexity is achievable by ordering the terms in the even-odd pattern. It would be interesting to prove (or disprove) that there exist other possible orderings giving the nearly optimal gate complexity. It would also be of interest to better understand the empirical behavior of these algorithms.
2 Complexity of the first-order algorithm and ordering robustness
In this section, we consider the problem of simulating a lattice Hamiltonian using the first-order product formula. We prove in Section 2.1 that, by ordering Hamiltonian terms in the even-odd pattern, the first-order formula can simulate an -qubit lattice system for time using elementary gates. However, this particular ordering of terms is not essential. We proceed to show in Section 2.2 that the ordering-robustness property holds in the first-order case: no matter how terms from the Hamiltonian are ordered, the first-order algorithm always has the same asymptotic complexity .
2.1 Analysis of the first-order algorithm
We first consider the case where the Hamiltonian is the sum of two terms , where and are time-independent Hermitian operators. The analysis for such a Hamiltonian will later be applied to the simulation of a lattice Hamiltonian. The evolution of for time is given by
whereas the first-order product formula is
Throughout this paper, we assume that , for otherwise we absorb the minus sign into the Hermitian operators , , and . We quantify the error by computing the spectral-norm distance . The conventional approach to error analysis uses the Baker-Campbell-Hausdorff formula or naive Taylor expansion, which seems unlikely to produce an asymptotically tight bound. Instead, our analysis uses an integral representation of the error operator that we develop based on Descombes and Thalhammer’s “local error analysis” .
The following lemma will be useful in our analysis.
Let , let be an operator, and let be a continuous operator-valued function of . Then the inhomogeneous operator equation
has a unique solution given by Lagrange’s variation-of-constants formula
See for example [33, Proposition 4.8 and Theorem 4.9]. ∎
To apply this lemma, we take the first time derivative of and obtain
Applying Lemma 1 with initial condition , we find an integral representation of the error
To see how well the first-order formula approximates the ideal evolution , it suffices to study the commutator .
We differentiate this commutator with respect to and get
Invoking Lemma 1 with initial condition , we find an integral representation for the commutator
Altogether, we have
which completes the local error analysis for the first-order formula.
We now apply the above analysis to the simulation of lattice Hamiltonians. Without loss of generality, we assume that is even. Let
be an -qubit lattice Hamiltonian as defined in Definition 1. We collect the terms of into and , where
Note that () can further be decomposed into terms such that no two terms have overlapping supports. At the same time, the support of each term from overlaps with any term from by only qubit.
We instantiate and in (14). The result is that
The left-hand side of the above equation is the approximation error of the first-order product formula, when the terms are ordered in the even-odd pattern
The analysis of the right-hand side requires more effort.
We expand and according to their definitions. We fix an arbitrary term in . For this particular term, the commutator is non-zero only when . If we adopt the convention , we find that
The above expansion is sufficient for the analysis of the first-order formula. However, we will make another step of simplification that is instructive for the analysis of higher-order formulas.
We note that terms in all have non-overlapping supports and therefore commute with each other, so
Observe that for , the exponentials and do not overlap with the commutator . Therefore, these exponentials can commute through the commutator and cancel in pairs. The result is that
Altogether, we have
We have thus established the following theorem.
Theorem 2 (Local error bound for the first-order algorithm).
Let be a lattice Hamiltonian. Order its terms in the even-odd pattern , so that the evolution of for time is represented by and the first-order formula is . Assume that . Then
The above error bound works well when is small. To simulate for a longer time, we divide the entire evolution into segments, and within each segment, we simulate the evolution using the first-order formula. The above analysis shows that the error of the entire simulation can be bounded by
Upper bounding the right-hand side by and solving for , we find that it suffices to take . This implies that the gate complexity of the first-order product-formula algorithm is if we simulate with constant accuracy.
2.2 Ordering-robustness of the first-order algorithm
Our error bound in the previous section works when terms of the lattice Hamiltonian are ordered in the even-odd pattern. However, this choice is not necessary. In this section, we show that the first-order product-formula algorithm has asymptotic gate complexity with respect to any ordering of the lattice terms.
Our analysis relies on the following swapping lemma.
Lemma 3 (Local error bound for swapping exponentials).
Let be nonnegative and let , be Hermitian operators. Then
Corollary 4 (Local error bound for swapping lattice terms).
Let be nonnegative and be a lattice Hamiltonian with . Then for any ,
Let be a lattice Hamiltonian. We now simulate this Hamiltonian using the first-order formula, but allow terms to be ordered arbitrarily as
where is a permutation on elements. Then the spectral-norm error is upper bounded by
A bound for the second term is given in Theorem 2. To handle the first term, we transform to by swapping neighboring exponentials. Every time two exponentials and are swapped, we use the triangle inequality and Corollary 4 to bound the error. The total number of swaps of exponentials and with is at most , incurring error . We have therefore established the following theorem.
Theorem 5 (Ordering robustness of the first-order algorithm).
Let be a lattice Hamiltonian, and order its terms according to a permutation , so that the evolution of for time is represented by , and the first-order formula is . Assume that . Then
Thus, we find that the asymptotic gate complexity of the first-order algorithm is always , no matter how the lattice terms are permuted. We call this phenomenon ordering robustness. The above theorem shows that the first-order product-formula algorithm is ordering-robust. Whether a similar property holds for a general higher-order formula remains an open question.
3 Canonical product formulas and order conditions
In this section, we introduce notation and terminology that is useful for studying higher-order formulas. We study canonical product formulas (denoted by ), which include well-known constructions such as the Lie-Trotter-Suzuki formulas. We use to denote the error of the derivative of a formula, from which we construct an integral representation of the product-formula error. To obtain the correct dependence on , we further represent as the product of and another operator that we denote by . Finally, we study the order conditions of , , and . We show that the order conditions of a product formula are automatically satisfied. (Similar order conditions are sketched in , but we discuss them here for completeness.) Whenever possible, we follow the notation and terminology of .
3.1 Canonical product formula and the operator
Similar to the first-order case, it is instructive to study a setting where the Hamiltonian is the sum of two Hermitian terms . The evolution of for time is given by the unitary operator , which may then be simulated using a specific product-formula algorithm, such as the Lie-Trotter formula or the Suzuki formulas introduced in Section 1. We will not analyze these formulas case-by-case. Instead, we consider canonical product formulas, a universal concept that includes well-known constructions.
Definition 2 (Canonical product formula).
Let be a Hamiltonian consisting of two terms , where and are Hermitian operators. We say that an operator-valued function is a canonical product formula for if it has the form
where and are real coefficients. The parameter denotes the number of stages, and is the th-stage operator for . We let be an upper bound on the coefficients, i.e.,
Finally, we say that the product formula has order for some integer if
We call an ()-formula if we need an explicit description of the parameters.
Although common constructions of product formulas involve stages where exponentials can be ordered both as and as
, we can achieve such orderings by padding with identity operators. In particular, we now show in detail how some well-known constructions of product formulas can be recast in the canonical form.
Example 1 (First-order formula).
The first-order formula may be represented as a -stage canonical formula by setting , whereas its reversed version is a -stage canonical formula with the choice .
Example 2 (Second-order formula).
The second-order formula may be represented as a -stage canonical formula by setting , whereas its reversed version is a -stage canonical formula with the choice .
Example 3 (()th-order formula).
The ()th-order Suzuki formula defined in (3) is an ()-formula, where , , and .
To study the order condition of a product formula, we need the following lemma.
Let be an operator-valued function that is infinitely differentiable. Let be a nonnegative integer. The following two conditions are equivalent.
Asymptotic scaling: .
Derivative condition: .
Condition implies by Taylor’s theorem. Assuming Condition holds, we must have that
for some (and for sufficiently small). Let be the first integer such that . We use Taylor’s theorem to find such that
We combine the above inequalities and divide both sides by . Taking the limit gives us a contradiction. ∎
By definition, a product formula has order for some integer if
holds for any . Invoking Lemma 6, we find an equivalent order condition
3.2 Integral representation of error and the operator
As in Section 2, we seek an integral representation of the product-formula error . To this end, we differentiate and rewrite the derivative as , where
Recall that is accurate up to order . Therefore, and we obtain the initial value problem . The solution of this problem is given by the variation-of-constants formula in Lemma 1, from which we obtain an integral representation of the product-formula error
We now determine an order condition for the operator . Since is at least first-order accurate, we have and therefore . By taking derivatives iteratively, one can show that
Conversely, if higher-order derivatives of satisfy the above condition, we must have
3.3 Correct -dependence and the operator
Our goal is to prove a variant of Theorem 2 for higher-order formulas: for a th-order formula with terms ordered in the even-odd pattern, we want the product-formula error to scale like . The integral representation established in (42) is a first step toward this goal. Combining (42) with (43) and using Taylor’s theorem, it can be shown that the product-formula error scales like , achieving the desired dependence on . However, this approach fails to give the correct -dependence. Indeed, a naive application of Taylor’s theorem will produce products of the operators and , which prevents the error from scaling linearly in .
A solution to this problem is to rewrite the integrand using the product formula and another operator. Specifically, we let be the operator such that
In quantum simulation, the product formula is unitary and therefore . However, we will see that has significantly richer structure than it might seem. Analyzing the combinatorial structure of will be the central topic of the next section. For now, we shall focus on its order condition.
We claim that (43) is equivalent to the order condition
By the general Leibniz rule
so (46) implies (43). We prove the converse by induction. For , we have and . Therefore, implies that . Assume that has been proved for all . We apply the general Leibniz rule to compute the th-order derivative of and find
Therefore for all .
3.4 Theorem of order conditions
We now summarize all the product-formula order conditions determined above in the following theorem.
Theorem 7 (Order conditions for canonical product formulas).
Let be a Hamiltonian consisting of two terms , where and are Hermitian operators. Let be an integer and let <