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.


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 shorttime 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 spatially local terms or when these terms satisfy Liealgebraic 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 2 × 6 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.

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 pth order product formula is OððΔ 0 sÞ pþ1 Þ, where s is a (short) time parameter and Δ 0 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 Δ 0 ( max l kH l k, 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 some of the complexity improvements in Table 1.
To obtain our results, we introduce the notion of effective Hamiltonians that are basically the H l 's restricted to act on a lowenergy subspace. The relevant norms of these effective operators is bounded by Δ 0 . 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 1 Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, USA. ✉ email: sahinoglu@lanl.gov; somma@lanl.gov nonlocal 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 Δ 0 . 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 "Methods" and proven in detail in Supplementary Information.

Product formulas and effective operators
For a time-independent Hamiltonian H ¼ P 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 2 R, a pth order product formula is a unitary W p ðsÞ ¼ e ÀisqH lq Á Á Á e Àis2H l 2 e Àis1H l 1 ; (1) where each s j 2 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 jsj ¼ P q j¼1 js j j and also assume jsj ¼ OðjsjÞ. The pth order product formula satisfies kUðsÞ À W p ðsÞk ¼ OððhjsjÞ pþ1 Þ, where h ¼ max l kH l k 4 . 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 qMr ¼ OðMrÞ, 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Þ ψ j i may be better under the additional assumption that ψ j i is supported on a low-energy 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 Δ 0 ! Δ, the corresponding effective operator is X ¼ Π Δ 0 XΠ Δ 0 , which is also Hermitian. We also define the unitaries UðsÞ ¼ e ÀisH and W p ðsÞ by replacing the H l 's by H l 's in W p (s). Note that h ¼ max l kH l k Δ 0 and UðtÞ ψ j i ¼ UðtÞ ψ j i. Then, using the known error bound for product formulas, we obtain kðUðsÞ À W p ðsÞÞ ψ j ik ¼ OððΔ 0 sÞ pþ1 Þ. This error  between the best-known complexity 23 and the complexity of lowenergy simulation for pth order product formulas.

Order
Previous result Low-energy simulation 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 τ/ε.
B. Şahinoğlu and R.D. Somma bound is a significant improvement over the general case if Δ 0 ( 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 in W p ðsÞ. We will address this issue and show that the improved error bound is indeed attained by W p (s) for a suitable Δ 0 .

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 ¼ P 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 neighbor 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. Table 2 summarizes the relevant parameters for the simulation of local Hamiltonians with product formulas.
Main result Theorem 1 Let H ¼ P L l¼1 H l be a k-local Hamiltonian as above, H l ≥ 0, Δ ≥ 0, 0 ≤ J|s| ≤ 1, and W p (s) a pth order product formula as in Eq. (1). Then, where Δ 0 ¼ Δ þ β 1 Jlog ðβ 2 =ðJjsjÞÞ þ β 3 J 2 Njsj 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 Δ 0 ! Δ, this error is OððΔ 0 jsjÞ 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, unlike H l , the evolution under each H l allows for leakage or transitions from the lowenergy subspace to the subspace of energies higher than Δ 0 . 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 Δ 0 . Thus, this effective norm depends on Δ and must also depend on s, as the support on highenergy states can increase as s increases, resulting in the linear contribution to Δ 0 in Thm. 1.
The log ðβ 2 =ðJjsjÞÞ factor in Δ 0 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 Δ 0 ! 1 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.
Each term of Δ 0 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ÞjsjÞ 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 Njsj 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 ðJlog ð1=ðJjsjÞÞjsjÞ 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 0 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 condition H l ≥ 0 The error bounds for product formulas used in Thm. 1 depend on the norm of the effective Hamiltonians H l . The assumption H l ≥ 0 will then assure that kH l k Δ 0 , 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 Δ 0 ) 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 is an example, where H l is a sum of terms like 1 − X i X j − Y i Y j − Z i Z j ≥ 0. 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.

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 1. 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 e ÀiHt=r % U r ¼ P K k¼0 ðÀiHt=rÞ k =k! 7 . 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 1norm ∥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 Oðτ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 high-energy and low-energy subspaces via a simple transformation. Whether such "high-energy" simulation improvement is possible or not remains open. In addition, 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.

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 Π > Λ 0 be the projector into the subspace spanned by eigenstates of energies greater than Λ 0 . Then, for a state ϕ j i that satisfies Π Λ ϕ j i ¼ ϕ j i, we consider a question on the support of e ÀisHl ϕ j i on states with energies greater than Λ 0 . 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 ¼ P L l¼1 H l be a k-local Hamiltonian of constant degree as above, H l ≥ 0, and Λ 0 ! Λ ! 0. Then, 8 s 2 R and ∀ l, kΠ > Λ 0 e ÀisHl Π Λ k e Àα1ðΛ 0 ÀΛÞ=J e α2JjsjN À 1 ; 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 ÀisHl . While the local interaction term could generate support on arbitrarily high-energy states, that support is suppressed by a factor that decays exponentially in Λ 0 À Λ.
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 Δ 0 . Lemma 2 Let H ¼ P L l¼1 H l be a k-local Hamiltonian of constant degree as above, H l ≥ 0, and Δ 0 ! Λ 0 ! Λ ! 0. Then, 8 s 2 R and ∀l, kΠ Λ 0 ðe ÀisHl À e ÀisHl ÞΠ Λ k e Àα1ðΛ 0 ÀΛÞ=J ðe α2JjsjN À 1Þ (5) and kΠ > Λ 0 e ÀisHl Π Λ k 3e Àα1ðΛ 0 ÀΛÞ=J ðe α2JjsjN À 1Þ ; where α 1 and α 2 are positive constants.

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

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

CODE AVAILABILITY
The code for the simulation results in Fig. 1 is available from the authors.