Hamiltonian simulation in the low-energy subspace

We study the problem of simulating the dynamics of spin systems when the initial state is supported on a subspace of low energy of a Hamiltonian $H$. This is a central problem in physics with vast applications in many-body systems and beyond, where the interesting physics takes place in the low-energy sector. We analyze error bounds induced by product formulas that approximate the evolution operator and show that these bounds depend on an effective low-energy norm of $H$. We find improvements over the best previous complexities of product formulas that apply to the general case, and these improvements are more significant for long evolution times that scale with the system size and/or small approximation errors. To obtain these improvements, we prove exponentially decaying upper bounds on the leakage to high-energy subspaces due to the product formula. Our results provide a path to a systematic study of Hamiltonian simulation at low energies, which will be required to push quantum simulation closer to reality.


I. INTRODUCTION
The simulation of quantum systems is believed to be one of the most important applications of quantum computers [1]. Many quantum algorithms for simulating quantum dynamics exist [2][3][4][5][6][7][8][9][10][11], with applications in physics [12,13], quantum chemistry [14][15][16], and beyond [17]. While these algorithms are deemed efficient and run in time polynomial in factors such as system size, ongoing work has significantly improved the performance of such approaches. These improvements are important to explore the power of quantum computers and push quantum simulation closer to reality.
Leading Hamiltonian simulation methods are based on a handful of techniques. A main example is the product formula, which approximates the evolution of a Hamiltonian H by short-time evolutions under the terms that compose H [4,5,18,19]. Each such evolution can be decomposed as a sequence of two-qubit gates [12] to build up a quantum algorithm. Product formulas are attractive for various reasons: they are simple, intuitive, and their implementations may not require ancillary qubits, which contrasts other sophisticated methods as those in Refs. [7,8]. Product formulas are also the basis of classical simulation algorithms including pathintegral Monte Carlo [20].
Recent works provide refined error bounds of product formulas [21][22][23][24]. These works regard various settings, such as when H is a sum of spatiallylocal terms or when these terms satisfy Lie-algebraic properties. Nevertheless, while these improvements are important and necessary, a number of shortcomings remain. For example, the best-known complexities of product formulas scale poorly with the norm of H or its terms, which can be very large or unbounded, even when the evolved quantum system does not explore high-energy states. These complexities may be improved under physically-relevant assumptions on energy scales. In fact, numerical simulations of few spin systems suggest that product formulas applied to low-energy states lead to much lower errors than that of the worst-case. Figure 1, for example, shows these errors for a 2x6 spin-1/2 Heisenberg model, suggesting that a complexity improvement is possible under a low-energy assumption on the initial state. Simulation results for related models present similar features. Nevertheless, our inability of simulating larger quantum systems with classical computers efficiently demands for analytical tools to actually demonstrate strict improvements on complexities of product formulas that apply generally.
To this end, we investigate the Hamiltonian simulation problem when the initial state is supported on a low-energy subspace. This is a central problem in physics that has vast applications, including the simulation of condensed matter systems for studying quantum phase transitions [25], the simulation of quantum field theories [13], the simulation of adiabatic quantum state preparation [26,27], and more. We analyze the complexities of product formulas in this setting and show significant improvements with respect to the best known complexity bounds that apply to the general case. i,j XiXj + YiYj + ZiZj, where Xi, Yi, and Zi are the Pauli operators for the ith spin, and H = H1 + H2, where H1 and H2 are the interaction terms represented by blue and red bonds, respectively. (a) The evolution operator for time s, U (s) = e −isH , is approximated by the first order product formula W1(s) = e −isH 1 e −isH 2 . The plot shows the largest approximation errors when acting on various low-energy subspaces associated with increasing energies, labeled by n = 1, 50,150,200, and in the worst case. (b) Similar results for when the evolution operator U (s) = e −isH is approximated by the second order product formula W2(s) = e −isH 1 /2 e −isH 2 e −isH 1 /2 .

A. Overview
Our main result is that, for a local Hamiltonian on N spins H = l H l with H l ≥ 0, the error induced by a p-th order product formula is O((∆ s) p+1 ), where s is a (short) time parameter and ∆ is an effective low-energy norm of H. This norm depends on ∆, which is an energy associated with the initial state, but also depends on s and other parameters that define H. The best known error bounds for product formulas that apply to the general case depend on the H l 's [23]. (Throughout this paper, . refers to the spectral norm.) Thus, an improvement in the complexity of product formulas is possible when ∆ max l H l , which can occur for sufficiently small values of ∆ and s. Such values of s appear in low-order product formulas (e.g., first order) or, for larger order, when the overall evolution time t is sufficiently large and/or the desired approximation error ε is sufficiently small. We summarize Order Previous result Low-energy simulation TABLE I. Improvements of low-energy simulation: Comparison between the best-known complexity [23] and the complexity of low-energy simulation for p-th order product formulas. Results show the Trotter number for constant ∆ and local Hamiltonians on N spins with constant degree and strength bounded by J, and τ = |t|J. ε is the approximation error. TheÕ notation hides polylogarithmic factors in τ /ε. some of the complexity improvements in Table I. To obtain our results, we introduce the notion of effective Hamiltonians that are basically the H l 's restricted to act on a low-energy subspace. The relevant norms of these effective operators is bounded by ∆ . One could then proceed to simulate the evolution using a product formula that involves effective Hamiltonians and obtain an error bound that matches ours. A challenge is that these effective Hamiltonians are generally non-local and difficult to compute. Methods such as the local Schrieffer-Wolff transformation [28,29] work only at the perturbative regime and numerical renormalization group methods for spin systems [30,31] have been studied only for a handful of models, while a general analytical treatment does not exist. Thus, efficient methods to simulate time evolution of effective Hamiltonians are lacking. We address this challenge by showing that evolutions under the effective Hamiltonians can be approximated by evolutions under the original H l 's with a suitable choice of ∆ . This result is key in our construction and may find applications elsewhere.
Our main contributions are based on a number of technical lemmas and corollaries that are given in the Methods section and proven in detail in Supplementary Information.

B. Product formulas and effective operators
For a time-independent Hamiltonian H = L l=1 H l , where each H l is Hermitian, the evolution operator for time t is U (t) = e −itH . Product formulas provide a way of approximating U (t) as a product of exponentials, each being a short-time evolution under some H l . For p > 0 integer and s ∈ R, a p-th order product formula is a unitary where each s j ∈ R is proportional to s and 1 ≤ l j ≤ L. The number of terms in the product may depend on p and L, and we assume L = O(1), q = O(1).
(The more general case is analyzed in Supplementary Information.) We define |s| = q j=1 |s j | and also assume |s| = O(|s|). The p-th order product for- . One way to construct W p (s) is to apply a recursion in Refs. [18,19]. These are known as Trotter-Suzuki approximations and satisfy the above assumptions.
By breaking the time interval t into r steps of sufficiently small size s, product formulas can approximate U (t) as U (t) ≈ (W p (s)) r . We will refer to r as the Trotter number, and this number will define the complexity of product formulas that simulate U (t) within given accuracy. Note that the total number of terms in the product formula is actually qM r = O(M r), where M is the number of terms in the product decomposition of each e −isj H l j .
Known error bounds for product formulas that apply to the general case grow with h and can be large. However, error bounds for approximating the evolved state U (t) |ψ may be better under the additional assumption that |ψ is supported on a lowenergy subspace. We then analyze the case where the initial state satisfies Π ≤∆ |ψ = |ψ , where Π ≤Λ is the projector into the subspace spanned by eigenstates of H of energies (eigenvalues) at most Λ ≥ 0. We assume H l ≥ 0. Our results will be especially useful when ∆/h vanishes asymptotically, and ∆ will specify the low-energy subspace.
The notion of effective operators will be useful in our analysis. Given a Hermitian operator X and ∆ ≥ ∆, the corresponding effective operator is X = Π ≤∆ XΠ ≤∆ , which is also Hermitian. We also define the unitariesŪ (s) = e −isH andW p (s) by replacing the H l 's byH l 's in W p (s). Note that h = max l H l ≤ ∆ and U (t) |ψ =Ū (t) |ψ . Then, using the known error bound for product formulas, we obtain (U (s) −W p (s)) |ψ = O((∆ s) p+1 ). This error bound is a significant improvement over the general case if ∆ h, which may occur when ∆ h. However, product approximations of U (t) require that each term is an exponential of some H l , which is not the case inW p (s). We will address this issue and show that the improved error bound is indeed attained by W p (s) for a suitable ∆ .

C. Local Hamiltonians
We are interested in simulating the time evolution of a local N -spin system on a lattice. Each local interaction term in H is of strength bounded by J and involves, at most, k spins. We do not assume that these interactions are only within neighboring spins but define the degree d as the maximum number of local interaction terms that involve any spin. Next, we write H = L l=1 H l , where each H l is a sum of M local, commuting terms [32] and LM ≤ dN . Each e −isH l in a product formula can be decomposed as products of M evolutions under the local, commuting terms with no error.
These local Hamiltonians appear as important condensed matter systems, including gapped and critical spin chains, topologically ordered systems, and models with long-range interactions [33][34][35][36]. For example, for a spin chain with nearest neighbour interactions, L = 2 and each H l may refer to interaction terms associated with even and odd bonds, respectively. In this case, h = O(N ). We will present our results for the case k = O(1) and d = O(1) in the main text, which further imply L = O(1) [32]. Nevertheless, explicit dependencies of our results in k, d, L, and other parameters that specify H can be found in Supplementary Note 4. H l be a k-local Hamiltonian as above, H l ≥ 0, ∆ ≥ 0, 0 ≤ J|s| ≤ 1, and W p (s) a p-th order product formula as in Eq. (1). Then, where ∆ = ∆ + β 1 J log(β 2 /(J|s|)) + β 3 J 2 N |s| and the β i 's are positive constants, β 2 ≥ 1.
The proof of Thm. 1 is in Supplementary Note 3 and we provide more details about it in the next section, but the basic idea is as follows. There are two contributions to Eq. (2) in our analysis. One comes from approximating the evolution operator with a product formula that involves the effective Hamiltonians and, as long as ∆ ≥ ∆, this error is O((∆ |s|) p+1 ), as explained. The other comes from replacing such a product formula by the one with the actual Hamiltonians H l , i.e., W p (s). However, unlikē H l , the evolution under each H l allows for leakage or transitions from the low-energy subspace to the subspace of energies higher than ∆ . In Supplementary Information, we use a result on energy distributions in Ref.
[37] to show that this leakage can be bounded and decays exponentially with ∆ . Thus, this effective norm depends on ∆ and must also depend on s, as the support on high-energy states can increase as s increases, resulting in the linear contribution to ∆ in Thm. 1.
The log(β 2 /(J|s|)) factor in ∆ only becomes relevant when |s| 1. This term appears in our analysis due to the requirement that both contributions to Eq. (2) discussed above are of the same order. Thus, as s → 0, we require ∆ → ∞ to make the error due to leakage zero, which is unnecessary and unrealistic. This term plays a mild role when determining the final complexity of a product formula, as the goal will be to make s as large as possible for a target approximation error. It may be possible that this term disappears in a more refined analysis.
Let r = t/s be the Trotter number, i.e., the number of steps to approximate U (t) as (W p (s)) r .
Thus, for overall target error ε > 0, it will suffice to satisfy (U (s) − W p (s))Π ≤∆ = O(εs/t). This condition and Thm. 1 can be used to determine r as follows.
Each term of ∆ in Thm. 1 can be dominant depending on s and ∆. First, we consider the first two terms, and determine a condition in s to satisfy ((∆ + J)|s|) p+1 = O(εs/t), by omitting the log factor. Then, we consider another term and determine a condition in s to satisfy (J 2 N |s| 2 ) p+1 = O(εs/t). These two conditions alone can be satisfied with a Trotter number Last, we reconsider the second term with log, and we require (J log(1/(J|s|))|s|) p+1 = O(εs/t). As the first two conditions are satisfied with a value for s that is polynomial in N and tJ/ε, this last condition only sets a correction to the first term in r in Eq. (3) that is polylogarithmic in |t|J/ε. Thus, the overall complexity of the product formula for local Hamiltonians is given by Eq. (3), where we need to replace O byÕ to account for the last correction. Note that the number of terms in each W p (s) is constant under the assumptions and r is proportional to the total number of exponentials in (W p (s)) r .
We give a general result on the complexity of product formulas that provides r as a function of all parameters that specify H in Thm. 2 of Supplementary Note 4.
The error bounds for product formulas used in Thm. 1 depend on the norm of the effective Hamil-toniansH l . The assumption H l ≥ 0 will then assure that H l ≤ ∆ , which is sufficient to demonstrate the complexity improvements in Eq. (3).
In general, H l ≥ 0 can be met after a simple shifting H l → H l + a l , and the assumption seems irrelevant. However, this shifting could result in a value of ∆ (or ∆ ) that scales with some parameters such as the system size N . In this case, the error bound in Thm. 1 would be comparable to that of the worst case (without the low-energy assumption) and would not provide an advantage.
Nevertheless, for many important spin Hamiltonians, the assumption H l ≥ 0 is readily satisfied. The Heisenberg model of Fig. 1 More general (anisotropic) Heisenberg models as well as the so-called frustration-free Hamiltonians that are ubiquitous in many-body physics also satisfy the assumption [38,39], where our results directly apply. For this class of models, the ground state energy is zero. This class contains interesting low-lying states in the subspace where, e.g., ∆ = O(1).
We provide more details on potential complexity improvements for the general case (H l 0) at the end of Supplementary Note 3.

III. DISCUSSIONS
The best previous result for the complexity of product formulas (Trotter number) for local Hamiltonians of constant degree is O(τ 1+1/p N 1/p /ε 1/p ), with τ = |t|J [23]. Our result gives an improvement over this in various regimes. Note that, a general characteristic of our results is that they depend on ∆, which is specified by the initial state. Here we assume that ∆ is a constant independent of other parameters that specify H. The comparison for this case is in Table I. For p = 1, we obtain a strict improvement of order N 1/3 over the best-known result.
For higher values of p, the improvement appears for larger values of τ /ε that may scale with N , e.g., τ /ε =Ω(N p−2+1/(p+1) ). In Supplementary Note 5, we provide a more detailed comparison between our results and the best previous results for product formulas as a function of ∆ and other parameters that specify H.
A more recent method for Hamiltonian simulation uses a truncated Taylor series expansion of Here, r is the number of "segments", and U (t) is approximated as (U r ) r . A main advantage of this method is that, unlike product formulas, its complexity in terms of ε is logarithmic, a major advantage if precise computations are needed. The complexity of this method for the low-energy subspace of H can only be mildly improved. A small ∆ allows for a truncation value K that is smaller than that for the general case [7]. Nevertheless, the complexity of this method is dominated by r, which depends on a certain 1-norm H 1 of H that is independent of ∆. Furthermore, quantum signal processing, an approach for Hamiltonian simulation also based on certain polynomial approximations of U (t), was recently considered for simulation in the low-energy subspace [40]. While the low-energy constraint may also result in some mild (constant) improvement, the overall complexity of quantum signal processing also depends on H 1 . For local Hamiltonians where k, d = O(1), and for constant ∆, the overall complexity of these methods isÕ(τ N 2 ), where we disregarded logarithmic factors in τ , N , and 1/ε. Our results on product formulas provide an improvement over these methods in various regimes, e.g., when ε is constant.
The obtained complexities are an improvement as long as the energy ∆ of the initial state is sufficiently small. As we discussed, the assumption H l ≥ 0 was used and, while our results readily apply to a large class of spin models, it may be in conflict with ensuring small values of ∆ in some cases. It will be important to understand this in more detail (see Supplementary Note 3), which may be related to the fact that, for general Hamiltonians (H l 0), an improvement in the low-energy simulation could imply an improvement in the high-energy simulation by considering −H instead. Indeed, certain spin models possess a symmetry that connects the highenergy and low-energy subspaces via a simple transformation. Whether such "high-energy" simulation improvement is possible or not remains open. Additionally, known complexities of product formulas are polynomial in 1/ε. This is an issue if precise computations are required as in the case of quantum field theories or QED. Whether this complexity can be improved in terms of precision as in Refs. [6][7][8]41] is also open.
Our work is an initial attempt to this problem. We expect to motivate further studies on improved Hamiltonian simulation methods in this setting by refining our analyses, assuming other structures such as interactions that are geometrically local, or improving other simulation approaches.

A. Leakage to high-energy states
A key ingredient for Thm. 1 is a property of local spin systems, where the leakage to high-energy states due to the evolution under any H l can be bounded. Let Π >Λ be the projector into the subspace spanned by eigenstates of energies greater than Λ . Then, for a state |φ that satisfies Π ≤Λ |φ = |φ , we consider a question on the support of e −isH l |φ on states with energies greater than Λ . This question arises naturally in Hamiltonian complexity and beyond, and Lemma 1 below may be of independent interest. A generalization of this lemma will allow one to address the Hamiltonian simulation problem in the low-energy subspace beyond spin systems.
Lemma 1 (Leakage to high energies). Let H = L l=1 H l be a k-local Hamiltonian of constant degree as above, H l ≥ 0, and Λ ≥ Λ ≥ 0. Then, ∀ s ∈ R and ∀ l, where α 1 and α 2 are positive constants.
The proof is in Supplementary Note 1. It follows from a result in Ref.
[37] on the action of a local interaction term on a quantum state of low-energy, in combination with a series expansion of e −isH l . While the local interaction term could generate support on arbitrarily high-energy states, that support is suppressed by a factor that decays exponentially in Λ − Λ.
Another key ingredient for proving Thm. 1 is the ability to replace evolutions under the H l 's in a product formula by those under their effective low-energy versions (and vice versa) with bounded error. This is addressed by Lemma 2 below, which is a consequence of Lemma 1. The proof is in Supplementary Note 2, where we also provide tighter bounds that depend on ∆ .

B. Relevance to the main result
The consequences of these lemmas for Hamiltonian simulation are many-fold and we only sketch those that are relevant for Thm. 1. Consider any product formula of the form W = q j=1 e −isj H l j . Then, there exists a sequence of energies Λ q ≥ . . . ≥ Λ 0 = ∆ such that the action of W on the initial low-energy state |ψ can be well approximated by that of W Λ = q j=1 Π ≤Λj e −isj H l j on the same state. Furthermore, each Π ≤Λj e −isj H l j in W Λ can be replaced by Π ≤Λj e −isjH l j and later by e −isjH l j within the same error order, as long as Λ q ≤ ∆ .
In particular, for sufficiently small evolution times s j and ∆ h, the resulting effective norm satisfies ∆ h for local Hamiltonians. This is formalized by several corollaries in Supplementary Note 3. Starting from W , we can construct the product for-mulaW = q j=1 e −isjH l j . Lemmas 1 and 2 imply that both product formulas produce approximately the same state when acting on |ψ , for a suitable choice of ∆ as in Thm. 1. IfW is a product formula approximation ofŪ (s) = e −isH , it follows that U (s) |ψ =Ū (s) |ψ ≈W |ψ ≈ W |ψ .

V. DATA AVAILABILITY STATEMENT
All relevant data used for Fig. 1 are available from the authors.

VI. CODE AVAILABILITY STATEMENT
The code for the simulation results in Fig. 1 is available from the authors.
and Z i are the Pauli operators for the ith spin, and H = H 1 + H 2 , where H 1 and H 2 are the interaction terms represented by blue and red bonds, respectively. (a) The evolution operator for time s, U (s) = e −isH , is approximated by the first order product formula W 1 (s) = e −isH1 e −isH2 . The plot shows the largest approximation errors when acting on various low-energy subspaces associated with increasing energies, labeled by n = 1, 50, 150, 200, and in the worst case. (b) Similar results for when the evolution operator U (s) = e −isH is approximated by the second order product formula W 2 (s) = e −isH1/2 e −isH2 e −isH1/2 . Table I: Improvements of low-energy simulation: Comparison between the best-known complexity [23] and the complexity of low-energy simulation for p-th order product formulas. Results show the Trotter number for constant ∆ and local Hamiltonians on N spins with constant degree and strength bounded by J, and τ = |t|J. ε is the approximation error. TheÕ notation hides polylogarithmic factors in τ /ε.

Supplementary Information
In the following, we let H be a k-local Hamiltonian of N spins, where each interaction term involves at most k > 0 spins. We will write H = L l=1 H l , where each H l is a sum of at most M k-local and commuting interaction terms. The H l 's may be obtained via graph coloring [32], where a graph can be constructed with vertices that are labeled according to the subset of spins in each interaction term and the edges connect vertices associated with the same spins, but more efficient constructions may be possible. Indeed, in many interesting examples such as spins on the square lattice, H is already given in the desired form. The degree of H, i.e., the highest number of interaction terms that act non-trivially on any spin, is d > 0. The strength of each local interaction term is bounded by J > 0, hence H l ≤ JM and H ≤ JM L. Throughout this paper, . refers to the spectral norm (largest eigenvalue for positive semidefinite and Hermitian operators). The total number of local terms in H is then upper bounded by M L and dN . We will assume that H contains exactly M L terms and thus N ≤ M L ≤ dN with no loss of generality (e.g., we can add or subtract trivial terms to H l and each spin appears, at least, in one term). Furthermore, following the coloring procedure described above, we may assume L ≤ kd + 1 [32].

Supplementary Note 1: Proof of Lemma 1
We employ Theorem 2.1 in Ref.
[37] that, for an operator A, states Here, λ = 1/(2gk), where g is an upper bound on the sum of the strengths of the interactions associated with any spin. In our case, we take g = dJ and λ = 1/(2Jdk). If E A is the subset of local interaction terms in H that do not commute with A, R is the sum of the strengths of the terms in E A . For any H l , we note that (H l ) n is a sum of, at most, M n terms, each of strength bounded by J n and containing, at most, kn spins. For each such term, R ≤ Jdkn, and we obtain We now consider the Taylor series expansion of the exponential, The triangle inequality and Eq. (9) imply where α = eJ.
To prove the second result, we use Lemma 1 together with Eq. (25) and standard properties of the spectral norm, and obtain In particular, when k and d are constants, we have αM ≤ eJdN , and Eq. (33) implies where α 1 = 1/(2dk) and α 2 = ed are also positive constants.

Supplementary Note 3: Approximation errors for product formulas
We consider generic product formulas of q > 1 terms, where s = s 1 , . . . , s q , s j ∈ R, and 1 ≤ l j ≤ L. We also define where Λ = (Λ 1 , . . . , Λ q ), and Λ j ≥ 0 for all j. Using Lemmas 1 and 2, we now prove a number of results (corollaries) on the approximation errors for these product formulas that will be required for the proof of Thm. 1. First, we will prove that W (s) produces approximately the same state as W Λ (s) when the initial state is supported on the low-energy subspace and for a suitable choice of Λ. Next, we will show that the approximation error from replacing W Λ (s) byW Λ (s) is of the same order as that of the first approximation for a suitable choice of ∆ and effective Hamiltonians. A similar result is obtained if we further replacē W Λ (s) byW (s). Combining these results we show that the state produced by W (s) approximates that produced byW (s) for a suitable choice of ∆ .
In particular, ∆ ≥ Λ q = ∆ + 1 λ (α|s|M + q log(q/δ)). We use the identitȳ The triangle inequality and Corollaries 3, 2, and 1 imply To prove the previous corollaries, we often used (e α|sj |M − 1) ≤ e α|sj |M , but we note that a different upper bound such as (e α|sj |M − 1) ≤ α|s j |M e α|sj |M may provide improved results (e.g., a smaller value of ∆ in Cor. 4), especially if α|s j |M 1. However, while this other upper bound could be used to improve the constants hidden by the O notation in the complexity of our method, we were unable to improve its asymptotic scaling from using (e α|sj |M − 1) ≤ α|s j |M e α|sj |M .

Proof of Thm. 1
For some ∆ ≥ ∆ ≥ 0, letŪ (s) = e −isH be the evolution operator with the effective Hamiltonian and W p (s), p ≥ 1, be the corresponding p-th order product formula obtained by replacing H l →H l in W p (s) of Eq. 1. Since the condition H l ≥ 0 implies H l ≤ ∆ , we obtain where (∆ ) = γ(L∆ |s|) p+1 is an upper bound of the error induced by product formulas using effective operators and γ = O(1) is a constant [23]. This error bound grows with ∆ . It does not exploit any structure of the effective Hamiltonians so it may be possible to improve it under further constraints. Additionally, and then The other contribution to the error is due to Cor. 4, which can be turned around to obtain a bound on the error that depends on ∆ . Let λ = 1/(2Jdk), α = eJ, and q > 1 be the number of terms in the product W p (s). Then, Cor. 4 implies and This error bound decreases when ∆ increases. It is now valid for all ∆ ≥ 0 but can be irrelevant (larger than 1) when, for example, ∆ ≤ ∆. We assume that our product formula is such that |s| ≤ κL|s| for a constant κ ≥ 1 and let α = κα. Then Thus, for any ∆ ≥ ∆, the triangle inequality implies (U (s) − W p (s))Π ≤∆ ≤ (∆ ) + 5δ(∆ ). For given |s|, we can search for ∆ ≥ ∆ that minimizes the overall error bound. Let that ∆ be ∆ min . The global minimum for the error is then (∆ min ) + 5δ(∆ min ) and, for all ∆ ≥ ∆, it satisfies (∆ min ) + 5δ(∆ min ) ≤ (∆ ) + 5δ(∆ ) .
Then, we can fix any value of ∆ ≥ ∆ and obtain a bound for the overall error from computing (∆ ) + 5δ(∆ ). In particular, we choose and, assuming J|s| ≤ 1, q > 1, we obtain The constraint in J|s| is to avoid errors larger than 1: it is sufficient for δ(∆ ) ≤ 1 and for ∆ ≥ ∆. Therefore, if J|s| ≤ 1, where ∆ was determined in Eq. (71) andγ = γ + 5 is a constant. While our choice of ∆ in Eq. (71) does not correspond to the global minimum of the overall error, an exact calculation show that it is not far from ∆ min and provides the same asymptotic complexity for our method. Nevertheless, if one is interested in complexity of product formulas, we will analyze four different cases as follows. In the first case, we assume Here, H 1 is the 1-norm of H, given by L l=1 H l in our case, and H ind−1 is the so-called induced 1norm of H. The latter is defined as follows. We write H = N j1,...,j k =1 h j1,...,j k , where each h j1,...,j k includes the k-local interaction terms of qubits labeled as j 1 , . . . , j k in H. Then, H ind−1 := max l max j l N j1,...,j l−1 ,j l+1 ,...,j k =1 h j1,...,j k .
That is, we fix certain qubit j l and consider all the interaction terms that contain that qubit. For a degree d Hamiltonian with k-local interaction terms (not necessarily geometrically local), each of strength at most J, H ind−1 ≤ dJ. Furthermore, H 1 can be upper bounded as H 1 ≤ JM L ≤ JdN . As a result, for a k-local Hamiltonian as above, the best known upper bound for the Trotter number is To compare our main result with Eq. (115), we express Eq. (111) in terms of d and N . We recall that the total number of local terms in H is M L ≤ dN , L ≤ dk + 1, and we assume that k = O(1), M = O(N ), and q = O(L) = O(d) for the p-th order product formula. Then, in the asymptotic limit, The results for various values of p and k = O(1) are in Table III.  By setting a constraint on the initial state, our low-energy simulation result provides an advantage in certain regimes where d may grow with N . In the following, we will assume that, e.g., ∆ = O(d 2 J) and fix the value of t/ε, to simplify the expressions. Under these assumptions, for p = 1 and d = O(N 1/4 ), the terms in both columns of Table III are of order N 3/2 . For p = 2 and d = O(N 1/6 ), the terms in both columns of Table III are of order N 3/4 . For p = 3 and d = O(N 1/8 ) the terms in both columns of Table III are of order N 1/2 . For general p ≥ 1 and d = O(N 1 2(p+1) ), both complexities are comparable and of order N 3/(2p) . When d is fixed and N grows, our complexities may be worse than those obtained in Ref. [23]. One reason for this is because the effective norms may grow large in this case and the error bound that we use for product formulas using effective operators do not exploit any structure such as the locality of interactions. Better error bounds may be possible in this case, resulting in improved complexities. However, even if N is large, our results regain an advantage at certain values of t, in particular if we scale t/ε with N . Doing so will set a bound on the effective norm so that the second term in our complexities stops dominating.
The special case where k = O(1), d = O(1), and ∆ is also a constant independent of other parameters that specify H can be directly obtained from Table 2 and is given in Table 1 in the main text.