Locality and Digital Quantum Simulation of Power-Law Interactions

The propagation of information in nonrelativistic quantum systems obeys a speed limit known as a Lieb-Robinson bound. We derive a new Lieb-Robinson bound for systems with interactions that decay with distance r as a power law, 1/rα. The bound implies an effective light cone tighter than all previous bounds. Our approach is based on a technique for approximating the time evolution of a system, which was first introduced as part of a quantum simulation algorithm by Haah et al., FOCS’18. To bound the error of the approximation, we use a known Lieb-Robinson bound that is weaker than the bound we establish. This result brings the analysis full circle, suggesting a deep connection between Lieb-Robinson bounds and digital quantum simulation. In addition to the new Lieb-Robinson bound, our analysis also gives an error bound for the Haah et al. quantum simulation algorithm when used to simulate power-law decaying interactions. In particular, we show that the gate count of the algorithm scales with the system size better than existing algorithms when α > 3D (where D is the number of dimensions).

Lieb and Robinson's original proof applies only to shortrange interactions, i.e., those that act over a finite range or decay at least exponentially in space.However, interactions in many physical systems, such as trapped ions [28,29], Rydberg atoms [30], ultracold atoms and molecules [31,32], nitrogen-vacancy centers [33], and superconducting circuits [34], can decay with distance r as a power law (1/r α ) and, hence, lie outside the scope of the original Lieb-Robinson bound.Thus, understanding the fundamental limit on the speed of information propagation in these systems holds serious physical implications, including for the applications mentioned above.Despite many efforts in recent years [4][5][6][7], a tight Lieb-Robinson bound for such long-range interactions remains elusive.
In this paper, we derive a new Lieb-Robinson bound for systems with power-law decaying interactions in D dimensions.While our bound is not known to be tight, it has four main benefits compared to the best previous bound for such systems [6]: (i) It is tighter, resulting in the best effective light cone to date [Eq.(17)].(ii) The bound applies at all times, and not just asymptotically in the large-time limit.(iii) The framework behind the proof is conceptually simpler, with an easy-to-understand interpretation based on physical intuition.(iv) Our approach is potentially applicable to studying a wider variety of quantities, including connected correlators [35,36] and higher-order correlators (for instance, the out-of-time-ordered correlator [37,38] and the full measurement statistics of boson sampling [26,39]) as we discuss in Sec.VI.
In contrast to the previous long-range Lieb-Robinson bounds [4][5][6][7], which all relied on the so-called Hastings-Koma series [4], our approach is based on a generalization of the framework Haah et al. [27] (HHKL) introduced as a building block for their quantum simulation algorithm.The essence of their framework is a technique for decomposing the time evolution of a system into evolutions of subsystems, with an error bounded by the Lieb-Robinson bound for short-range interactions [1].We extend the HHKL framework to longrange interactions and to a more general choice of subsystems.Remarkably, these modifications enable us to derive a tighter Lieb-Robinson bound for long-range interactions than the one we use in the analysis of the decomposition [5].
Additionally, we return to the original motivation of Haah et al.'s framework: the digital simulation of lattice-based quantum systems.We generalize the HHKL algorithm to simulate systems with power-law decaying interactions.The algorithm scales better as a function of system size than previous algorithms when α > 3D, and the speed-up becomes more dramatic as α is increased.
The structure of the paper is as follows.In Sec.II, we state our main results and summarize the proof of the new Lieb-Robinson bound.In Sec.III, we lay out the precise mathematical framework for the proof and generalize the technique for decomposing time-evolution unitaries [27] to power-law decaying interactions and to more general choices of subsystems.After that, we present two applications of the unitary decomposition in Sec.IV and Sec.V, which can be read independently of each other.Specifically, in Sec.IV, we use the unitary decomposition to derive the improved Lieb-Robinson bound for long-range interactions.In Sec.V, we analyze the performance of the HHKL algorithm from Ref. [27] when ap-FIG. 1.A demonstration of the unitary decomposition in Lemma 1. Panel (a): the three disjoint regions A, B, C in D = 1 and D = 2 dimensions with A convex and compact.Panel (b): Lemma 1 allows the evolution of the whole system to be approximated by a series of three evolutions of subsystems.The horizontal axis lists the sites in each of the three sets A, B, C (not necessarily according to their geometrical arrangement, particularly in higher dimensions).Each box is an evolution for time t of a Hamiltonian supported on the sites the box covers.These evolutions can be forward (white fill) or backward (orange fill, with dagger) in time.
plied to simulating long-range interacting systems.We conclude in Sec.VI with an outlook for the future.

II. SUMMARY OF RESULTS
In this section, we summarize our main results for the case of a one-dimensional lattice.Without loss of generality, we assume that the distance between neighboring sites is one.The unitary decomposition technique in Sec.III is generalized from a similar result for short-range interactions in Ref. [27].We use it to approximate the evolution of a long-range interacting system ABC by three sequential evolutions of its subsystems AB, B, and BC (see Fig. 1).We assume that the interaction strength between any two sites in the system is bounded by 1/r α , with r being the distance between the sites and α a nonnegative constant.This restriction on the Hamiltonian norm also sets the time unit for the evolution of the system.
There are two sources of error in the approximation: one due to the truncation of the Hamiltonian of the system ABC (we ignore the interactions that connect A and C), and the other due to the Hamiltonians of the subsystems AB, B, and BC not commuting with each other.For a fixed value of α, if the distance between the two regions A and C (see Fig. 1a) is large enough, namely α, the two error sources have the same scaling with .To estimate the error, for example from the truncation, we sum over interactions connecting sites in A and C, and obtain a total error of O 1/ α−2 (in one dimension) for the approximation in the unitary decomposition (as shown in Appendix A 1).
In Sec.IV, we use the unitary decomposition to prove a Lieb-Robinson bound for long-range interactions that is stronger than previous bounds, including the one we use in the proof of the unitary decomposition.The subject of such a FIG. 2. A step-by-step construction of the unitary Ũ such that Ũ † OX Ũ ≈ U † T OX UT .Each box represents an evolution of the subsystem covered by the width of the box for a fixed time.The colors of the boxes follow the same convention as in Fig. 1.In panel (a), the unitary UT is written as a product of evolutions of the same system in M = 5 consecutive time slices.(b) The evolution in the last (bottom) time slice is decomposed using the method in Fig. 1, with the choice of subsystems A, B, C such that X is contained in A. The evolutions of the subsystems B and BC (hatched boxes) therefore commute with OX and cancel out with their counterparts from U † T , resulting in (c).In panel (d), we repeat the procedure for the secondfrom-bottom time slice, but note the different choice of A, B, C from panel (b).This difference is necessary to ensure that the evolutions of B and BC commute with the evolution(s) from the previously decomposed time slice(s).We then commute them through OX again and remove them from the construction of Ũ in panel (e).Repeatedly applying the unitary decomposition for the other time slices, we obtain the unitary Ũ in panel (f), which is supported on a smaller region than the original unitary UT .With a proper choice of the size of B, we can make sure that Y lies outside this region, and, therefore, Ũ commutes with OY .
bound is usually the norm of the commutator T O X U T evolved under a long-range Hamiltonian for time T and another operator O Y supported on a set Y that is at least a distance R away from the support X of O X .Here, we briefly explain the essence of the proof using a one-dimensional system with fixed α and large enough R, T α as an example.The strategy is to use the aforementioned unitary decomposition to construct another unitary Ũ such that (i) will be approximately zero, up to the error of our approximation.For fixed α, we consider M ∝ T equal time slices and use the unitary decomposition to extract the relevant parts from the evolution U T in each time slice.Each time we decompose a unitary, we choose the subsystems A, B, C so that only A overlaps with the supports of the unitaries from the previous time slices (see Fig. 2), and therefore the evolutions of B and BC can be commuted through O X to cancel their counterparts from U † T (Fig. 2b and Fig. 2d): (1) The remaining evolutions that contribute to the construction of Ũ are supported entirely on a ball of radius ∼ M around X, where is the size of B and is chosen to be the same in all time slices.By choosing ∼ R/M and M < R so that Y lies outside this ball, the commutator norm [O X (T ), O Y ] is at most the number of time slices multiplied by O 1/ α−2 , which is the decomposition error per time slice.Therefore, we obtain a Lieb-Robinson bound for long-range interactions in one dimension: where c lr,α is a constant that may depend on α, but not on T, R. Setting the commutator norm to a small constant yields the causal region inside the effective light cone: T R α−2 α−1 .For comparison, the previous best Lieb-Robinson bound produces a light cone T R α−2 α [6].Our bound is therefore tighter in the asymptotic limit of large R and large T , while its proof is substantially more intuitive than in Ref. [6].A more careful analysis (Sec.IV) shows that our light cone also becomes linear in the limit α → ∞, where the power-law decaying interactions are effectively short-range.Moreover, our bound works for arbitrary time T , while the bound in Ref. [6] applies only in the long-time limit.We provide a more rigorous treatment as well as a bound for D-dimensional systems in Sec.IV.Section V then then discusses the original motivation for the unitary decomposition-digital quantum simulation-in the case of long-range interactions that decay as a power law.For α > 2D, our analysis shows that the HHKL algorithm [27] requires only O T n(T n/ε) 2D α−D log T n ε twoqubit gates to simulate the evolution of a system of n sites arranged in a D-dimensional lattice for time T with an error at most ε.For large α, the gate count of the algorithm scales with n significantly better than other algorithms.

III. FRAMEWORK
In this section, we present the technique for approximating the time evolution of a system by evolutions of subsystems.We later use this technique to derive a stronger Lieb-Robinson bound (Sec.IV) and an improved quantum simulation algorithm (Sec.V) for systems with long-range interactions.
We consider n sites arranged in a D-dimensional lattice Λ ⊂ N D of size L = O n 1/D and D ≥ 1. Recall that, without loss of generality, we assume the spacing between neighboring lattice sites is one.This assumption sets the unit for distances between sites in the lattice.We shall embed the lattice Λ into the real space R D .The intersection X ∩ Λ therefore contains every lattice site in a subset X ⊂ R D .The system evolves under a (possibly) timedependent Hamiltonian H Λ (t) = ı,  h ı,  (t), with h ı,  (t) being the interaction between two sites ı,  ∈ Λ.Without ambiguity, we may suppress the time-dependence in the Hamiltonians.We say a system has power-law decaying interactions if h ı,  ≤ 1 ı−  α , where • denotes both the matrix and the vector 2-norms, for some nonnegative constant α and for all ı = .[Note that h ı, ı may have arbitrarily large norm.]For readability, we denote by H X = ı, ∈X h ı,  the terms of H Λ that are supported entirely on a subset X ∩ Λ, and by U X t1,t2 ≡ T exp −i t2 t1 H X dt the evolution unitary under H X from time t 1 to t 2 , where T is the time-ordering operator.We also denote by dist (X, Y ) the minimum distance between any two sites in X and Y , by X c = R D \ X the complement of X in real space, by ∂X the boundary of a compact subset X, by Φ(X) the area of ∂X, and by XY the union X ∪ Y .In the following, we keep track of how errors scale with time, distance, and α, while treating the dimension D as a constant.
We now describe how to approximate the evolution of the system to arbitrary precision by a series of evolutions of subsystems using a technique we generalize from Ref. [27].
Lemma 1.Let A, B, C ⊂ R D be three distinct regions with non-empty interiors such that A ∪ B ∪ C = R D .Let A be both compact (closed and bounded) and convex.We have for all α > D + 1.Here, v, c 0 ∈ R + are positive constants, γ is a constant that can be chosen arbitrarily in the range (0, 1), and = dist (A, C) is the distance between sets A and C.
We emphasize that this lemma applies to arbitrary sets A that are both convex and compact.The sets we focus on include D-balls and hyperrectangles in R D .The former geometry is relevant in the proof of our new Lieb-Robinson bound, the latter in the analysis of the HHKL algorithm for longrange interactions.
Lemma 1 allows us to approximate the evolution of a long-range interacting system ABC by that of subsystems AB, B, BC (Fig. 1).The features of the function ξ α ( ) are better understood by considering two limiting cases of physical interest.First, when α is finite and (the distance between A and C) is large compared to α, the function ξ α ( ) behaves like which decays only polynomially with .In the second limit, as α → ∞ for a large but finite , we recover from ξ α ( ) the exponentially decaying error bound e −γ -a trademark of finite-range interactions [1,27].
The proof of Lemma 1, while more general, bears close resemblance to the corresponding analysis for short-range interactions in Ref. [27].However, there are two key differences.First, in order to make the approximation in Lemma 1, some interactions between sites separated by a distance greater than are truncated from the Hamiltonian.While such terms vanish in a system with short-range interactions, here they contribute O Φ(A)/ α−D−1 to the error of the approximation.In addition, instead of the original Lieb-Robinson bound [1] which applies only to systems with short-range interactions, we use Gong et al.'s generalization of the bound for longrange interactions [5].The result is an approximation error that decays with polynomially as O Φ(A)/ α−D−1 , in addition to the exponentially decaying error that exists already for short-range interactions.Nevertheless, the error can always be made arbitrarily small by choosing to be large enough.
In Sec.III A below, we present the proof of Lemma 1.After that, we demonstrate the significance of Lemma 1 with two applications: a stronger Lieb-Robinson bound for long-range interacting systems (Sec.IV) and an improved error bound for simulating these systems (Sec.V).Both sections are selfcontained, and readers may elect to focus on either of them.

A. Error bound on the unitary decomposition
Here, we will outline the proof of Lemma 1. Similar to Ref. [27], we begin with an identity: Our aim is to approximate W t by U C 0,t † U B 0,t † U BC 0,t , from which Lemma 1 will follow.For that, we look at the generator of W t [27], i.e., a Hamiltonian G t such that for all time.Exact differentiation of W t yields [40,41] where H X:Y = i∈X,j∈Y h ij (t) denotes the sum of terms supported across disjoint sets X and Y , and δ trunc , δ overlap are error terms we now define and evaluate.Note that the first term in Eq. ( 9) is the generator of U C 0,t † U B 0,t † U BC 0,t -the unitary with which we aim to approximate W t .
In contrast to the approximation for short-range interacting systems in Ref. [27], there are two sources of error in Eq. ( 9).The first error term δ trunc arises after we discard H A:C from Eq. (7).For the short-range interactions in Ref. [27], this error vanishes when the distance between A and C is larger than the interaction range.However, in our case, there is a nontrivial truncation error associated with ignoring long-range interactions between A and C: (10) for α > D + 1, where c tr is a constant [Eq.(A5)], = dist (A, C) is the distance between A and C. The factor of 1/l α in the bound comes from the requirement that the twobody interactions decay as a power law 1/r α , while the term D is due to the sum over all sites in the D-dimensional set C. Another factor of Φ(A) arises after summing over the volume of A, which we assume to be a compact and convex set.The detailed evaluation of the norm is presented in Appendix A 1.
The other error, which we define to be δ overlap , is the result of the approximation used between Eqs. (8) and (9).In the former equation, the operator evolves under H AB + H C , whereas in the latter, it evolves under the reduced Hamiltonian H B + H C , thus incurring the error: To understand why δ overlap is small, recall that H B:C is the sum of terms h b, c that are supported on two sites b ∈ B and c ∈ C. Since the strengths of such terms decay as 1/r α (with r the distance between the sites b and c), the main contribution to H B:C -and thus to δ overlap -comes from the terms where b and c are spatially close to each other.But since the sets A, C are separated by a large distance , if the site b is close to C, then it must be far from A. Thus, the evolution of h b, c for a short time under H AB can be well-approximated by evolution under H B alone.In Appendix A 2, we make this intuition rigorous using Gong et al. [5]'s generalization of the Lieb-Robinson bound to systems with long-range interactions.
In the end, we obtain the following bound on δ overlap : where c ov is a constant [Eq.(A24)] and γ ∈ (0, 1) is a free parameter.The bound has contributions from two competing terms: one that decays polynomially with and another that decays exponentially.The polynomially decaying term is dominant for fixed α and large , whereas the exponentially decaying term prevails as α → ∞ for fixed .The errors δ trunc and δ overlap in approximating the generator G W combine to give an overall error in approximating From this, we obtain the error bound in Lemma 1, with c 0 = max{c tr , c ov }/v.
Before discussing applications of Lemma 1, we pause here to note that the Lieb-Robinson bound in Gong et al. [5] used in the above analysis is not the tightest-known bound for longrange interactions [6].Our use of this bound, however, does not lead to a suboptimal error bound in Lemma 1.For finite α, the error bound is dominated by the polynomially decaying term 1/ α−D−1 , which arises from the truncation error δ trunc rather than δ overlap .Therefore, this error term would not benefit from a tighter Lieb-Robinson bound.In the limit α → ∞, on the other hand, we shall see later that the lemma already reproduces the short-range Lieb-Robinson bound, which is optimal up to a constant factor.Thus, we expect that using stronger Lieb-Robinson bounds would produce no significant improvement for the error bound in Lemma 1.

IV. A STRONGER LIEB-ROBINSON BOUND
In this section, we will use Lemma 1 to derive a stronger Lieb-Robinson bound for long-range interactions.The first generalization of the Lieb-Robinson bound to power-law decaying interactions was given by Hastings and Koma [4].However, their bound diverges in the limit α → ∞, where the power-law decaying interactions are effectively short-range.Later, Gong et al. [5] derived a different bound that, in this limit, does indeed converge to the Lieb-Robinson bound for short-range interactions.While we used this bound in Sec.III to prove Lemma 1, we will also show that by using this lemma, we can in turn derive a Lieb-Robinson bound for longrange interactions that is stronger than the one in Gong et al.In fact, our bound produces a tighter effective light cone than even the strongest Lieb-Robinson bound for long-range interactions known previously [6].
Recall that the subject of a Lieb-Robinson bound is the commutator norm where O X , O Y are two operators supported respectively on two sets X, Y geometrically separated by a distance R, and U Λ 0,T is the time-evolution unitary of the full lattice Λ under a power-law decaying Hamiltonian, as defined above.
To compare different bounds, we analyze their effective light cones, which, up to constant prefactors, predict the minimum time it takes for the correlator C(T, R) to reach a certain value.For example, the original Lieb-Robinson bound [1] produces a linear light cone T R for short-range interactions.For long-range interactions, Hastings and Koma [4] first showed that C(T, R) ≤ ce vT /R α for some (α-dependent) constants c, v.By setting C(T, R) equal to a constant, the bound gives an effective light cone T log R in the limit of large T and R. Gong et al. [5] later achieved a tighter light cone that is linear for short distances and becomes logarithmic only for large R. Shortly after, Foss-Feig et al. [6] derived a bound with a polynomial light cone: Equation ( 14) was the tightest light cone known previously.
In the remainder of this section, we use Lemma 1 to derive a Lieb-Robinson bound for long-range interactions that produces an effective light cone tighter than the one in Ref. [6], while also using a much more intuitive approach.In addition, our bound works for all times, unlike the bound in Ref. [6], which applies only in the long-time limit.
Theorem 1 (Lieb-Robinson bound for long-range interactions).Suppose O X is supported on a fixed subset X.For α > 2D, we have Here R = dist (X, Y ) is the distance between the supports of O X and O Y , c lr , clr , v are constants that may depend only on D [defined in Appendix C], and ξ α is given by Eq. (3).
Before we prove Theorem 1, let us analyze the features of the bound.Although the general bound in Eq. ( 15) looks complicated, it can be greatly simplified in some limits of interest.For example, for finite α, in the limit of large vT > α and large R such that R/(vT ) α, the term algebraically decaying with R/(vT ) in ξ α (Rα/(vT )) dominates the exponentially decaying one [see also Eq. ( 3) and Eq.(C16)].Therefore, the Lieb-Robinson bound in this limit takes the form: where c lr,α is finite and may depend on α [Eq.(C18)].We can immediately deduce the effective light cone given by our bound for a finite α: which is tighter than Eq. ( 14) (as given by Ref. [6]).In particular, for α close to 2D, the exponent in Eq. ( 17) can be almost twice that of Ref. [6] (the larger the exponent, the tighter the light cone).
On the other hand, in the limit α → ∞, vT is finite and therefore always less than α.Hence our bound converges to the short-range bound C(T, R) ≤ 2c lr e vT −γR .We note that FIG. 3. A construction of the unitary Ũ which results in an improved Lieb-Robinson bound for long-range interactions in Theorem 1.The horizontal axes list the sites in each subset.Here Br denotes a D-ball of radius r centered on X, and Sr = B r+ \ Br a D-dimensional shell of inner radius r and outer radius r + , for some parameter to be chosen later.(See Fig. 7 in Appendix C for an illustration of the sets.)The evolution unitaries are represented by boxes with the same color convention as in Fig. 1.We first divide the interval [0, T ] into M = 5 equal time slices (upper panel).Note that because we consider OX (T ) in the Heisenberg picture, the vertical axis is therefore backward in time so that the bottom time slice will correspond to the first unitary applied on OX .The evolution during each time slice is approximated by three evolutions of subsystems using Lemma 1 (lower panel).The bottom two unitaries have their supports outside X and therefore commute with OX .They cancel with their Hermitian conjugates from Ũ † in Ũ † OX Ũ .Repeating the argument for higher time slices, we can eliminate some unitaries (hatched boxes) from the construction of Ũ .Finally, we are left with Ũ consisting only of unitaries (white boxes) that are supported entirely on the Dball B r 0 +5 of radius r0 + 5 .Therefore, Ũ commutes with OY , whose support lies in the complement B c r 0 +5 of B r 0 +5 . in this limit, the exponent of the light cone in Eq. ( 17) also converges to one, which corresponds to a linear light cone, at a linear convergence rate [see Eq. (C20) for details].These behaviors are naturally expected since a power-law decaying interaction with very large α is essentially a short-range interaction.
As mentioned earlier, we derive Theorem 1 by constructing a unitary Ũ such that (i) We note that Ũ does not necessarily approximate U Λ 0,T .It then follows from the two requirements that the commutator norm C(T, R), defined in Eq. ( 13), is upper bounded by the error of the approximation in (i).
We also note that the assumption on the norms of the interactions being bounded excludes several physical systems whose local dimensions are unbounded, e.g.bosons [see Ref. [42,43] for discussions of information propagation and Lieb-Robinson bounds in these systems].However, our Lieb-Robinson bound may still apply if the dynamics of the systems can be restricted to local Hilbert subspaces which are finite-dimensional.Examples of such situations include trapped ions in the perturbative regime [29] and noninteracting bosons [26].
To construct Ũ , we use Lemma 1 to decompose the unitary U Λ 0,T into unitaries supported on subsystems, each of which either contains X or is disjoint from X.The unitaries of the latter type can be commuted through O X to cancel out with their Hermitian conjugates from U Λ 0,T † .The remaining unitaries form Ũ , which is supported on a smaller subset than U Λ 0,T .In particular, with a suitable decomposition, the support of Ũ can be made to not contain Y , and, therefore, Ũ commutes with O Y .The step-by-step construction of the unitary Ũ has also been briefly described earlier in Sec.II and in Fig. 2, using the specific case of a one-dimensional system with a finite α.This construction immediately generalizes to higher dimensions and to arbitrary α, including the α → ∞ limit.The construction of Ũ for arbitrary D is summarized in Fig. 3.
We note that there is more than one way to decompose the unitary U Λ 0,T in the construction of Ũ .Different constructions of Ũ result in different approximation errors, each of which provides a valid bound on the commutator norm C(T, R).Therefore, the goal is to find a construction of Ũ with the least approximation error.In Appendix C, we present the construction that results in the bound in Theorem 1.Although we have evidence suggesting that the construction is optimal, we do not rule out the existence of a better construction.

V. BETTER PERFORMANCE OF DIGITAL SIMULATION
In this section, we generalize the algorithm in Ref. [27] to simulating long-range interactions.In general, the aim of quantum simulation algorithms is to approximate the time evolution unitary U Λ 0,T using the fewest number of primitive, e.g.two-qubit, quantum gates.Here, we show that in addition to the stronger Lieb-Robinson bound presented in the previous section, Lemma 1 can also be used to perform error analysis for the HHKL algorithm (Ref.[27]) in the case of interactions that decay as a power law, therefore improving the theoretical gate count of digital quantum simulation for such interactions.
Using the best known rigorous error bounds, simulations based on the first-order Suzuki-Trotter product formula [44] use O T 2 n 6 /ε gates to simulate the evolution U Λ 0,T of a time-dependent Hamiltonian on n sites up to a fixed error ε.(In this section, the big O is with respect to n, T, and 1/ε.)The generalized (2k)th-order product formula uses O n 2 (T n 2 ) 1+1/(2k) /ε 1/2k quantum gates.While this scaling asymptotically approaches O T n 4 as k → ∞, it suffers from an exponential prefactor of 5 2k [45].More advanced algorithms, e.g., those using quantum signal process-FIG.4. A demonstration of the HHKL decomposition [27] of the evolution of a fixed time interval for a system with m = 10 blocks, each consisting of sites.As before, each box represents a unitary (white) or its Hermitian conjugate (orange) supported on the covered sites.Using Lemma 1, the HHKL decomposition approximates the evolution of the whole system [panel (a)] by three unitaries supported on subsystems [panel (b)].By applying Lemma 1 repeatedly [panels (c) and (d)], the evolution of the whole system is decomposed into a series of evolutions of subsystems, each of size at most 2 .ing (QSP) [46] or linear combinations of unitaries (LCU) [47], can reduce the gate complexity to O T n 3 log(T n/ε) .Our error analysis below (Lemma 2) reveals that, when α is large, the number of quantum gates required by the HHKL algorithm to simulate long-range interactions scales better as a function of the system size than previous algorithms.
The HHKL algorithm itself uses either the QSP algorithm or the LCU algorithm as a subroutine to simulate the dynamics of a subset of the sites for one time step.Although the QSP algorithm does not handle time-dependent Hamiltonians, LCU can be applied to time-dependent Hamiltonians.Our results assume that (i) the local terms h ı, ı (t) have bounded norms for all ı ∈ Λ, and (ii) the Hamiltonian H Λ (t) varies slowly and smoothly with time so that h |X| ≡ max t ∂H X t /∂t exists and scales at most polynomially with |X| for all subsets X ⊂ Λ.These restrictions allow portions of the system to be faithfully simulated using LCU (or QSP, for a timeindependent Hamiltonian).
A. HHKL-type algorithm for simulating long-range interactions Although Ref. [27] focused on simulating short-range interactions, their (HHKL) algorithm can also be used to simulate long-range interactions.Here, we analyze the performance of their algorithm in simulating such systems.In the HHKL algorithm [27], the evolution of the whole system is decomposed, using Lemma 1, into elementary unitaries, each evolving a subsystem of at most (2 ) D sites, where is again a length scale to be chosen later.For a fixed time t, the algorithm simply simulates each of these elementary unitaries using one of the existing quantum simulation algorithms.In particular, we shall use LCU or (for a time-independent Hamiltonian) QSP due to their logarithmic dependence on the accuracy.
In this section, we assume α is finite and analyze the gate count in the limit of large system size n α.As a consequence, the block size can also be taken to be much larger than α.For simplicity, we will not keep track of constants that may depend on α.Recall that in this limit, the error bound in Lemma 1 is at most where we have assumed t = O (1).Using Lemma 1, we obtain the error bound for the first step of the HHKL algorithm, which can be summarized by the following lemma.

Lemma 2 (HHKL decomposition).
There exists a circuit that approximates U Λ 0,T up to error O T n/ α−D , where ≤ n 1/D /2 is a free parameter.The circuit has depth at most 3 D T and consists of O T n/ D elementary unitaries, each of which evolves a subsystem supported on at most (2 ) D sites for time t = O (1).
Proof.We now demonstrate the proof by constructing the circuit for a one-dimensional lattice (Fig. 4).A generalization of the proof to arbitrary dimension follows the same lines and is presented in Appendix D.
First, we consider M ∝ T equal time intervals then naturally decomposes into M consecutive simulations of U Λ tj ,tj+1 .We then divide the system into m consecutive disjoint blocks, each of size = n/m (Fig. 4).Denote by L k (k = 1, . . ., m) the set of sites in the k-th block.Using Lemma 1, we can approximate This approximation can be visualized using the top two panels of Fig. 4. Repeated application of Lemma 1 yields the desired circuit (bottom panel of Fig. 4), with each elementary unitary evolving at most 2 sites for time t.
To obtain the error estimate in Lemma 2, we count the number of times Lemma 1 is used in our approximation.In each of the M time slices, we use the lemma O (m) = O (n/ ) times, each of which contributes an error of O 1/ α−2 [see Eq. ( 18) with Φ = O (1) in one dimension].Therefore, with M ∝ T , the error of using the constructed circuit to simulate as given in Lemma 2.
The error bound for the approximation in Lemma 2 leads to an upper bound on the gate complexity of digital quantum simulation, as stated in the following theorem.
Theorem 2. For α > 2D, there exists a quantum algorithm for simulating U Λ 0,T up to error at most ε with gate complexity This gate complexity can be achieved by applying the HHKL algorithm [27] for long-range interactions, as described above.First, the evolution of the whole system U Λ 0,T is approximated by O T n/ D elementary unitaries as provided in Lemma 2. Each of these elementary unitaries is then simulated using one of the existing algorithms, e.g., LCU, with error that we require to be at most ε D /T n.If the Hamiltonian is time-independent, one can also use the QSP algorithm to simulate the elementary unitaries.
In the decomposition of the evolution, the accuracy of the approximation can be improved by increasing the block size . By Lemma 2, to achieve an overall error at most ε, we need When simulating the elementary unitaries, since each is an evolution of at most (2 ) D sites for time t = O (1), the LCU algorithm with error at most ε D /T n uses O 3D log T n ε h D two-qubit gates [45].Recall that we assume h D scales at most polynomially with D .With the block size from Eq. ( 22), we find the total gate complexity of simulating the O T n/ D elementary unitaries is The scaling of G D as a function of the system size n is significantly better than existing algorithms for large α.In particular, at T = n, this HHKL algorithm for long-range interactions requires only O n 2+ 4D α−D log n gates, while algorithms such as QSP or LCU use O n 4 log n gates or more.Therefore, the algorithm provides an improvement for FIG. 5.The gate count for simulating the dynamics of a onedimensional Heisenberg chain [Eq.( 25)] of length n, with α = 4, T = n, and ε = 10 −3 .We compare the gate count of the HHKL algorithm (orange square) to the QSP algorithm (blue circle).Note that the HHKL algorithm based on Lieb-Robinson bounds also uses the QSP algorithm as a subroutine.We obtain the scatter points using the approach described in Appendix E and fit them to a powerlaw model (solid lines).The asymptotic scalings of the gate count obtained from the power-law fits (n 3.29 for HHKL, n 4.00 for QSP) agree well with our theoretical predictions (see Table I).
α > 3D.However, the gate complexity of the algorithm depends polynomially on 1/ε, in contrast to the logarithmic dependence achieved by QSP and LCU, and by the HHKL algorithm for systems with short-range interactions.While this poly(1/ε) scaling is undesirable, in practice, the total error of the simulation is often set to a fixed constant (for example, see Ref. [48]) and effectively the dependence of ε only contributes a prefactor to the gate complexity of the algorithm.
As an example, in Fig. 5, we estimate the actual gate count of the HHKL algorithm in simulating a Heisenberg chain [Eq.(25)] and compare it with the gate count of the QSP algorithm (up to the same error tolerance).Because of the poly(1/ε) overhead, the HHKL algorithm based on Lieb-Robinson bounds uses more quantum gates for simulating small systems, but eventually outperforms the QSP algorithm when the system size n is large.
It is also worth noting that, in the limit α → ∞, the gate complexity becomes O (T n log (T n/ε)), which coincides (up to a polylogarithmic factor) with the result for short-range interactions in Ref. [27].This behavior is expected, given that a power-law decaying interaction with α → ∞ is essentially a nearest-neighbor interaction.However, we caution readers that at the beginning of this section, we have assumed that α is finite so that n α.Hence, the gate count in Eq. ( 24) is technically not valid in the limit α → ∞.Nevertheless, the error bound in Lemma 1 reproduces the estimate for short-range interactions in Ref. [27], and therefore, repeating the argument of this section in the limit α → ∞ should also reproduce the gate count for simulating short-range interactions in Ref. [27].

B. Numerical evidence of potential improvement
Up to now, we have seen that Lieb-Robinson bounds can improve the error bounds of quantum simulation algorithms, as demonstrated by the HHKL algorithm.We now provide numerical evidence hinting at the possibility of further improving the error bounds.
Although the HHKL algorithm outperforms previous ones when α > 3D, it remains an open question whether there is a faster algorithm for simulating long-range interactions.We also note that the gate complexities are only theoretical upper bounds, and these algorithms may actually perform better in practice [49].
As an example, we compute the empirical gate count of a Suzuki-Trotter product formula simulation of a onedimensional long-range interacting Heisenberg model where B j ∈ [−1, 1] are chosen uniformly at random and σ j = (σ x j , σ y j , σ z j ) denotes the vector of Pauli matrices on the qubit j.Specifically, we consider a simulation using the fourth-order product formula (PF4).We use a classical simulation to determine the algorithm's performance for systems of size n = 4 to n = 12 for time T = n, and extrapolate to larger systems.For each n, we search for the minimum number of gates for which the simulation error is at most ε = 10 −3 .We plot in Appendix F this empirical gate count, which appears to scale only as O n 3.64 with the system size n.We list in Table I the gate counts of several popular algorithms for comparison.The theoretical gate complexity of PF4 is O n 5.75 [44], while the QSP and LCU algorithms both have complexity O n 4 log n .These numerics show that the PF4 algorithm for simulating long-range interacting systems performs better in practice than theoretically estimated; in fact, it even performs almost as well as the HHKL algorithm based on Lieb-Robinson bounds [which scales as O n 3.33 log n by our earlier analysis].Whether other quantum simulation algorithms, including the HHKL algorithm, can perform better than suggested by the existing bounds remains an important open question.

VI. CONCLUSION & OUTLOOK
To conclude, we derived an improved bound on how quickly quantum information propagates in systems evolving under long-range interactions.The bound applies to powerlaw interactions with α > 2D, such as dipole-dipole interactions in 1D (often realizable with nitrogen-vacancy centers [33] or polar molecules [32]), trapped ions in 1D [28,29], and van-der-Waals-type interactions between Rydberg atoms [30] in either 1D or 2D.For finite α > 2D, our Lieb-Robinson bound gives a tighter light cone than previously known bounds-including the one used in the proof of Lemma 1.As of yet, we are not aware of any physical systems

Algorithm
Scaling with n = T Scaling with ε Empirical PF4 O n 3.64 -Our HHKL bound O n 3.33 log n O log(1/ε)/ε 0.67 PF4 bound [45] O n 5.75  O 1/ε 0.25 QSP bound [46] O n 4 log n O (log(1/ε)) LCU bound [47] O n 4 log n O (log(1/ε)) TABLE I.A comparison between the gate complexity of several quantum simulation algorithms for simulating one-dimensional power-law systems at T = n and α = 4.Our analysis shows that the HHKL algorithm performs at least as well as the empirical gate count of PF4, while having a similar poly(1/ε) scaling with the error ε.It is not known whether the empirical gate count of PF4 can scale with ε better than suggested by the best proven bound (the third row).
that saturate the Lieb-Robinson bounds for power-law interactions, including the new bound.In the limit α → ∞, our bound asymptotically approaches the exponentially decaying bound for short-range interactions.Our bound gives a linear light cone only in this limit, however, and it remains an open question whether there exists a stronger bound with a critical α c such that the light cone is exactly linear for α ≥ α c [50].
Currently, there are no known methods for quantum information transfer that are faster than linear for α ≥ D + 1.It is possible, therefore, that a stronger bound exists with a finite α c ≥ D + 1.It is our hope that the present work, as well as the techniques that we use, will help motivate the search for such stronger bounds.Our technique immediately extends the HHKL algorithm in Ref. [27] to the digital quantum simulation of the above systems.Our error bounds indicate that the gate complexity of the algorithm is better than that of other state-of-the-art simulation algorithms when α is sufficiently large (α > 3D), and matches that of the short-range algorithm when α → ∞.
However, the empirical scaling of other algorithms-such as product formulas-indicates that this gate complexity may only be a loose upper bound to the true quantum complexity of the problem.While a matching lower bound for the gate complexity of the HHKL algorithm is provided in Ref. [27] for Hamiltonians with short-range interactions, we do not know of any techniques that could provide a corresponding bound for long-range interactions.In addition to improving the quantum gate complexity, our results may also aid in the design of better classical algorithms for simulating long-range interacting quantum systems.In particular, while we still expect the classical gate complexity to be exponential in the simulation time, there may be room for a polynomial improvement.
While the use of Lieb-Robinson bounds to improve the performance of quantum algorithms is a natural extension of Haah et al., the opposite direction-using quantum algorithms to improve Lieb-Robinson bounds-is new.The connection from quantum simulation algorithms to Lieb-Robinson bounds that we have established opens another avenue for the condensed matter and atomic/molecular/optical physics communities to potentially benefit from future advances in quantum algorithms.In addition to proving a stronger Lieb-Robinson bound, the tools we developed may help to answer other questions regarding both short-range and long-range interacting systems.Using the same unitary construction as Theorem 1, we can generalize the bounds on connected correlators [35,36] to long-range interacting systems.Our results can also provide a framework for proving tighter bounds on higher-order commutators, such as out-of-time-order correlators [37,38].Previous methods used to derive Lieb-Robinson bounds-due to their use of the triangle inequality early in their proofs-have not been able to capture the nuances in the growth of such correlators.In addition to the more intuitive proof of the Lieb-Robinson bounds, our framework can be used to provide an alternative, simpler proof of the classical complexity of the boson-sampling problem [39], which generalizes the result in Ref. [26] to long-range interactions and also to more general Hamiltonians with arbitrary local interactions [51].By taking advantage of the unitary decomposition in Lemma 1, we obtain a longer time interval within which the sampler in Ref. [26] is efficient [51].Moreover, by generalizing from two-body to many-body interactions, our technique may find applications in systems whose Hamiltonians include interaction terms between three or more sites, e.g.many-body localized systems in the l-bit basis [52].
Note added: Shortly after we submitted our paper to arXiv, Else et al. [53] posted their work on a different Lieb-Robinson bound for power-law decaying interactions.For Hamiltonians consisting of at most two-body interactions, the bound in Ref. [53] and our bound both apply in the same regime, α > 2D [54].Within this regime, our bound results in a strictly tighter light cone than Ref. [53].However, the bound in Ref. [53] also applies to Hamiltonians consisting of k-body interactions, for any integer k.Generalizing our framework to cover such k-body interactions would be an interesting future direction.use the following lemma, which generalizes a similar lemma in Ref. [27] to arbitrary, time-dependent Hamiltonians.Lemma 3. Let Ω ⊂ Λ be a subset of sites.Let H Ω (t) = i,j∈Ω h ij (t) be the terms of H Λ (t) supported entirely on Ω.Let O X (τ ) be an observable supported on a subset X at a fixed time τ .We have: where Proof.To prove the lemma, we shall move into the interaction picture of H Ω (t) and treat ds be respectively the Hamiltonian and the evolution operator in the interaction picture.We have: Thus, Lemma 3 follows.
By substituting Λ → AB, Ω → B, O X → H B:C , τ → t, and noting that operators supported on disjoint subsets commute, we can show using Lemma 3 that where the observables are, in general, evaluated at different times s ≤ t.We note that while it is necessary to keep track of s, t for completeness, one should pay more attention to the supports of the operators, as they carry useful information about the locality of the system.In fact, let us pause for a moment to discuss why the right hand side of Eq. (A11) should be small when , the distance between A and C, is large.Whenever the supports of h a, b (s) FIG. 6.An illustration of h ab and h bc in a one-dimensional lattice.For short-range interactions, the sets {a, b } and {b, c} are separated by a distance of the same order as the size of B (upper figure).The contributions from these terms to δoverlap are bounded using a Lieb-Robinson bound.However, for long-range interactions, {a, b } and {b, c} can be geometrically close to each other (lower figure).In such cases, the norms of h ab and h bc decay as |b − a| −α and |c − b| −α and, therefore, their contributions to δoverlap are small.and h b, c (t) are far from each other, we can bound their commutator norm using a Lieb-Robinson bound for long-range interactions.We use the bound by Gong et al. [5]: where r = dist { a, b }, { b, c} is the distance between the supports, γ ∈ (0, 1) is a constant that can be made arbitrarily close to 1, and c, v are constants that depend only on D.
However, in contrast to short-range interacting systems, here b, b run over all possible sites in B, so in principle the distance between the supports of h a, b (s) and h b, c (t) can be small (Fig. 6).Fortunately, if that is indeed the case, then although the Lieb-Robinson bound is trivial, the assumption that h a, b and h b, c fall off as b − a −α and c − b −α , respectively, makes the summand in Eq. (A11) small.
Let us now evaluate the sum in Eq. (A11).In the following, we shall consider b = b , since the estimation for the case b = b follows a similar, but less complicated argument.Using Gong et al.'s Lieb-Robinson bound for long-range interactions [5]: where γ ∈ (0, 1) is a constant that can be chosen arbitrarily close to 1, while c, v are finite and bounded constants for all α, and is the distance between the supports of h b, c (t) and h a, b (s) (see Fig. 6).Since each term of δ overlap contributes a sum of an algebraically decaying as 1/r α and an exponential decaying as e −γr terms, it is convenient to evaluate their contributions separately.
First, let us find the contribution from the algebraically decaying part.It is straightforward to find out their contributions to δ overlap when r takes one of the four allowed values.Depending on which value r takes, we use either Lemma 5 or Lemma 7 in Appendix G to evaluate the sums over b and b .For example, the contribution from the terms where r = b − b is at most where λ 3 , λ 4 are constants that arise after we use Lemma 7 in Appendix G twice to evaluate the sums over b and b consecutively, and the sums over a, c have been bounded in the previous section (see Eq. (A2)).The constant λ 5 absorbs both λ 3 , λ 4 and the constants from the sums over a, c.
On the other hand, if r = b − c , we use Lemma 7 to evaluate the sum over b and Lemma 5 for the sum over b: where λ 6 , λ 7 come from the uses of Lemma 7 and Lemma 5 respectively.The constant λ 8 absorbs both λ 6 , λ 7 and the constants from the sums over a, c.Repeating for the other values of r, we find that the contribution from the algebraically decaying terms in Equation (A13) to δ overlap is at most for some constant λ 9 .
Next, let us find the contribution from the exponentially decaying term in Equation (A13).If r = b − b , we have where we have applied Lemma 8 in Appendix G to obtain the first inequality, Lemma 8 twice again and Lemma 7 to get the second inequality, and then Lemma 5 and Lemma 6 for the sums over a, c similarly to Appendix A 1. The constants λ 10 , λ 11 arise from the applications of the lemmas and are absorbed into a constant λ 12 .We note that the constant γ in the last three lines are different from the one in the first line (see Lemma 8 for details).However, they both are constants that can be chosen arbitrarily between 0 and 1.Therefore, we denote them by the same constant γ for convenience.
Repeating the argument for other choices of r in Eq. (A14), we find that the contribution from the exponentially decaying terms to δ overlap is still at most the right hand side of Eq. (A22).
Combining Eq. (A21) and Eq.(A22), we have for a constant λ 13 .Since D−1 ≤ (D−1)!ε D−1 e ε for any arbitrary small positive constant ε, we can upper bound where we have absorbed ε into the definition of γ and c ov .This completes the estimation of δ overlap .
placed by Ũ † 0 O X Ũ0 , we can approximate for some Ũ1 supported entirely on B r0+2 .The error of this approximation is at most c 0 e vt Φ(B r0+ )ξ α ( ).
By induction to all k = 2, . . ., M − 1, we can construct Ũ = Ũ0 Ũ1 . . .ŨT −1 such that where the overall error is at most and where we have replaced the surface area 2 ) r D−1 and M by T /t.Also by induction, the unitary Ũ is supported entirely on B r0+M .Therefore the lemma follows.
We are now ready to prove our Lieb-Robinson bound in Theorem 1.Without loss of generality, we assume the origin is in X.Since X = O (1), there exists r 0 = O (1) such that X is a subset of B r0 .By Lemma 4, there exists a unitary Ũ supported entirely on a D-ball B r with r = r 0 + M such that If we choose the number of time slices (M ) and the block size ( ) such that r ≤ R + r 0 , the set Y will lie outside the support B r of Ũ † O X Ũ , and therefore Ũ † O X Ũ will commute with O Y .Note that for a fixed value of M , the error should decrease with a larger value of .Therefore, to prove the strongest bound, we should choose as large as possible, i.e., = R/M , and hence M = R/ .Substituting the value of M and t = T /M into Eq.(C11), we obtain the error bound in terms of alone: where b 2 = b 1 r D−1 0 is a finite constant.We note that the above bound is valid for all values of ≤ R. The tightest bound can therefore be obtained by choosing a value for that minimizes the above expression.Our intuition and numerical evidence suggest that this happens when ∼ Rα vT , so in the below analysis, we aim to choose as close to this value as possible.
To proceed, we consider two regimes of time T , when vT ≥ α and when vT < α.In the former regime, we choose = Rα vT ≤ R and substitute into Eq.(C12) to get where we have used e α −1 α ≤ 1 for all α ≥ 1, 1 − α vT ≤ 1 and 1 + R ≤ 2R.In particular, if Rα vT > x 0 , where x 0 is the larger solution of x α−D−1 = e γx , the algebraically decaying term in ξ α dominates the exponentially decaying one, and therefore Combining Eqs.(C13) and (C16), we obtain a bound on the commutator norm: where The light cone implied by the bound is In the limit α → ∞, the exponent of the light cone converges to one at a rate given by On the other hand, if vT < α, we simply choose = R. Equation (C12) then becomes Therefore, we arrive at the Lieb-Robinson bound in Theorem 1 with clr = b 2 .
Similar to the D = 1 case, we first break the unitary into O (T ) unitaries exp(−iHt) for some t = O (1).We then use an algorithm consisting of D steps to break the simulation of exp(−iHt) into simulations of Hamiltonians on smaller hypercubes of size at most 2 .In the first of the D steps, we cut the D-dimensional lattice into L/ layers, each with the same thickness , a parameter to be chosen later.In this step, the cross section of the cut is L D−1 .Therefore, by Lemma 1, each time a new layer is generated, we accumulate an error of For T L/ layers of the first step, the accumulated error will be ε Next, for each of the O (T L/ ) layers of D −1 dimensions, we break them again into L/ layers of D − 2 dimensions.Using Lemma 1 with a cross section L D−2 , we find the error of the second step which decreases with faster than the error of the first step.More explicitly, in the kth of the D steps, the error is Appendix E: Estimation of the actual gate count In this section, we describe how we estimate the actual gate count of the HHKL algorithm and the QSP algorithm in simulating one-dimensional power-law systems.
The direct implementation of the QSP algorithm requires computing a sequence of rotation angles on a classical computer, which is prohibitive for large-size Hamiltonian simulation.Instead, we use a suboptimal approach described in Ref. [49].To simulate H = L j=1 β j H j for time t and accuracy ε, where L is the number of terms in the Hamiltonian, β j ≥ 0 and H j are both unitary and Hermitian, we divide the entire evolution into r segments.We choose r sufficiently large so that each segment is short enough for the classical preprocessing.Specifically, we choose r = j β j t τ max (E1) and τ max = 1000 [55].Within each segment, we choose q to be the smallest positive integer satisfying 4( j β j t/r) q so that the overall error is at most ε.This gives M = 2(q − 1) phased iterates within each segment [49].
The number of elementary operations of each phased iterate is log(L) + 4L + 8L.Here, the first term corresponds to the reflection along an L-dimensional state |0 ; the second term costs the preparation/unpreparation of an L-dimensional state; and the third term is the cost of selecting L two-body operators.We thus estimate the gate complexity of the QSP algorithm as log(L) + 12L rM .
Next, in order to determine the gate count of the HHKL algorithm, we need an estimate for the error of the unitary decomposition in Lemma 1. Recall that for D = 1, the error given by our analysis is b/ α−2 , where b is a constant that can be estimated numerically by computing the actual error for small values of and extrapolating for larger .
Since simulating the evolution of a generic system is classically intractable even for a moderate system size, we study only the one-dimensional Heisenberg model given in Eq. ( 25) and restrict our calculation to the single-excitation subspace.In Fig. 8, we plot the error of the unitary decomposition in Lemma 1 at several different values of (for system size n = 300 and evolution time t = 0.01).The scaling of the error agrees well with our prediction.By fitting the data to b/ α−2 , we obtain an estimate b = 1.62 × 10 −3 .
Recall that there are T /t time slices in the HHKL algorithm.In each time slice, there are n/ blocks of size and 2n/(2 ) blocks of size 2 .To meet the total error at most ε, we need to choose (see also Eq. ( 20)) By multiplying the number of blocks by the gate count for using QSP to simulate a single block, we arrive at the total gate count presented in Fig. 5.
Appendix F: Numerical performance of the product formula This section includes the numerical performance of the fourth-order product formula (PF4) used to simulate the evolution of the system given in Eq. ( 25) for time T = n.We plot this numerical performance as well as the theoretical estimates for the gate counts of the PF4, LCU, QSP, and HHKL algorithms in Fig. 9. FIG. 9.The empirical gate count of PF4 (purple dots) from n = 4 to n = 12, extrapolated to larger system sizes (solid, purple), for simulating the dynamics of the Hamiltonian in Eq. ( 25) for time T = n at a fixed error tolerance.The error bars are smaller than the size of the markers and hence not visible in the plot.Also shown in dashed lines are the slopes of the gate counts of several advanced algorithms for comparison.These slopes represent the scaling of the gate counts as functions of y-intercepts, which represent a constant multiplicative factor, should be ignored.

Appendix G: Mathematical tools
This section contains a collection of mathematical results omitted from the previous sections.In Appendix G 1, we present the upper bounds on standard sums we use in the proof of Lemma 1 in Appendix A. In Appendix G 2, we show how we estimate the sum over the convex set A in Eq. (A3) by parameterizing the elements of the set by their distance to the boundary of A. We also note that we use the same notation "λ" for constants that appear in different lemmas.

Standard sums
In this section, we present upper bounds on a a few standard sums used in the previous sections.Specifically, we use Lemma 5 to bound Eq. (A2), Eq. (A18), Lemma 6 to bound Eq. (A15), Eq. (A22), Lemma 7 to bound Eq. (A18), Eq. (A22), and Lemma 8 to bound Eq. (A22).Lemma 5. Let Λ be a D-dimensional lattice and r be the coordinates of sites in Λ.For α > D + 1 and R > √ D, there exists a constant λ that may depend on D but not on R, α such that: In particular, it implies that the sum r∈Λ converges for all α > D.
Proof.The proof of this bound is straightforward.For simplicity, we first assume none of the coordinates of r is zero.Since 1 x α is a decreasing function of x for all α > 0, we can always bound the sum over such r by an integral where denotes the sum over r with no zero coordinate and g(D) ≡ 2π ). Next, consider r with exactly one zero coordinate.These sites lie on D hyperplanes, each of dimension (D −1).Therefore the contribution from them can be evaluated using the above integral with D → D − 1: By repeating this argument for the sums over r with different number of zero coordinates, we arrive at Lemma 6.Let Λ be a D-dimensional lattice and r be the coordinates of sites in Λ.For all R > 0, there exists a constant λ that may depend on β, D but not on R such that: where β is a positive constant.In particular, it also implies that the sum r∈Λ converges.
Proof.The proof of this lemma follows the same idea as of Lemma 5.However, note that the function x β e −x is a decreasing function of x only when x ≥ x 0 for some x 0 that depends only on β.Therefore, if R ≥ x 0 , we follow the exact same lines as in the proof of Lemma 5.For example, if none of the coordinates of r is zero, we can bound for some constants λ 1 , λ 2 that depend only on β, D.
On the other hand, if R < x 0 , we consider The lemma should follow if we can argue that λ can be chosen independently of R. Indeed, since 1 ≤ R < x 0 and from the previous calculation, we know that the sum over r converges to a constant that depends only on β, D. This concludes the proof of Lemma 6.
where λ is a constant independent of a, c, α.
Proof.A proof of the lemma is presented in Ref. [4].
Proof.Without loss of generality, assume c = 0. Let = c − a = a be the distance between c and a.We need to Let B µ be a D-ball of radius µ centered around c for some arbitrary constant µ ∈ (0, 1).We shall divide the sum over b into two regimes, corresponding to b inside and outside B µ .
In the first regime where b is inside B µ , we can show using the triangle inequality that a − b ≥ (1 − µ) .Therefore, the sum over these b can be bounded by where we have used the fact that b∈B µ b β e b converges and is bounded by a constant λ which may depend only on D, β.
In the second regime, we bound a where the last sum is bounded using Lemma 6 and noting that µ < 1. Combining Eq. (G12), Eq. (G13), we arrive at a bound Note that if we choose µ = 1 2−γ , then γ = γ 2−γ takes on a value between 0 and 1, which can be arbitrarily close to 1.

Parameterizing a convex set
In this subsection, we show how we evaluate the sum over a in Eq. (A3).First, we parameterize a convex set by the distance to its boundary.The following lemma simplifies a sum over every site in a convex set to a sum over the above distance, multiplied by the boundary area of the set.Lemma 9. Let A ⊂ R D be a compact and convex set in R D with non-empty interior.Let C ⊂ R D be another subset dis-joint from A, and let = dist (A, C) be the smallest distance between elements of the two sets.Furthermore, we denote by a = dist ( a, C) the minimal distance from a given lattice site a in A to C. For a decreasing function f : R → R, we shall have where η is a constant that may depend only on D and Φ(A) is the boundary area of A.
Proof.Let us divide the set A ∈ R D into disjoint subsets S µ = { a ∈ A : µ ≤ dist ( a, ∂A) ≤ µ + 1} for µ = 0, 1, . . .Note that the assumption that the interior of A is non-empty implies that dist ( a, ∂A) is not uniformly zero.Roughly speaking, S µ contains the sites in A whose distances to the boundary ∂A are between µ and µ + 1.Therefore, a ≥ + µ for all a ∈ S µ .We then obtain where |S µ ∩ Λ| is the number of lattice sites that lie within S µ .
Let A µ = { a ∈ A : dist ( a, ∂A) ≥ µ} be a subset of A containing sites at least a distance µ from the boundary of A. Clearly, S µ = (A µ \ A µ+1 ) ∪ ∂A µ+1 and ∂S µ = ∂A µ ∪ ∂A µ+1 .Roughly speaking, S µ is a shell with the outer surface A µ , the inner surface A µ+1 and a unity thickness.The number of lattice sites in S µ will be bounded by ηΦ(S µ ) = η(Φ(A µ ) + Φ(A µ+1 )) (see Appendix G 2 a for the definition of the constant η).Since A is compact and convex, Φ(A µ+1 ) < Φ(A µ ) < Φ(A) (see Appendix G 2 b).Therefore, we arrive at the lemma.

a. The number of lattice sites in a compact region
In this subsection, we shall provide an upper bound on the number of lattice sites inside a compact set A ⊂ R D .We use this bound in Eq. (G18) to estimate the number of lattice sites in the set |S µ ∩ Λ| by its boundary area.Let A > = a ∈ A ∩ Λ : dist (a, ∂A) > 1 3 be the set of lattice sites that are at least a distance 1 3 away from the boundary ∂A, and let A ≤ = A \ A > be the other lattice sites of A.
First, note that for every lattice site a in A > , there exists a D-ball B 1/4 ( a) of radius 1 4 that contains no other lattice site and B 1/4 ( a) ⊂ A. Therefore, the number of lattice sites in Next, to count the lattice sites in A ≤ , we note that for every a ∈ A ≤ , we can select a point f ( a) ∈ ∂A on the boundary such that f ( a)− a ≤ 1 3 .We now argue that f ( a)−f ( b) ≥ Therefore, a D-ball B 1/6 (f ( a)) around f ( a) ∈ ∂A shall contain no f ( b) of any other lattice site b ∈ A ≤ .Therefore, the number of lattice sites in A ≤ is at most η 2 Φ(A), where Φ(A) = |∂A| is the boundary area of A and η 2 is the area of a (D − 1)-dimensional disk of radius 1/6.In summary, the number of lattice sites in A is therefore at most η 1 V (A) + η 2 Φ(A).In particular, for a shell A whose volume V(A) can be upper bounded by η 3 Φ(A), the number of lattice sites will be at most ηΦ(A), where η = η 1 η 3 + eη 2 .
b. Convex sets in R D are shrinkable In the proof of Lemma 9 [see the discussion after Eq. (G18)], we used the fact that Φ(A µ ) < Φ(A).In this section, we will show that this property of A-which we term shrinkability-holds if A belongs to the class of convex and compact sets in R D .The formal definition is as follows: Definition 1 (Shrinkable set).A compact set A ⊂ R D with boundary ∂A is shrinkable if, for all r > 0, A r = { a ∈ A : dist ( a, ∂A) ≥ r}, we have that Φ(A r ) = |∂A r | ≤ |∂A| = Φ(A).
In other words, a set is shrinkable if the surface area of the boundary of A r ⊆ A is no larger than that A. In this section, we will prove that convexity is a sufficient condition for shrinkability.Recall that a set is compact if it is both closed and bounded, whereas convexity is usually defined as follows: Definition 2. A set A is convex if for any x, y ∈ A and any θ such that 0 ≤ θ ≤ 1, we have θx + (1 − θ)y ∈ A.
Examples of convex sets include D-balls and hyperrectangles, which are also shrinkable.To prove this holds in general, we will first show that if A is convex, then A r is also convex (or empty) for all r > 0. To do this, we formulate an equivalent definition of a convex set as an intersection of halfspaces.From this definition, it follows that halfspaces are convex sets.A folk lemma [56] states that a closed set A is convex iff for some countable index set I. In other words, A is equivalent to the intersection of all halfspaces that contain it.Since convexity is preserved under arbitrary intersection, this implies that A is convex.The converse follows from the separating hyperplane theorem-see [56] for details.
With this equivalent definition of convexity in hand, we will prove that A r is also convex.To show that A is shrinkable, we must show that Φ(A r ) = |∂A r | ≤ |∂A| = Φ(A).Following a standard technique in the literature, we define the nearest-point projection of R D onto a convex set and then show that it is a contraction.The following lemma implies that such a mapping is well-defined.Proof.Since A is compact, the continuous function d x (y) = x − y must achieve its minimum value on A. Now suppose that minimum value of d x occurs at a point y ∈ A. We will show that y is unique.Assume for the sake of contradiction that there exists some point ỹ ∈ A such that d x (y) = d x (ỹ), but y = ỹ.Then the set of points x, y, and ỹ form an isosceles triangle, with y ỹ as the base.Dropping an altitude from x intersects this line segment at the midpoint m such that x − m < x − y = x − ỹ .But m is a convex combination of y and ỹ, i.e. m = 1 2 (y + ỹ) ∈ A, so we have reached a contradiction.Thus, y must be unique, and, therefore, p A (x) is well-defined.
The projection function p A (x) can be interpreted as generalizing the concept of the orthogonal projection into an affine subspace.It is also well-known that the nearest point projection p A is a contraction mapping.
Lemma 12.Given a nearest-point projection p A : R D → A onto a convex set A, it holds for all x, y ∈ R D that p A (x) − p A (y) ≤ x − y .
Proof.While the lemma can be proved for all x, y ∈ R D , for our purposes, we only need to consider x, y / ∈ A. Assume that p A (x) = p A (y). Then consider the hyperplanes H x and H y that pass through p A (x) and p A (y) respectively, and are perpendicular to the line segment p A (x)p A (y). (See the geometric diagram in Fig. 10.) We prove by contradiction that x (y) and p A (y) (p A (x)) lie on opposite sides of H x (H y ).Suppose without loss of generality that x and p A (y) lie on the same side of H x .Then the point where the altitude from x intersects the line segment p A (x)p A (y) would lie in A, contradicting the fact that p A (x) is the nearest-point in A to x.Thus, x (y) must lie on the opposite side of H x (H y ) from p A (y) (p A (x)).Then, as shown FIG. 10.The nearest-point projection of two points x and y onto a compact set A (oval).Also depicted are the line segment connecting the two image points pA(x) and pA(y), as well as the two hyperplanes orthogonal to it.in Fig. 10, the points x and y must fall outside the rectangular strip between the two hyperplanes.From this we conclude that p A (x) − p A (y) ≤ x − y .
The above result proves that the projection p A (x) is indeed a contraction.Since contraction mappings do not increase lengths, we can use this fact to demonstrate that the boundary of A r is less than that of A. Proof.Consider the projection p Ar : A → A r .Note that for r > 0, we have that A r = {x ∈ A | d(x, A c ) ≥ r} is entirely contained in the interior of A, which implies that A r ∩ ∂A = ∅.Thus, our situation satisfies the assumption we made in the proof of Lemma 12.
Under the action of p Ar , any point in R D outside of A r will get mapped to ∂A r .In particular, since the map is onto, ∂A will get mapped to ∂A r , i.e. p(∂A) = ∂A r .Using the fact that p Ar is contractive, we have that from which we conclude that A is a shrinkable set.This provides the final step in our proof of Lemma 9. Note that we do not require an explicit formula for the surface area of the boundary of a D-dimensional convex set.In general, one may use the Cauchy-Crofton formula to calculate this quantity-for more details, see Theorem 5.5.2 of Ref. [57].
, which is dominated by the error in the first step for all k > 1.Therefore, the error of cutting the D-dimensional lattice of size L into L D / D subsystems is still O T L D / α−D .To meet a fixed total error ε, we need to choose ∝ T L D /ε 1 α−D .The geometrical constraint < L requires α > 2D.Finally, simulating each of the O T L D / D subsystems using the LCU algorithm up to ε D /(T L D ) accuracy requires O 3D log T L D /ε D quantum gates.Therefore, the overall gate complexity of the algorithm is

1 3
for all distinct lattice sites a = b in A ≤ .Indeed, since a, b are distinct lattice sites, the least distance between them is 1, i.e. a − b ≥ 1.Using a triangle inequality, we can show that

Definition 3 .
A halfspace H is given by the points {x ∈ R D | a T x ≥ b}, where a ∈ R D \{0}.

Lemma 10 .
If a compact set A ⊂ R D is convex, then A r = { a ∈ A : dist ( a, ∂A) ≥ r} is convex (or empty) for all r > 0. Proof.Write A as the intersection of half-spaces H k = {x ∈ R D | a T k x ≥ b k }, for k ∈ I. Then A r is the intersection of the half-spaces given by H r k = {x ∈ R D | a T k x ≥ b k +r}.By the converse of the above lemma, A r is convex (or empty).

Lemma 11 .
Given a non-empty, compact and convex set A ⊆ R D and a point x ∈ R D , there exists a unique point p A (x) ∈ A such that p A (x) = arg min y∈A x − y .