Complexity of Digital Quantum Simulation in the Low-Energy Subspace: Applications and a Lower Bound

Digital quantum simulation has broad applications in approximating unitary evolution of Hamiltonians. In practice, many simulation tasks for quantum systems focus on quantum states in the low-energy subspace instead of the entire Hilbert space. In this paper, we systematically investigate the complexity of digital quantum simulation based on product formulas in the low-energy subspace. We show that the simulation error depends on the effective low-energy norm of the Hamiltonian for a variety of digital quantum simulation algorithms and quantum systems, allowing improvements over the previous complexities for full unitary simulations even for imperfect state preparations due to thermalization. In particular, for simulating spin models in the low-energy subspace, we prove that randomized product formulas such as qDRIFT and random permutation require smaller Trotter numbers. Such improvement also persists in symmetry-protected digital quantum simulations. We prove a similar improvement in simulating the dynamics of power-law quantum interactions. We also provide a query lower bound for general digital quantum simulations in the low-energy subspace.


Introduction
Quantum computers possess the potential to simulate the dynamics of given quantum systems more efficiently compared to their classical counterparts [1][2][3].A main approach for simulating Hamiltonian evolution is the digital quantum simulation, which maps the quantum systems into qubits and approximates the evolution using digitalized quantum gates.The most standard digitalization approaches are product formulas, including the Suzuki-Trotter formulas [4][5][6][7][8][9][10], which decompose the target Hamiltonian into a sum of non-commuting terms and perform a product of exponentials of each individual term.Alternative principles have also led to advanced digital quantum simulation approaches such as linear combinations of unitaries [11,12], qubitization [13], truncated Taylor series [14], quantum signal processing [15,16], etc.Nevertheless, product formulas are still the most popular and promising approaches due to their relatively simple implementations, especially on near-term devices [17,18].
Recent works provided refined analysis for standard product formulas [19] and proposed improved digital quantum simulation algorithms including randomized compiling [20][21][22], random permutation [23], symmetry transformation [24], and multiproduct formula [12,25,26] based on product formulas.While most of these simulation algorithms are important and efficient for approximating Hamiltonian evolution, the resource requirement depends heavily on the norm of (the terms of) the Hamiltonian.However, a significant class of simulation tasks in quantum physics do not explore high-energy states, suggesting a potential to improve the resource demand under physically relevant assumptions on energy scales.As pointed out by Ref. [27], standard product formulas require smaller complexity to simulate the dynamics of k-local spin models when the states are restricted within a low-energy subspace.Hitherto, a general and complete theoretical framework for digital quantum simulations in the low-energy subspace still requires further exploration.In addition, Ref. [28] investigated eigenvector-dependent Trotter error which is closely related to simulating evolutions with quantum states chosen from a restricted subspace of the full Hilbert space, which is further explored by the concurrent work [29].Ref. [30] studied quantum simulation with low-dimensional initial state following spectral analysis.
To this end, we systematically investigate the Hamiltonian simulation problem in the low-energy subspace.We analyze the Trotter number and gate complexity of randomized and symmetry-protected digital quantum simulation and prove the significant improvements compared to the best-known results for simulating states chosen from the full Hilbert space.We show that similar improvement can be obtained for quantum models such as the power-law model.In general, we provide a lower bound for Hamiltonian simulation in the low-energy subspace with logarithmic dependence on ∆ and linear dependence on t.We further prove that the improvement above persists even when the quantum states are imperfectly prepared due to thermalization.Our results pave the way to more practical and efficient digital quantum simulations, which is an essential step toward the implementation of quantum simulation.After we posted the original version of our work, Ref. [31] proposed improved bounds of low-energy simulation.

Our results
We begin with some basic concepts for digital quantum simulations.Given a target Hamiltonian H = L l=1 H l consisting of L terms, we consider simulating an n-qubit quantum state ρ undergoing a unitary evolution V (t) = e −iHt for some time t using a unitary U .To achieve this goal, we split t into r (referred to as the Trotter number or the step complexity) identical Trotter steps such that the step size satisfies δ = t/r ≪ 1, and perform a product formula W (δ) = e −iδqH lq • • • e −iδ 1 H l 1 in each time step, where l 1 , . . ., l q are chosen from 1, . . ., L. A standard instance of product formula is the Lie-Trotter formula where u p = (4 − 41/(p−1) ) −1 [4][5][6].For large enough r, we have V (t) ≈ U (t) = S p (δ) r .The main challenge for implementing product formulas is to decide a proper r to guarantee that the simulation error is bounded below the error threshold ϵ.In the general case, the simulation error for a unitary simulation U is given by U ρU † − V ρV † .For technical simplicity, we consider the spectral norm ∥•∥ as the error metric in this paper.As U ρU † − V ρV † ≤ 2∥U − V ∥, we consider the spectral distance between the unitaries U and V in the following analysis.We remark that the trace norm is also usually employed to measure the distance in density matrices.To upper bound the distance between the output states of two unitary evolutions U ρU † and V ρV † , the diamond norm which is defined as the maximal trace distance between the output states of two unitary evolutions U(•) = U • U † and V(•) = V • V † given any identical pure state as input.It is further proved that the diamond norm between two unitary channels can be bounded by (see e.g.Lemma 3.4, [21]) Therefore, it is enough to consider the spectral norm ∥U − V ∥ in the following.
Given a Trotter number r, the Lie-Trotter formula provides a second-order simulation approximation error O(λ 2 H t 2 /r), where λ H = L l=1 ∥H l ∥ ≤ Lλ h and λ h = max l ∥H l ∥ [19].Therefore, choosing the Trotter number r 1 = O(t 2 λ 2 H /ϵ) is sufficient to ensure the error below ϵ.Similarly, choosing the Trotter number r p = O((λ H t) 1+1/p /ϵ 1/p ) is sufficient to ensure the error below ϵ for p-th order Lie-Suzuki-Trotter formulas.
It is sometimes possible to improve the complexity of product formulas with further information or assumptions in the form of Hamiltonians and states.In this paper, we consider digital quantum simulation in the low-energy subspace.Without loss of generality, we assume the Hamiltonian H is composed of L positive semi-definite terms H l ≥ 0. 1 For the simulation task, we only care about the low-energy quantum states that are below spectrum ∆ ≪ ∥H∥.We denote the projector to this subspace as Π ≤∆ .The formal definition for the low-energy simulation is given as follows: Definition 1.Given a Hamiltonian H with the above assumptions and a quantum state ρ with energy lower than ∆ ≤ ∥H∥, our goal is to find a simulation channel U(•) such that the simulation error is below some threshold ϵ, i.e., U(ρ) − V ρV † ≤ ϵ, where V = e −iHt is the target evolution of H with simulation time t.
In the special case when U(•) in Definition 1 is a unitary channel for some unitary U , we only need to consider the norm ∥(U − V )Π ≤∆ ∥ when estimating the upper bound for the simulation error between the simulation unitary U and the target evolution V .A specific well-studied model is the spin Hamiltonian H which has k-local interactions and strengths for each individual interaction term bounded below by J. Assuming each H l contains at most M interaction terms and the number of interaction terms acting on any spin bounded by d > 0, Ref. [27] showed that choosing a Trotter number ) for a p-th order Suzuki-Trotter formula suffices to achieve ϵ accuracy, where n is the number of qubits.Compared to the simulation complexity for Suzuki-Trotter formulas for the general case simulation, we can obtain an improvement on n dependence for p = 1 as λ H = O(n) for k-local Hamiltonians.Ref. [27] also showed that other improvements are possible for p ≥ 2 when choosing the parameters d, M , and J of the H within some regimes.This result indicates that standard product formulas applied to low-energy states have smaller errors than the worst-case error for k-local spin Hamiltonians.

Model
Approach Unitary simulation Low-energy simulation Table 1: A summary of low-energy simulation complexity (Trotter number) for product formulas (PF), qDRIFT, random permutation (rand.perm.), and symmetry protection (sym.prot.)approaches, with comparisons to the state-of-the-art results.The models include k-local Hamiltonians and power-law models.The results in this work are marked in the bold font.Here, Õ(•) omits the polylogarithmic dependence on the parameters.For k-local Hamiltonians, we assume that d, L, and J to be O(1) except for the random permutation approach.For the symmetry protection approach, 0 ≤ θ ≤ 1 is a constant depending on the structure of the Hamiltonian, see (21).In the last line, α is the decay factor for the power-law interactions.Here, g is the interaction strength and g = O(1), O(log n), and O(n 1−α/D ) for α > D, α = D, and α < D, where D is the dimension of the system.For random approaches, we only list the results considering the simulation errors in expectation.We derive the Trotter number required for ensuring the random fluctuation in the following sections.In addition, we provide corresponding gate complexities for the results in the above table.We can observe that for all the approaches, low-energy simulation improves power dependence on n, t, ϵ in a trade-off or decoupling way and n or spectral norm scaling is partially substituted by ∆.
In this work, we conduct a systematic study of quantum simulation under the lowenergy subspace assumption by proving the robustness of low-energy simulation under imperfect state preparation due to thermalization, proving the lower bound for low-energy simulations, analyzing the complexity required for implementing various digital quantum simulation approaches [20,23,24], and simulating different models such as k-local Hamiltonians and power-law models.We summarize our results in the following Table 1 with comparisons to the previous results.We also carry out extensive numerical experiments to benchmark the low-energy simulation using various algorithms in physically relevant models.

Robustness against imperfect state preparations due to thermalization
Concerning experimental implementations of quantum simulations, it is hard to guarantee that the support for the input states is restricted exactly in the subspace below the threshold ∆.In particular, the state preparation procedures on these near-term devices are imperfect, which produces a Gaussian tail on the higher energy states in the energy spectrum.Quantitatively, the probability distribution for eigenstates of energy E higher than ∆ satisfies a Gaussian distribution ∝ e −(E−∆) 2 /σ 2 , where σ ≪ ∆ is the variance [32][33][34].In this work, we prove that such imperfect state preparation causes an O((σ/∆) 2 ) multiplicative extra error in Section 2.1, which is ignorable compared to the original lowenergy simulation error.This result indicates that the low-energy improvement persists even under imperfect state preparations.

Randomized quantum compiling (qDRIFT)
. The first approach we analyze is the qDRIFT algorithm proposed by Ref. [20].We recall that the first-order Lie-Trotter algorithm deterministically cycles through every term H l in the Hamiltonian H in each step.The qDRIFT algorithm [20], however, approximates the target evolution using a channel constructed by averaging products Here, each unitary U i is a short-time evolution on a single term H l sampled according to the probability distribution H t 2 /ϵ) steps suffice to ensure the error within ϵ in expectation, which is explicitly independent of L.
In this work, we consider the performance of qDRIFT in the low-energy subspace.We denote H = Π ≤∆ ′ HΠ ≤∆ ′ for some carefully chosen ∆ ′ ≥ ∆ to be fixed later.We then introduce the parameter λ H = L l=1 H l where H l = Π ≤∆ ′ H l Π ≤∆ ′ .We now consider the low-energy version of qDRIFT where we evolve the system under H l λ H / H l in each time step with probability p l = H l /λ H . Assuming the Hamiltonian H is not ill-divided in the sense that λ H / H l = O(L) for all l = 1, . . ., L, we obtain the following complexity upper bound.

Theorem 1. Let H = L
l=1 H l be an n-qubit k-local Hamiltonian with parameters M , J, and d defined as above (LM J = O(n)).By choosing the Trotter number we can ensure that the expected simulation error in the low-energy subspace below ∆ for the qDRIFT algorithm is bounded by Moreover, if we choose the Trotter number then with probability at least 1 − χ, we can guarantee that the simulation error for a single product U is bounded by ϵ.
We leave the full proof for Theorem 1 in Section 3.1.Theorem 1 provides the step complexity (Trotter number) required for simulating low-energy states of k-local spin Hamiltonians both in the expectation case and in the case when one needs to ensure the performance of every single simulation unitary U with high confidence.Compared to the step complexity for implementing the qDRIFT algorithm in the full Hilbert space, the above theorem indicates that we can expect an improvement concerning the dependence on system size n as long as the number of terms in each H l M = Ω(1), which can be satisfied by most of the spin models.Although the step complexity (1) in expectation and the step complexity (2) to control random fluctuations are no better than that for the low-energy simulation using standard product formulas, we remark that this deficiency originates from the disadvantage of qDRIFT in simulating k-local spins even in the full Hilbert space [21].
Random permutation.Except for the qDRIFT algorithm, we consider the random permutation approach in simulating k-local Hamiltonians [23].While the qDRIFT approach implements random evolution on each term of the Hamiltonian, the random permutation approach randomly permutes the sequence of the terms in higher-order Suzuki-Trotter formulas.In particular, for any permutation in the permutation group σ ∈ S L , we denote where u p = (4 − 4 1/(p−1) ) −1 .In each time step δ, we consider averaging over all possible σ ∈ S L .For the simulation in low-energy subspace, we estimate the simulation error using In practice, it is difficult to implement all the permutations σ and average the results as the number of all permutations L! increases exponentially with L. We consider a random sampling procedure similar to the qDRIFT algorithm.We sample by implementing a random evolution p (δ) a randomly chosen permuted Suzuki-Trotter formula and σ i ∼ S L i.i.d.random permutations.In expectation, the evolution in each step is 1 L! σ∈S L S σ p (δ), which is the same as the standard random permutation approach.Concerning implementing this approach in the low-energy subspace, we obtain the following complexity result.
we can ensure that the expected simulation error in the low-energy subspace below ∆ for the random permutation approach is bounded by ϵ.For the random sampling implementation of the random permutation algorithm, if we choose the Trotter number then with probability at least 1 − χ, we can ensure that the simulation error for a single random sequence is bounded by ϵ.
We leave the technical details for Theorem 2 in Section 3.1.Theorem 2 shows the improvement over the step complexity of random permutation in the full Hilbert space [23] concerning the dependence on systems size n.Compared to the step complexity in (4) with that for standard product formulas in the low-energy subspace [27], we can obtain an L 1/2p reduction in the step complexity.This originates from the reduction from implementing the random permutation in the Hilbert space [23].
Based on the random permutation approach, Ref. [35] proposed an approach to double the order of approximation via the randomized product formulas.In particular, given r ≥ (5 p/2−1 + 1 2 )Lλ h t and p-th order contributing product formula, this algorithm provides a (2p + 1)-th order simulation for the ideal evolution.For this approach, we also prove the following corollary on the Trotter numbers required for low-energy simulations.
we can ensure that the expected simulation error in the low-energy subspace below ∆ for the doubling order approach is bounded by ϵ.And the Trotter complexity r prob to ensure the algorithm converges remains the same as that in Theorem 2.

Symmetry protection.
Other than the randomized digital quantum simulation approaches, another approach to improve the performance of digital quantum simulations is to implement a symmetry transformation in each time step [24].Given a Hamiltonian H, we assume that it is invariant under a group of unitary transformations denoted by C. For each unitary C chosen from C, we explicitly have [H, C] ≡ 0. The group C represents a symmetry group of the system.According to Ref. [24], we "rotate" each implementation of the Lie-Trotter formula S 1 (δ) via a symmetry transformation C µ ∈ C in each time step.The simulation for V = e −iHt reads V (t) ≈ r µ=1 C † µ S 1 (δ)C µ .Assuming the simulation error coherent, the digital quantum simulation of an evolution V (t) may end up with the time evolution of a different Hamiltonian H eff .We denote κ = H eff − H and rewrite the simulation under the assumption ∥κ∥ is small In the last line, we use the Baker-Campbell-Hausdorff (BCH) formula.The simulation error for this symmetry protection approach can thus be represented as For the symmetry protection approach, the Trotter number required for the k-local spin model is provided by Ref. [24] r = max where θ ∈ [0, 1] is a constant that depends only on the structure of the Hamiltonian H and the properties of the symmetry group C.An intuitive definition of θ is given by considering the simulation error v 0 of the Lie-Trotter formula.The scaling of the ratio between E Cµ∈C [C † µ v 0 C µ ] and ∥v 0 ∥ is n θ .Quantitatively, we have There are several examples provided with a specific value of θ.For example, when we draw symmetry transformations randomly, the behavior of the error would be analogous to a random walk, which results in θ = 0.5 under specific settings.The rigorous proof for this intuition for the localized Heisenberg model was provided in Ref. [24].Following a similar derivation in Ref. [26], one can explicitly obtain a construction where θ achieves the maximal value of 1. Now, we consider the performance of this symmetry protection approach in the lowenergy subspace, which is also an open question in Ref. [24].We estimate the simulation error by computing the following error term We prove the following theorem concerning the low-energy simulations of k-local Hamiltonians using the symmetry protection approach.

Theorem 3. Let H = L l=1 H l be an n-qubit k-local Hamiltonian with parameters M , J, and d (LM J = O(n)). By choosing the Trotter number
we can ensure that the symmetry protection approach provides an approximation within error ϵ in the low-energy subspace below ∆.Here, the value of 0 ≤ θ ≤ 1 is a constant depending on the structure of the Hamiltonian [24].
The proof for Theorem 3 is provided in Section 3.2.The above theorem provides the step complexity for the symmetry protection approach in the low-energy simulation.Compared with the complexity required for full Hilbert space simulations [24], the λ H = O(n) dependence is replaced by L∆ in two terms and O(n 1/2 ) in the other two terms, which enables an improvement concerning the n dependence.When θ > 0, the required Trotter number is smaller than that for implementing the standard first-order Lie-Trotter product formula in the low-energy subspace.
Power-law model.Theorem 1, Theorem 2, and Theorem 3 consider different digital quantum simulation approaches in simulating the low-energy quantum states of k-local spin models.We also consider power-law models, which can be regarded as a specific extension of k-local spin models at k = 2.In this case, there is no bound on the degree d compared to the k-local spin model, and additional constraints on the interaction strength J for each term.We consider a D-dimensional power-law interaction model consisting of n qubits.A power-law interaction H = i,j∈Σ H i,j with exponent satisfies: where i, j ∈ Σ are the qubit sites, Λ ⊆ R D is a D-dimensional square lattice, and H i,j is an operator supported on two sites i, j.Some notable examples include Coulomb interactions (α = 1), dipole-dipole interactions (α = 3), and van der Waals interactions (α = 6).According to Lemma H.1 of Ref. [19], we have Here, g is an upper bound of the strengths of the interactions associated with a single spin qubit.We can observe that every term of the power-law Hamiltonian is 2-local.In the following, we assume that each part H l contains at most M interaction terms.We obtain the following complexity result for simulating such a power-law model in the low-energy subspace.

Theorem 4. Consider simulating a D-dimensional power-law Hamiltonian H = L
l=1 H l with interaction strength g in the low-energy subspace ∆.Assume that each H l contains at most M interaction terms.Then, the Trotter number required to ensure the simulation error below ϵ is The corresponding gate complexity is given as: The proof for the above theorem is provided in Section 3.3.Compared with the best previously known complexity result for simulating the power-law model [19], the step complexity in (9) shows an improvement concerning the n dependence at p = 1.Notice that the best result for full Hilbert space simulations is obtained by simulating each term H i,j separately using one term H l .For p > 1, we can also obtain some advantage for specific choices of parameters M and ∆ in this case as we assume that H l 's consist of more than one term of H i,j .

Lower bound for low-energy simulations
We now consider the necessary number of queries to the Hamiltonian H to simulate any given state ρ in the low-energy subspace ≤ ∆.In Section 4, we show the following theorem as the lower bound on the quantum query complexity of simulations in the low-energy subspace:

Theorem 5 (Informal, see Theorem 6 for the formal version). There exists a Hamiltonian
H such that simulating some state with constant error within the low-energy threshold ≤ ∆ for some chosen scaled time τ requires at least Ω max{τ, log(∆) log log(∆) } queries to H.
The above theorem indicates that, similar to full space simulations [10], simulating Hamiltonians in the low-energy subspace also requires a linear number of queries to H in the simulation time τ .In addition, we also obtain a logarithmic dependence on the threshold ∆.Although this is far from the upper bound, where we require poly(∆) queries, we remark that this gap also exists in the full space simulation [36].The main idea to prove this theorem is to show a Hamiltonian of which exact low-energy subspace simulations for any time τ > 0 enable one to compute the parity of a string, which requires a linear number of queries to the length of the string [37,38].

Open questions
Our paper leaves several open questions for future investigations: • Ref. [39] developed a theory of average error for Hamiltonian simulation with random inputs and showed that improvement is possible when considering the average-case instead of the worst-case error.It is interesting to explore if we can obtain further improvement if we consider the simulation error for low-energy random states.
• In this paper, we study the simulation of evolution V on a low-energy quantum state ρ.It is also meaningful to study the complexity for estimating the expected value Tr(OV ρV † ) for some observable O instead of the full state V ρV † for a low-energy state ρ.
• Can we propose some Hamiltonians that are widely considered in the field of quantum information such that we can obtain improvement in the step complexity of digital quantum simulation under the low-energy assumption?
• Throughout the paper, the norm of the Hamiltonian is assumed to be bounded and its dimension is assumed to be finite.A natural open question is whether we can generalize our results concerning simulations in low-energy subspace to Hamiltonians with unbounded norm or infinite-dimensional Hilbert space.
Roadmap.The rest of the paper is organized as follows.In Section 2, we provide the general framework for low-energy simulations.In particular, we recap some of the results in [27] and prove several important lemmas for our results.We clarify the difference between our work and previous works.We further prove the robustness of the improvement originates from the low-energy assumption in imperfect state preparation due to thermalization.In Section 3, we consider the applications of low-energy settings for different digital quantum simulation approaches and quantum models.In Section 3.1, we consider randomized approaches including qDRIFT and random permutation for simulating k-local Hamiltonians, and provide the proofs for Theorem 1, Theorem 2, and Corollary 1.In Section 3.2, we consider the symmetry protection approach and provide the proof for Theorem 3. We show that the low-energy setting can provide improvements in power-law Hamiltonians in Section 3.3 and prove Theorem 4. Finally, we prove the query lower bound in Theorem 5 for simulating dynamics in the low-energy subspace in Section 4.

The General Framework
Notations.
We consider simulating the Hamiltonian H = L l=1 H l composed of n particles.We assume that each term H l contains at most M interaction terms and each interaction term acts on at most k qubits.For each qubit, we assume that the strength of interaction between this qubit and the rest qubits is bounded by g measured by the spectral norm.Without loss of generality, we assume that H l ≥ 0 for any l.We also assume the number of terms acting on each qubit is bounded by d and the strength of each term is bounded by J. Thus, we have g ≤ dJ.
Throughout this paper, we use the notation Õ(•) which omits the polylogarithmic dependence on the parameters.Given positive constant Λ, we denote Π ≤Λ the projector to the subspace spanned by eigenstates of H with energy smaller than Λ.We further denote Π >Λ = I − Π ≤Λ as the projector to the orthogonal subspace.We also summarize the notations used throughout this paper in Table 2.

Symbol Definition
Symbol Definition  Simulating k-local Hamiltonians in the low-energy subspace.For k-local Hamiltonians with bounded interaction strengths, the following lemma holds as the backbone of the framework for analyzing the simulation performance in the low-energy subspace.
Lemma 1 (Theorem 2.1 of [40]).Given H = L l=1 H l defined above with parameters g, L, and k, and any operator A, the following inequality holds where R A is the strengths of interaction terms acting on A and λ = (2gk) −1 .Here, Λ ′ ≥ Λ ≥ 0 are two positive values.
Based on Lemma 1, we can obtain some corollaries when A is an evolution of some terms H l .We list these lemmas in Appendix A. Using these lemmas, Ref. [27] obtained the following result concerning the performance of arbitrary product formulas of length q in the low-energy subspace.Given H = L l=1 H l an n-qubit k-local Hamiltonian with parameters M , J, and d (LM J = O(n)), the simulation error between the target evolution V (δ) and a product formula W (δ) of p-th order error is bounded by [27]: where α is some constant, and λ = (2Jdk) −1 .We can further deduce the following lemma regarding the Trotter number required for low-energy simulations based on the above error decomposition.
Lemma 2 (Eq.( 111) of [27]).Let H = L l=1 H l be an n-qubit k-local Hamiltonian with parameters M , J, and d (LM J = O(n)).By choosing the Trotter number we can ensure that the p-th order product formula W (δ) provides an approximation within error ϵ in the low-energy subspace below ∆.
Furthermore, as the implementation of the p-th order Suzuki-Trotter formula requires O(2 • 5 p/2−1 L) gates [6,10], the gate complexity required is Notice that LM J = O(n), the step complexity has an explicit n 1 2 1+ 1 2p+1 term.For p = 1, the step complexity obtained by Lemma 2 achieves a reduction concerning the dependence on n compared to the performance for simulating high energy states.For general p ≥ 1 and d = O(n 1 2(p+1) ), we can also obtained a similar improvement.As an illustration for the reduction of simulating dynamics of Hamiltonian assuming input states from the low-energy subspace, we perform numerical experiments to benchmark the product formulas in simulating k-local Hamiltonians.In particular, we study the following homogeneous Heisenberg model without external fields: where X i , Y i , and Z i are Pauli operators acting on the i-th spin, and the summation is over every adjacent spin pair.Since the ground energy E 0l < 0 for every H l , we implement the shift H l → H l − E 0l I for every H l .We empirically pick the threshold ∆ = 4 to make sure that the corresponding low-energy subspace is neither (close to) the full Hilbert space nor an empty subspace.
As shown in the left-hand side of Figure 1, the Hamiltonian in (13) is decomposed into three interaction terms: H = H 1 + H 2 + H 3 , respectively represented by yellow, green, and red bonds.Here, a bond represents all the interactions between a pair of spins.In the right-hand side of Figure 1, we plot the simulation error ϵ as a function of single-step evolution time δ with respect to different Trotter numbers r.The time evolution operator V (δ) is approximated by the second-order product formula S 2 (δ).We can observe that the projection to the low-energy subspace significantly reduces the error compared to the full Hilbert space for the same parameter configuration.We further remark that the scaling of the numerical simulation regarding parameters fits perfectly with the theoretical results.

Low-energy simulation with imperfect state preparations due to thermalization
In the previous analysis, we focused on simulating the dynamics of low-energy input states for a given Hamiltonian.We assumed the quantum state is chosen in the low-energy subspace under a threshold ∆.However, the state preparation on near-term devices is not perfect, resulting in a Gaussian tail distribution over higher energy states on the energy spectrum.It is natural to ask if the improvement obtained by low-energy simulations persists under the imperfect state preparation scenario.In the following, we provide an affirmative answer to this question.
In real quantum systems and experiments, the actually prepared states usually deviate from the ideal states due to thermalization [32][33][34], which results in a Gaussian tail on the spectrum distribution beyond the energy threshold ∆.We focus on pure states.and suppose the perfect input state |ψ p ⟩ has energy ∆, and denote the actual input state to be |ψ⟩.After thermalization deviation, the mean energy and variance of H in the state |ψ⟩ is In the case when H is local and |ψ⟩ has finite correlation length, σ 2 ψ will scale as O(n) [32].Assuming imperfect state preparation, we adapt E ψ = ∆ and σ 2 ψ ≪ ∆ from current analog state preparation schemes [33,34].Quantitatively, the support of the imperfect input state is a Gaussian distribution with parameter σ ψ around ∆.By using the error representation in (10) and (11), we can compute the simulation error Err imp as: where Recall that for standard Suzuki-Trotter formulas obtained by the Yoshida method [41], p is an even number.Therefore, we can compute Err imp as As σ ψ ≪ ∆ ≤ ∆ ′ and p a constant, we can discard the lower-order terms for ∆ ′ and obtain We compare the above error decomposition with (10) and observe that we only get an additional O(∆ ′ (p−1) σ 2 ψ ) term when the state preparation is imperfect.Therefore, we conclude that for imperfect state preparation due to thermalization, the improvement obtained by the low-energy simulation assumption persists.We further claim that this robustness persists for all the analyses on any following low-energy performance for various digital quantum simulation approaches.This is because the low-energy simulation error for any following algorithms has a polynomial scaling on ∆ ′ .We can thus employ similar calculations as (14) and prove the robustness.

Applications
Intuitively, the resource requirement for different simulations reduces when the evolved quantum system does not consider high-energy states.Yet, Lemma 2 only provides strict proof for standard product formulas and k-local Hamiltonians.In this work, we explore other digital quantum simulation algorithms and Hamiltonians.In this section, we provide the technical details for our results in applying low-energy analysis for different digital quantum simulation algorithms and Hamiltonians.

Randomized product formulas
In this subsection, we analyze the performance of randomized simulation approaches in the low-energy subspace.In particular, we provide the proof for Theorem 1, Theorem 2, and Corollary 1.

qDRIFT algorithm.
We begin with recapping the qDRIFT algorithm in the below Algorithm 1.

Algorithm 1: The qDRIFT algorithm
evolution time t, and number of steps r 1: At each time interval t/r: evolve a random term in Hamiltonian U i = e −i(t/r)X i , where X i is randomly chosen as λ H H l / H l with probability H l /λ H Output: The unstructured (randomly generated) product formula U = U r . . .U 1 We now analyze the performance of the qDRIFT algorithm in the low-energy subspace and prove Theorem 1.We still consider the k-local Hamiltonians and parameters d, M , and J defined above.We carefully pick ∆ ′ ≥ ∆ a value to be determined later.We denote Given a random sequence U r U r−1 . . .U 1 generated by Algorithm 1, our goal is to estimate the error where V = e −iHt is the ideal evolution.To begin with, we decompose this error into three terms: where V = e −iHt and U i are obtained by running Algorithm 1 on H with corresponding H l 's, and choosing X i = λ H H l / H l in the i-th step with probability H l /λ H . Here, the first term is the projection error for the low-energy subspace estimation.The second term is the random fluctuation in the low-energy subspace.The third term is the deterministic bias in the low-energy subspace.
We first consider the deterministic bias error term.For the projected Hamiltonian H, λ H = L l=1 H l = O(L∆ ′ ) for some ∆ ′ ≥ ∆ to be fixed later.It is worthwhile to mention that for k-local Hamiltonians, λ H is explicitly independent of system size n while λ H = O(n) [19].According to Lemma 3.5 of Ref. [21], for Hamiltonian H = L l=1 H l and U i , the deterministic bias term is bounded by Here The next step is to provide a bound for the random fluctuation term U r . . .U 1 − (E[U i ]) r .For the randomized sequence generated by the qDRIFT algorithm projected into the lowenergy subspace U r . . .U 1 , we denote We can verify that this sequence {B k : k = 0, . . ., r} satisfies the martingale properties as follows: • Causality: every B k completely depends on the previous information of B k−1 , . . ., B 1 , i.e., the choice of U k−1 , . . ., U 1 .
• Status quo: the expectation value for B k equals to B k−1 conditioned on the previous history, i.e., The property for such martingale features similar properties to an unbiased random walk.Based on Freedman's inequality and its application to martingales [42][43][44][45][46][47], we have Lemma 10 in Appendix A. Starting from this lemma, we consider the martingale provided by the qDRIFT algorithm defined above.We define where the second inequality follows the convexity of mathematical expectation to arrive at the fourth line.
We invoke Lemma 10 with v bounded by r∥C k ∥ 2 and obtain that , where 2 n is the dimension of the Hilbert space.With probability at least 1 − χ, we have Finally, we provide the bound for the projection error (U r . . .U 1 − U r . . .U 1 )Π ≤∆ .Intuitively, we have This is because [U i , Π ≤∆ ′ ] = 0 for any i = 1, . . ., r and Π ≤∆ ′ Π ≤∆ ′ = Π ≤∆ ′ .However, this yields no bound on this term (U r . . .U 1 − U r . . .U 1 )Π ≤∆ .To address this issue, we consider the following two-phase decomposition of a qDRIFT algorithm.We write r = r 1 r 2 with r 1 and r 2 to be fixed later.We decompose the qDRIFT sequence of r steps into r 1 groups of "random product formulas", each of r 2 steps.The full qDRIFT sequence can be written as W r 1 . . .W 1 where W i = U i,r 2 . . .U i,1 for i = 1, . . ., r 1 with U i,j the random qDRIFT step for j = 1, . . ., r 2 .We thus decompose the simulation error as where V = e −iHt the second line follows from the fact that We then consider an arbitrary W i (denoted as W in the following), and employ the above error decomposition as: According to the above calculation, with probability 1 − χ we have the second and the third terms as: where λ H = O(L∆ ′ ).
We now provide a bound for the truncation error (U r 2 . . .U 1 − U r 2 . . .U 1 )Π ≤∆ .We assume that U i evolves on H i for some i = 1, . . ., r 2 and denote δ i , the error can be bounded as follows follows from Lemma 7 and Lemma 8: where α = eJ, λ = (2Jdk) −1 , and δ ′ = (δ ′ r 2 , . . ., δ ′ 1 ) with |δ ′ | = O(Lt/r 1 ).Now, we are ready to prove Theorem 1.We first consider the expectation error and prove (1).In this case, we do not need to guarantee the random fluctuation term.We compute the r 2 such that the projection error is of the same scale as the deterministic bias, i.e., Based on the previous bound on these two terms of error, we conclude that To achieve this, we can fix the value of ∆ ′ as follows:

The total error is thus
with λ H = O(L∆ ′ ).Now, we consider three cases separately to derive the final step complexity.
In the first case, ∆ is the dominant term of ∆ ′ .In this case, we require In the second case, O 1 λ α t r 1 M L is the dominant term of ∆ ′ .We require In the last case, is the dominant term of ∆ ′ .Except for a mild polylogarithmic correction, this term is the same as r 2 log r 2 /λ.In this case, we can prove that the corresponding step complexity is not the dominant term.Optimizing over all choices of r 2 , we have the step complexity minimized when r 2 = O(1).
In general, the Trotter step complexity is which finishes the proof for (1).The corresponding gate complexity is given that O(L) gates are required to implement a Suzuki-Trotter formula in each time step.Now, we consider the step complexity required to ensure qDRIFT converges and prove (2).If we want to further control the random fluctuation, at r 2 = O(1) we require To achieve this, we can fix the value of ∆ ′ as follows: . Following a similar procedure to the average-case performance, we derive the Trotter number to ensure a 1 − χ success probability for qDRIFT as which finishes the proof for Theorem 1.In addition, the corresponding gate complexity is .

Random permutation.
In the previous qDRIFT algorithm, we only consider applying a random evolution in each time step, which is a random version of the first-order Lie-Trotter formula.A randomized approach for higher-order formulas, known as random permutation, was proposed in Ref. [23].As shown in Algorithm 2, this algorithm randomly permutes the order of Hamiltonian terms within each block to S σ i p (δ = t/r) for a permutation σ i ∈ S L from the permutation group.

Algorithm 2: The random permutation algorithm
Input: Hamiltonian H = L l=1 H l , evolution time t, and number of steps r 1: At each time interval δ = t/r: evolve a random term in Hamiltonian U i = S σ i p (δ) with a randomly chosen sequence with probability 1/L! Output: The unstructured (randomly generated) product formula U = U r . . .U 1 We now consider the performance of this algorithm in the low-energy subspace and prove Theorem 2. Similar to the previous section, we can decompose the error into three terms: where U i = S σ i p (δ) with Hamiltonian H and U = e −iHt .The first term is the projection error for the low-energy subspace estimation.The second term is the random fluctuation in the low-energy subspace.The third term is the deterministic bias in the low-energy subspace.We introduce two new parameters We begin with the deterministic bias error term.According to the direct calculation [4,23], we can derive that (see, e.g.Eq. (C14) of Ref. [21]) According to Lemma 3.5 of Ref. [21], for Hamiltonian H = L l=1 H l and U i , the deterministic bias term is bounded by Here, U 1/r = e −iHt/r = e iHδ .
The next step is to provide a bound for the random fluctuation term U r . . .
We still use the Freedman inequality for martingales, and define For the random permutation approach, each C k is bounded by Here, the last equation follows from the fact [5,23] that L p .Therefore, we can derive from Lemma 10 that with probability at least 1 − χ the random fluctuation scales as Finally, we provide the bound for the projection error (U r . . .U 1 − U r . . .U 1 )Π ≤∆ .According to Lemma 9, we decompose the projection error as follows: where λ = (2Jdk) −1 , α = eJ, and q(p) is the length of the product formula in each step depending only on p. Now, we are ready to derive Theorem 2. Similar to the previous section, we still consider two cases: In the first case, we consider the average-case performance of the random permutation approach in the low-energy subspace and prove (4).In this case, we ignore the random fluctuation and compute the step complexity r such that the projection error is of the same scale as the deterministic bias, i.e., We can derive the bound on ∆ ′ and r as This finishes the proof for (4), which provides an O(L 1/p ) or O(L 1/(2p+1) ) reduction on the step complexity compared to standard product formulas [27].This means that as p increases, the random permutation will eventually lose its reduction-a similar observation to that in Ref. [23].The gate complexity requires an O(L) overhead from the Trotter number r.
In the second case, we want to further control the random fluctuation.Thus we require We can derive the bound on ∆ ′ and r to ensure that the simulation error is bounded by ϵ with probability at least 1 − χ as which finishes the proof for (5).The overhead for controlling randomized fluctuation will quickly disappear as the order p increases.Yet, the advantage of this randomized approach compared to deterministic product formulas would also disappear correspondingly.
Similarly, The gate complexity requires an O(L) overhead from the Trotter number r.

Doubling the order of product formula.
In Ref. [35], Cho, Berry, and Huang proposed an approach to double the order of digital quantum simulations via the randomized product formulas.We also consider the performance of this approach in the low-energy subspace and obtain Corollary 1.We refer to Ref. [35] for a detailed description of this algorithm.As the algorithm requires the implementation of some random well-designed unitaries according to some well-designed probability {p h,l }, we have to assume that U (l) h 's also satisfy the properties of exp(−iH l δ)'s.We still consider the error decomposition as follows: According to Theorem 2 of Ref. [35] and the assumption that r ≥ (5 p/2−1 + 1 2 )Lλ h t, we can first compute the deterministic bias as For the random fluctuation, we still define Therefore, we can derive that with probability at least 1 − χ, the random fluctuation scales as For the projection error, we have where λ = (2Jdk) −1 , α = eJ, and q(p) is the length in each step depending only on p.
To prove Corollary 1, we consider the step complexity required to ensure the averagecase performance and the convergence of the algorithm, respectively.In the first setting, we ignore the random fluctuation and compute the step complexity r such that the projection error is of the same scale as the deterministic bias, i.e., We can derive the bound on ∆ ′ and r as This finishes the proof for r exp in (6).For r prob , we follow a similar proof to Theorem 2 to obtain the same result.The gate complexity required correspondingly is To benchmark the performance of random Hamiltonian simulation approaches in the low-energy subspace, we further carry out extensive numerical experiments concerning simulating the one-dimensional localized homogeneous Heisenberg models of different system sizes n ∈ [4,8] without external fields as follows: where X i , Y i , and Z i are Pauli operators acting on the i-th spin, and the summation is over every adjacent spin pair.The Hamiltonian in (20) is decomposed into three interaction terms:  The solid lines denote errors in the full Hilbert space while the dashed lines denote errors in the low-energy subspace.We distinguish the total evolution time t or the order of the product formula p by the color of the lines.We also compare the numerical and the theoretical error scalings.(a) We plot the worst-case and the low-energy simulation errors of the qDRIFT algorithm as a function of step number r.(b) The same as the upper left subfigure, but for the random permutation approach at p = 2 and p = 4. (c) We plot the worst-case and the low-energy simulation errors of the qDRIFT algorithm as a function of system size n.(d) The same as the lower left subfigure, but for the random permutation approach at p = 2 and p = 4.
by the curly brackets.Since the ground energy E 0l < 0 for every H l , we implement the shift H l → H l − E 0l I for every H l .We empirically pick ∆ = 14 (the shifted ground energy of the largest system n = 8) such that for system sizes n ∈ [4,8], the low-energy subspace covers neither (nearly) the whole spectrum nor a vacant subspace.We also provide the theoretical scaling for simulation error in the figures.We remark that as we can not theoretically calculate the constant coefficient, we cannot compare the exact values of theoretical and numerical simulation error.In the remaining figures, we compare the scaling between these two cases instead.
In Figure 2, we plot the qDRIFT error and the random permutation error as a function of step number r and system size n in four subfigures.We fix the total evolution time t and randomly sample 10 sequences for each choice of parameters.For the qDRIFT approach, in order to distinguish between closely spaced error bars, we made a tiny shift to the abscissa of nearby data points.We can observe that the projection to the low-energy subspace reduces the error compared to the full Hilbert space for the same parameter configuration.We plot the numerical results for the qDRIFT approach in the left two figures of Figure 2. We estimate that ϵ = O(nr −1/2 ) for full Hilbert space simulations and ϵ = O(n 2/3 r −1/2 ) for low-energy simulations.The numerical performances for both cases have better scaling concerning the n dependence than theoretical bounds as the latter one gives ϵ = O(n 3/2 r −1/2 ) and ϵ = O(max{n 1/2 r −1/2 , n 3/2 r −3/2 }) for the two cases as shown in Ref. [21] and (16).We remark that we take the random fluctuation into consideration here as the variance is much larger than the average value for simulation errors.In the right two figures, we plot the numerical results for the random permutation approach.We estimate that ϵ = O(n p/2+1 r −p ) and ϵ = O(nr −p ) for worst-case and lowenergy simulations, which are again better than the theoretical bounds of ϵ = O(n p+1 r −p ) and ϵ = O(max{n 1/2 r −p−1/2 , n p+3/2 r −2p−3/2 }) in Ref. [23] and (19).We can conclude that numerical error for the low-energy subspace has a better performance than the theoretical bounds in numerics and provides a steady reduction in the simulation error with various parameters and settings.

Symmetry protections
In this subsection, we consider an alternative approach, symmetry protection [24].Given a Hamiltonian H, we consider the group of unitary transformations denoted by C such that [H, C] ≡ 0 for each unitary C chosen from C. As mentioned in Section 1.1, we implement the Lie-Trotter formula S 1 (δ) in each time step with a symmetry transformation C ∈ C. Explicitly, the simulation for V = e −iHt reads where V (t) = e −iHt is the ideal evolution.We briefly recap the results in Ref. [24] for the intuition of reducing simulation error using symmetry protection.Intuitively, U (t) can be represented as an evolution on a slightly different Hamiltonian H eff and where κ = H eff − H.
For k-local Hamiltonians, we define the following quantities [19]: Ref. [24] obtains the following result concerning the simulation error in each time step: Lemma 3 (Lemma 7 of [24]).For all δ such that βδ ≤ 2α and α 2 δ ≤ β + γ, the Trotter error in each step is represented by where and V(δ) ≤ ιδ 3 with ι = 5(γ + β)/6.Based on the above result, the Trotterization error for the full Hilbert space simulation is: Lemma 4 (Theorem 1 of [24]).The simulation error for the whole evolution using symmetry protection for k-local Hamiltonians is bounded by: where is the averaged first-order error term.
We then focus on v 0 .According to Ref. [24], in general ∥v 0 ∥ ∝ ∥v 0 ∥/r θ for θ ∈ [0, 1].Hence, it suffices to choose a step complexity of to bound the simulation error below ϵ.There are several examples given for the value this θ.For example, when we draw symmetry transformations randomly, the behavior of the error would be analogous to a random walk, which results in θ = 0.5.The rigorous proof for this intuition for localized Heisenberg model was provided in [24].There are also some instances when θ = 1 [24,26].Now, we consider the performance of this symmetry protection approach in the lowenergy subspace, which is also an open question in [24].Our goal is to compute the following error term We start with the following observation: Proposition 1. Suppose S is the symmetry group for H, then S is a subgroup of the symmetry group for H, where H = Π ≤∆ ′ HΠ ≤∆ ′ .
We then decompose the error into three parts: where S 1 (δ) is the Lie-Trotter formula for H = L l=1 H l .The first part is 0 according to the property of the projector.The second part, according to the results of the symmetry protection approach, gives: According to Lemma 9 and the fact that q = L for first-order Trotterizations, we have Analogous to the intuition from symmetry transformation, we assume for some θ ∈ [0, 1].Here, we remark that the simulation error for the symmetry-protection approach satisfies θ ∈ [0, 1] is a rigorously proved result from Ref. [24] as shown in (7).
We can derive the step complexity here as (λ = (2Jdk) −1 ): which finishes the proof for Theorem 3. The gate complexity of this approach depends on the property of each symmetry transformation C µ and varies from case to case.However, we can still observe that the step complexity above is strictly better than that provided by either standard Trotterizations or the symmetry-protected simulation in the full Hilbert space.
We conduct numerical experiments concerning simulating the one-dimensional Heisenberg models of different system sizes n ∈ [4,8] in (20).We consider three different schemes of simulations: the standard first-order Trotterization ("Standard"), the first-order Trotterization protected by a random set symmetry transformation ("Random-ST"), and the first-order Trotterization protected by the optimal set ("Optimal-SP").For the "Random-ST" approach, considering the SU (2)-symmetry of Heisenberg models, each step we sample a 2 × 2 matrix from norm distribution ensuring diagonal conjugation, then scale the determinant to 1 and finally tensorize the normalized matrix to the n-th power.We randomly repeat this process 10 times and plot the average value as well as error bars.For the   "Optimal-SP" approach, we employ C 0 = U ⊗n H , C µ = C µ 0 as [24], where U H denotes singlequbit Hadamard gate.
In Figure 3, we plot the symmetry protection error as a function of step number r and system size n in two subfigures.In our numerical simulations, we fix the total evolution time t = 1.0.We can observe that the projection to the low-energy subspace reduces the error compared to the full Hilbert space for the same parameter configuration.We estimate that full Hilbert space simulations and low-energy simulations respectively yield ϵ = O(nr −1 ) and ϵ = O(n 2/5 r −1 ) for the standard Trotter formula, ϵ = O(nr −3/2 ) and ϵ = O(n 1/5 r −3/2 ) for the random set symmetry transformation, ϵ = O(nr −2 ) and ϵ = O(n 4/5 r −2 ) for the optimal set symmetry protection.The numerical performances for all the three approaches conform to the scaling of theoretical bounds, which give ϵ = O(max{nr −(1+θ) , n 2 r −2 }) and ϵ = O(max{r −(1+θ) , r −2 , nr −(2+θ) , n 2 r −4 }) for the full Hilbert space and low-energy subspace simulations.

Simulating power-law models
Except for k-local Hamiltonians, we also show that simulating low-energy states provides improvement for power-law interactions H = i,j∈Σ H i,j with exponent satisfying: where i, j ∈ Σ are the qubit sites, Λ ⊆ R D is a D-dimensional square lattice, and H i,j is an operator supported on two sites i, j.Through direct computation [19], we have Here, g is an upper bound on the strengths of the interactions associated with a single spin.We can observe that every term of the power-law Hamiltonian is 2-local.As assumed by the Theorem 4, each part H l contains at most M k-local interaction terms.We consider an operator A, and denote E A as the subset of interaction terms in H that do not commute with A and R A as the strengths of the terms in E A .Then we can bound R by R ≤ 2g and R ≤ 2gn for any H l or H n l for arbitrary n.We consider the following term: Using the Hadamard formula [48], we write For ℓ = 0, we have ∥K ℓ ∥ = ∥A∥, and for ℓ > 0, we have where h X denotes a term acting a two-qubit string X and the sum X j |X j−1 ,...,X 1 denotes a summation over the non-zero terms in the commutator.After shifting each h X to positive semi-definite operator, we denote hX [40].Thus the norm of ∥K ℓ ∥ can be bounded by Now, we compute the value for each term.Notice that the sum of the norms of h X that do not commute with A is bounded by R A .The sum of the norms of h X that does not commute with another h Y is bounded by 2g.We have Therefore, we have where m := R A 2g .Thus, we have Formally, we can prove the following lemma: Lemma 5. Consider a Hamiltonian with 2-local terms and suppose the interaction on each qubit is bounded by strength g.For any 0 ≤ s ≤ 1 2g , we have e sH Ae −sH ≤ ∥A∥(1 − 2gs) −R A /2g , where g, R A , and E A are defined as above.
Following the relationship [40] that we obtain the following lemma: Lemma 6.Consider a Hamiltonian with 2-local terms and suppose the interaction on each qubit is bounded by strength g.We have where λ := (4g) −1 .
According to the above lemma, we can deduce that where we assume that Λ ′ − Λ − 2R ≥ 0.
Thus, we can deduce that for α = eJ a constant.Compared with Lemma 8 and Lemma 7 in Appendix A, we find that k is exactly replaced by a constant 2g.Thus, similar to the proof for Lemma 2, we can prove the step complexity in Theorem 4 as: It is straightforward to observe that the gate complexity required is Compared with the gate count required for full-energy simulation tasks in [19], our result gets rid of the explicit overhead of O(n 2 ) gates to implement each Trotter step.When we combine multiple terms in each H l , the overhead can be largely reduced.
Here, we perform numerical experiments to benchmark the product formulas in simulating low-energy states for power-law Hamiltonians.In particular, we study the following 2D power-law Heisenberg spin-1/2 model without external fields (as shown in the left-hand side of Figure 4): where X i , Y i , and Z i are Pauli operators acting on the i-th spin, and the summation is over every spin pair.Since the ground energy E 0l < 0 for every H l , we make a shift H l → H l − E 0l I for every term.We empirically pick ∆ = 15 such that the low-energy subspace contains eigenstates suitable.We plot the Trotter error as a function of single-step evolution time at different decay factors α in the right-hand side of Figure 4.The Hamiltonian in (22) is decomposed into three interaction terms: the horizontal, the vertical, and the remaining, respectively represented by yellow, green, and red arrows as shown in the left-hand side of Figure 4.The target evolution within each time step V (δ) is approximated by the second-order product formula S 2 (δ).We can observe that the projection to the low-energy subspace significantly reduces the error compared to the full Hilbert space for the same parameter configuration.Moreover, in the full Hilbert space, the error increases as α grows, while in the low-energy subspace, the opposite is true.This could be explained by the fact that the norm of the Hamiltonian decreases as α increases.Therefore, ∆ = 15 becomes less stringent compared to the whole spectrum of the Hamiltonian as α increases.For a fixed ∆, the proportion of the spectrum in the low-energy subspace increases for larger α, which results in a smaller improvement compared to the full Hilbert space simulation.We further remark that the scaling of the numerical simulation regarding parameters fits perfectly with the theoretical results.

A Lower Bound for Simulating in the Low-Energy States
In the above discussion, we show that Hamiltonian simulation in the low-energy subspace can provide improvement compared to full space simulation.It is natural to ask what is the lower bound in the resource requirement for low-energy simulations.Here, we provide a lower bound on the query complexity to simulate any Hamiltonian H for low-energy simulations.In particular, we prove the following theorem: Theorem 6 (Lower bound).Given any positive number D, we can construct a 2-sparse positive semi-definite Hamiltonian H with ∥H∥ = Θ(D) and zero minimal eigenvalues such that simulating some state within the low-energy threshold ∆ and accuracy ϵ for some chosen scaled time τ requires at least Ω max τ, log(∆/ϵ) log log(∆/ϵ) queries to H.
Proof.For the construction of the Hamiltonian, we adopt the idea of proving the no-fastforwarding theorem for computing parity using quantum computation [37,38].We start with the following Hamiltonian We can observe that the Hamiltonian is actually an J x operator acting on a spin-D/2 system.The operator norm of the Hamiltonian is exactly D 2 .There are two eigenstates corresponding to the eigenvalue −D/2 and D/2, the former within which is the ground state.In addition, we can find that ⟨0| e iH 1 t/2 |D⟩ = (sin(t/2)) D .
This Hamiltonian has been utilized as a component to construct lower bounds for Hamiltonian simulation through the following construction H 2 [10,36].In particular, we consider the Hamiltonian H 2 generated by an D-bit string x = x 1 , . . ., x D acting on the quantum state |i, j⟩ with j ∈ {0, . . ., D} and i ∈ {0, 1}.The nonzero entries of H 2 read: ⟨i, j| H 2 i ′ , j + 1 = ⟨i, j + 1| H 2 i ′ , j + 1 = (D − i)(i + 1) 2 , j = 0, . . ., D − 1, where i ⊕ i ′ = x i+1 .This Hamiltonian, however, has a minimal eigenvalue of −D/2.In order to make the Hamiltonian positive semi-definite, we consider the following Hamiltonian which adds an additional e −i∆t global phase to the evolution of any input states.As we can always add back the inverse of a global phase to the input state when we fix t, we simply ignore the effect of the term e −i∆t in the following analysis.while the overlap between e iHt |0, 0⟩ and |D, p(x) ⊕ 1⟩ is strictly zero.Notice that the Hamiltonian is 2-sparse and we can only query at most two values of x i by querying H once. Therefore, simulating the evolution of input state |0, 0⟩ for time t = π yields an unbounded-error algorithm for computing the parity of x, thus requiring Ω(D) queries when we measure the output state in computational basis and consider the probability of obtaining |D, p(x)⟩ [10].However, there are two issues with this technique in our setting.The first issue is that the state |0, 0⟩ is not in the low-energy subspace.It has energy ⟨0, 0| H |0, 0⟩ = D/2 due to the second term DI/2.To address this issue, we consider the input state as a mixed state of the minimal eigenstate |ψ min ⟩ and |0, 0⟩ as This state has energy exactly Tr(Dρ in ) = ∆.The second issue is that, in approximated simulation, we allow an ϵ error tolerance.Here, we use the techniques in [36] and consider the choice of D such that the overlap between the evolution of the second term The lower bound for the number of queries to H is thus Ω(D) = Ω max τ, log(∆/ϵ) log log(∆/ϵ) , which is what we expect.

Theorem 2 .
Let H = L l=1 H l be an n-qubit k-local Hamiltonian with parameters M , J, and d (LM J = O(n)).By choosing the Trotter number

Corollary 1 .
Let H = L l=1 H l be an n-qubit k-local Hamiltonian with parameters M , J, and d (LM J = O(n)).By choosing the Trotter number

4 Figure 1 :
Figure 1: Quantum simulation in full Hilbert space vs. low-energy subspace: Trotter error induced by product formulas for a 2 × 6 homogeneous Heisenberg spin-1/2 model.The solid lines denote simulation errors in the full Hilbert space while the dashed lines denote errors in the low-energy subspace.The plot distinguishes different Trotter numbers r by the color of the lines.

Figure 2 :
Figure2: Full Hilbert space vs. low-energy subspace: simulation error induced by the qDRIFT and random permutation algorithms for localized Heisenberg spin-1/2 models.The solid lines denote errors in the full Hilbert space while the dashed lines denote errors in the low-energy subspace.We distinguish the total evolution time t or the order of the product formula p by the color of the lines.We also compare the numerical and the theoretical error scalings.(a) We plot the worst-case and the low-energy simulation errors of the qDRIFT algorithm as a function of step number r.(b) The same as the upper left subfigure, but for the random permutation approach at p = 2 and p = 4. (c) We plot the worst-case and the low-energy simulation errors of the qDRIFT algorithm as a function of system size n.(d) The same as the lower left subfigure, but for the random permutation approach at p = 2 and p = 4.

Figure 3 :
Figure3: Full Hilbert space vs. low-energy subspace: Error induced by the symmetry protection approach for Heisenberg spin-1/2 models.The solid lines denote errors in the full Hilbert space while the dashed lines denote errors in the low-energy subspace.We distinguish the experiments with different schemes by the color of the lines.We also compare the numerical and the theoretical error scalings.(a) We plot the worst-case and the low-energy simulation errors as a function of step number r.(b) We plot the worst-case and the low-energy simulation errors as a function of system size n.

Figure 4 :
Figure 4: Full Hilbert space vs. low-energy subspace: Trotter error induced by the product formulas for a 3 × 3 power-law model.The solid lines denote errors in the full Hilbert space while the dashed lines denote errors in the low-energy subspace, and the plot distinguishes interaction power α by the color of the lines.

Table 2 :
The notation table.