Bounding the finite-size error of quantum many-body dynamics simulations

Finite-size error (FSE), the discrepancy between an observable in a finite system and in the thermodynamic limit, is ubiquitous in numerical simulations of quantum many body systems. Although a rough estimate of these errors can be obtained from a sequence of finite-size results, a strict, quantitative bound on the magnitude of FSE is still missing. Here we derive rigorous upper bounds on the FSE of local observables in real time quantum dynamics simulations initialized from a product state. In $d$-dimensional locally interacting systems with a finite local Hilbert space, our bound implies $ |\langle \hat{S}(t)\rangle_L-\langle \hat{S}(t)\rangle_\infty|\leq C(2v t/L)^{cL-\mu}$, with $v$, $C$, $c$, $\mu $ constants independent of $L$ and $t$, which we compute explicitly. For periodic boundary conditions (PBC), the constant $c$ is twice as large as that for open boundary conditions (OBC), suggesting that PBC have smaller FSE than OBC at early times. The bound can be generalized to a large class of correlated initial states as well. As a byproduct, we prove that the FSE of local observables in ground state simulations decays exponentially with $L$, under a suitable spectral gap condition. Our bounds are practically useful in determining the validity of finite-size results, as we demonstrate in simulations of the one-dimensional (1D) quantum Ising and Fermi-Hubbard models.

Introduction Numerical simulations are crucial to our understanding of many-body quantum matter, and are routinely applied in all fields of physics and in chemistry. Unfortunately, many numerical techniques popular in these fields incur significant FSEs when approximating properties of a large (potentially infinite) system by properties of a finite one. The most direct example is exact diagonalization (ED), which exactly solves the finite system numerically [1][2][3][4]. Accessible system sizes are limited since the Hilbert space dimension grows exponentially with system size; for the simplest case of interacting spin-1/2s, a state-of-the-art ground state calculation is limited to ∼ 45 spins [5]. FSEs also significantly affect other techniques, such as density matrix renormalization group (DMRG) [6][7][8][9], many tensor network algorithms [10,11], quantum dynamical typicality-based algorithms [12][13][14][15][16], and quantum Monte Carlo (QMC) [17], and they are a significant source of error for simulating quantum systems on quantum computers [18] and for analog quantum simulations using ultracold matter [19], trapped ions [20], and other platforms [21].
It is often difficult to characterize FSEs. The standard method to assess them is to calculate and compare observables for different system sizes, ideally using finitesize scaling [22]. Although useful, this method has limitations. One is that it offers no guarantees. Two different system sizes may have results that closely agree, but at larger sizes the physics changes and observables deviate [23]. Another is that one may not be able to study multiple system sizes that are sufficiently large to get a good estimate of the convergence.
In this paper, we derive rigorous upper bounds on the error of approximating observables in a large, possibly infinite, quantum many-body system by results in a smaller one. The bounds are applicable to arbitrary Hamiltonians for which a Lieb-Robinson (LR) bound exists. For quantum dynamics simulations starting from product initial states and evolving under locally interacting Hamiltonians with a finite local Hilbert space, the bound for a local observableŜ is where v, C, c, and µ are constants that can be computed explicitly and depend on the Hamiltonian, observable, and boundary condition. Such dynamics is explored in a wide variety of ultracold matter experiments, such as quantum quenches and slow ramps in Rydberg atoms [24][25][26][27][28][29], molecules [30][31][32], Fermi gases [33], atoms in optical lattices [34][35][36][37][38], and optical clocks [39]. This dynamics can probe fundamental phenomena, such as many-body localization [40][41][42][43][44], prethermalization [45,46], and generation of topological defects near critical points [47]. This bound is then extended to a large family of correlated initial states satisfying an exponential clustering condition. While our main focus is on dynamics, we also show that the FSE of local observables in a many-body ground state decays exponentially in system size, under a suitable spectral gap condition. The idea behind our bound is that locality -specifically that one piece of a system does not instantly affect far-away pieces -imposes strong constraints on quantum dynamics [48,49]. This can be seen by considering evolution under a Hamiltonian initiated from a product state (other scenarios can be understood by similar arguments). As illustrated in Fig. 1, an observable in a region X will be affected by FSEs only after a long enough time for information to propagate from the boundary to X. This idea is made precise by relating FSE to unequal FIG. 1. An illustration of our bounds. In locally interacting systems, information propagates no faster than the LR speed v, so it takes a finite amount of time tc ∼ dXB/v for the effect of the boundary, ∆Ĥ to affect the center site ob-servableŜ. Here dXB is the distance between the supports X (yellow square) and B (red crossed circles) ofŜ(t) and ∆Ĥ, respectively. time correlation functions, which can then be bounded by a LR bound [50], a direct consequence of locality. Although similar ideas of applying LR bounds to analyze the performance of some numerical algorithms have been employed in Refs. [51][52][53][54][55][56][57][58][59], the connection to FSE has not been made explicit, and the practical utility of the bounds for numerics was not demonstrated. This idea has also been applied to estimate FSE in a non-rigorous way, for example in Ref. [60].
Our FSE bound not only shows the convergence of finite-size approximations in principle, but is tight enough to be useful in practice, which we demonstrate in simulations of some prototypical models. For example, in the dynamics following a sudden change of parameters in a 1D transverse field Ising model (TFIM) with L = 21 sites, the error bounds for the transverse magnetization and nearest-neighbor correlations remain extremely small to times where they have evolved close to equilibrium. Furthermore, the bounds are reasonably tight: the time at which the error bound becomes significant is only 20-25% smaller than the time at which the actual FSE becomes noticeable. We similarly demonstrate this for the non-equilibrium relaxation of the Fermi-Hubbard model (FHM) from a checkerboard state, inspired by experiments and theory of Refs. [61,62]. The precision of these bounds is enabled by the major quantitative improvements offered by recent LR bounds [63,64].
In addition to their quantitative utility, these bounds provide insights into the convergence of numerical methods, and open the way to designing new algorithms. One immediate consequence of the bounds is to rigorously show that the FSE decays exponentially with the linear dimension of the system for periodic boundary condition (PBC), as well as for open boundary condition (OBC) provided that one measures observables only near the center of the system, as commonly employed in the DMRG community. If one instead averages the measurement over all sites in OBC, then our bound indicates that the error decays only algebraically. Similar behavior at finite temperature has been observed and analyzed in Ref. [65]. Furthermore, if one compares PBC to OBC with center site measurement, our error bound for PBC decays twice as fast with distance as the bound for OBC at early times, suggesting that PBC gives more reliable results at early times [66]. These insights may lead to new methods; one example is that they show why the movingaverage cluster expansion (MACE) method of Ref. [31] converges exponentially faster than alternative schemes.
A simple bound for both OBC and PBC Consider the dynamical evolution of a quantum many-body system on an infinite d-dimensional lattice, governed by a locally interacting HamiltonianĤ. For illustrative purpose, in Fig. 1 we draw the configuration for a 1D nearestneighbor interacting lattice model. Let |ψ be the initial product state,Ŝ be a local observable to be measured that acts on a finite region X (center point in Fig. 1), and let ∆Ĥ = jV j be the sum of all the interaction terms between the inner and outer parts of the system (red links in Fig. 1). If PBC are used, we further subtract from ∆Ĥ the interaction between the first and the last site (brown link in Fig. 1). LetĤ L and |ψ L denote the Hamiltonian and the initial state of the finite-size simulation, respectively (i.e. the restriction ofĤ and |ψ to the L-site inner system). DenoteĤ =Ĥ − ∆Ĥ, so that H decouples into two commuting terms, one acting only on the inner system, the other acting only on the outer system. The FSE of the observableŜ is where A ψ ≡ ψ|A|ψ , and we set = 1 throughout. SinceĤ decouples into two independent spatial regions (inner and outer) and |ψ is a product state, the first term in Eq. (2) can be rewritten as e iĤ tŜ e −iĤ t ψ . Inserted into Eq. (2), the two expectation values are taken in the same state |ψ , so their difference can be bounded by the operator norm | ψ|Â|ψ | ≤ Â . Using the unitary invariance of operator norm Â = ÛÂV for arbitrary unitary operatorsÛ ,V , we have whereÛ I (t) = e −iĤt e iĤ t is the evolution operator in the interaction picture, which satisfiesÛ I (0) = 1 and i∂ tÛI (t) =Û I (t)∆Ĥ(t), where ∆Ĥ(t) = e −iĤ t ∆Ĥe iĤ t . Now applying the fundamental theorem of calculus and the triangle inequality, we obtain a bound on the FSE: The integrand is the quantity bounded by LR bounds, so to upper bound the FSE, one can insert the relevant LR bound. We focus on locally interacting systems, but Eq. (4) applies equally to long-range interactions by substituting the corresponding LR bounds [58,[67][68][69][70][71][72][73] in those systems. For a locally-interacting system, the currently tightest LR bound is obtained by computing the series in Eq. (S19) of the Supplemental Material (SM) [74], which is based on Refs. [63,64], although this may not be efficiently computable in general. A slightly looser but efficiently computable method is discussed in Ref. [64], in which one solves a system of first order linear differential equations for a number of variables proportional to the system size. To see the qualitative features of the bound for large systems, we can insert the simple expression given in Eq. (3) of Ref. [64] into Eq. (4) to obtain where c j are constants independent of t and d Xj , D(Ŝ,V j ) is the distance between the operatorsŜ,V j in the commutativity graph (CG) as introduced in Ref. [64], and v is the LR speed. The distance in the CG is related to the distance in real space d Xj by D(Ŝ,V j ) = ηd Xj −µ, where η, µ are (straightforwardly determined) constants, and d Xj is the distance between X and j in real space. Therefore the rhs of Eq. (5) is bounded by ( vt d XB ) ηd XB −µ , where d XB = min j∈B d Xj . Despite its simplicity, the t-dependence of this bound generically agrees with the exact error to lowest order in t in OBC [74].
Besides its practical utility for bounding FSE in calculations, as demonstrated below, this result has qualitative implications. One is to rigorously support the common practice of measuring observables close to the center site in OBC numerics (e.g. in the DMRG community), rather than averaging over all sites. This minimizes the error bound, since the center size maximizes d XB . This choice yields our main result in Eq. (1) for the OBC case. Our bound allows one to extend this. For example, in dimension greater than one, we can minimize FSE by choosing an optimal cluster shape that minimizes the rhs of Eq. (5), and run simulations on the optimal shape.
An improved bound for PBC In the previous section we treated PBC in a way similar to OBC. But it turns out that the resulting bound in Eqs. (4) and (5) is qualitatively loose at small t for PBC. The reason for this can be intuitively understood as follows. The two terms in the rhs of Eq. (2) can be expanded in t. As we discuss in greater detail in the SM [74], the FSE forŜ(t) actually is only contributed by terms in e iĤtŜ e −iĤt whose spatial span is larger than L and terms in e iĤ L tŜ e −iĤ L t that wrap around the whole periodic system. The leading order of these terms is proportional to t L , where L is the length of the shortest non-contractible loop on the PBC commutativity graph, which is roughly twice as large as the exponent D(Ŝ,V j ) in Eq. (5). The SM [74] extends methods developed in Refs. [63,64] to derive a rigorous upper bound for |δ Ŝ (t) | that leads to this improved t L scaling. The main result is where the constant C p is given in Eq. (S52), v p is the LR speed in the p-th direction given in Eq. (S53), and L p is the size of the periodic system in the p-th direction in commutativity graph. L p is related to the real space system size L p by L p = η p L p − µ p for constant integers η p , µ p . We note that while this bound improves the smalltime exponent of the PBC bound by a factor of 2 compared to Eqs. (4) and (5), the timescale t c ≈ min p L p /2v p on which the bound exponentially grows is still approximately the same as Eq. (5). Besides its quantitative utility, Eq. (6) shows that in anisotropic systems where v p is different in each direction, one should choose L p ∝ v p in order to minimize the FSE. FSE in non-degenerate gapped ground states So far, we have been discussing FSEs of quantum dynamics simulations. We now derive a bound on FSE of local observables in non-degenerate ground states under a gap assumption. This result is interesting in its own right, and will also be useful for our subsequent generalization of the dynamics error bound to correlated initial states.
The HamiltonianĤ, observableŜ, boundary terms ∆Ĥ, etc. are the same as before. For convenience we suppose that the operatorŜ =Ŝ l has unit norm and acts nontrivially only within a region X = X l of diameter l which sits on the center of the finite-size cluster.
The difference now is that we consider the observable Ŝ = Tr[ρŜ] given by the ground state density matrix ρ. The numerical simulation approximates this thermodynamic quantity by the expectation value in the finite size ground state Tr[ρ LŜ ], whereρ L is the ground state density matrix ofĤ L . Our result relies on an assumption that the interpolated HamiltonainĤ(λ) ≡Ĥ − λ∆Ĥ is non-degenerate for all 0 ≤ λ ≤ 1 and has a uniform spectral gap min 0≤λ≤1 ∆(λ) = ∆ > 0. When this condition is satisfied, then analogously to Eq. (1) we have [74] Bounds for correlated initial states We now generalize our error bound to dynamics initiated from a class of (possibly mixed) initial statesρ, for which there exists a sufficiently good finite size approximationρ L satisfying Eq. (7). This includes non-degenerate gapped ground states (under the condition described above), but also includes translation invariant matrix product states (MPS) with a finite bond dimension [75], and finite temperature thermal states [76]ρ = e −βĤ /Tr[e −βĤ ] above a certain temperature, whereρ L = e −βĤ L /Tr[e −βĤ L ], i.e. the thermal state ofĤ L .
Given that we have an initial stateρ L satisfying Eq. (7), we can bound the dynamics FSE as is the reduced density matrix ofρ on the finite cluster, and in the second line we used the triangle inequality. The second term can be bounded using the same method as in Eq.
(2), since Tr[ρÂ] ≤ Â for any density matrix ρ. To bound the first term, we insert the expansion (8), and notice thatŜ l (t)−Ŝ l−1 (t) is an operator acting on X l , whose norm is bounded by Eqs. (4,5) to be Ŝ l (t) −Ŝ l−1 (t) ≤ C(2vt/l) ηl/2 for some constant C. For initial states satisfying Eq. (7), this implies where C 1 and C 2 are model-dependent constants that can be explicitly determined.
Example: 1D TFIM We test our dynamics error bounds in simulations of prototypical models for quantum many-body physics, starting with the TFIM, This is a canonical model for quantum phase transitions [77,78], and occurs in materials like CoNb 2 O 6 [79], cold atom [80][81][82] and trapped ion [83][84][85][86][87] experiments, and superconducting circuits [88,89]. We numerically study the dynamics of this model at the critical point J = h for several L, and calculate the exact evolution for L = ∞. Specifically, we study the dynamics of σ x (t) starting from |ψ(0) = | →→ . . . → . Analogous dynamics in the 2D TFIM has been explored in Rydberg atom experiments [27,28]. Fig. 2a shows σ x j (t) from L = 5 to 21 using PBC, along with the exact L = ∞ solution [90]. To obtain a FSE bound for σ x j (t) , we use the LR bound given in Eq. (S24) of Ref. [74] [obtained from the general bound Eq. (S20)], which, after inserting into Eq. (4), yields As Fig. 2a shows, this error bound provides a guarantee of the numerical calculations' accuracy out to interesting and useful timescales. For the L = 21-site calculation, the bound guarantees that the results are accurate (within 10 −2 ) up to times Jt ∼ 3.5, where the observable has nearly reached equilibrium. Furthermore, this time is in reasonable accord with the true time at which FSE becomes important (within 20%). We emphasize that the FSE bound never made use of the TFIM's exact solution. The bound Eq. (4) can be applied to any system including in dimensions greater than one. As we will now demonstrate in the 1D FHM, the bound still provides a useful guarantee of the finitesize results when no exact solution is available.
We numerically study the relaxation dynamics of a charge density wave state |ψ(0) = |202020 . . . 20 , analogous to previous theory [97] and experiments [98], where 2 (0) means a doubly occupied (empty) site. We run the finite-size simulations in OBC, and measure the density imbalanceM (t) = [N even (t) −N odd (t)]/L [97]. To get a FSE bound forM (t), we use the currently tightest LR bound, obtained by numerically summing the series in Eq. (S19) of the SM [74] and inserting the result into Eq. (4). Fig. 2b shows the results.
Our error bound can be compared to estimates of FSE obtained from comparing calculations of different sizes. For example, one can take the difference between the L = 10 and L = 12 as a rough estimate of the FSE of the L = 12 calculation. Our bound is comparable in its guaranteed timescale of convergence to this conventional estimate. For example, we can guarantee that the FSE in M (t) of the L = 12 result is less than 1% for Jt ≤ 1.2, comparable to the time Jt ∼ 1.6 where the L = 10 and 12 results differ noticeably.
Conclusions We have presented a rigorous upper bound on the FSE of local observables measured in numerical simulations of quantum dynamics starting from a large class of initial states. For product initial states, the bounds show an advantage of using PBC at early times. We also presented a generalization to simulations of local observables in non-degenerate gapped ground states. In all the cases we considered, the bounds decay exponentially in system size, and guarantee the accuracy of finite size dynamics simulation up to a time scale t c ∼ L/2v. These insights into FSE can motivate better algorithms.
The quantitative utility of the bounds is demonstrated in the 1D TFIM and 1D FHM. In both cases, the error bounds are extremely small up to timescales where there is interesting physics and even equilibration, and they are reasonably tight compared to the actual FSEs. We expect these bounds to provide useful tools to researchers going forward, providing FSE bounds on numerical calculations and suggesting new numerical methods that minimize this error.
We thank Miles Stoudenmire, Miroslav Hopjan, Bhuvanesh Sundar and Ian White for discussions, and Brian Neyenhuis for a careful reading of the manuscript. This work was supported in part by the Welch Foundation (C-1872), the National Science Foundation (PHY-1848304), and the Office of Naval Research (N00014-20-1-2695).
Supplemental Material for "Bounding the FSE of quantum many-body dynamics simulations" Zhiyuan Wang, Michael Foss-Feig, and Kaden R. A. Hazzard The supplemental material fills in technical details of the basic results in the main text. We first prove the ground state FSE bound in Eq. (7) under the gap assumption we mentioned there. Then we give the full derivation of the dynamics FSE bound for correlated initial state, Eq. (9). We then compare our dynamics error bounds to perturbation theory at early times and show that they are qualitatively tight for a large class of product initial states, i.e. they grow with time as t a with the correct exponent a. The remaining task is to give a detailed derivation of the improved error bound in PBC given in Eq. (6). An important intermediate step is to prove Eq. (S18), which is an analog of Eq. (4), expressing the FSE in terms of a LR commutator. We give three different methods to numerically upper bound the LR commutator that appears in the rhs of Eq. (S18). The first one is to evaluate the series in Eq. (S20), which combines the methods in Refs. [63,64]. This leads to the tightest bound, but is only efficiently computable in special cases. The second one is to numerically solve the differential equation (S30) and then calculate Eq. (S39). This method (which is based on Ref. [64]) is only slightly looser than the first one but is computationally efficient in general. The third method makes further simplifications which result in the simple analytic expression in Eq. (6), whose constants are given in Eqs. (S52,S53,S55). We emphasize that for any desired system, one may compute the LR bound and use it in the error formulas Eqs. (4) and (S18). In this way, as LR bounds are refined in the future or extended to more general systems (e.g. long-range interacting [71][72][73], bosonic [99], or continuum ones), these refinements can immediately be used in the FSE bounds.

FSEs in gapped non-degenerate ground states
In this section we prove the exponential decay of FSE in gapped non-degenerate ground states under the assumption stated in the main text. Recall thatĤ is the Hamiltonian in the thermodynamic limit, andĤ =Ĥ − ∆Ĥ is obtained fromĤ by removing the boundary links. The interpolated HamiltonianĤ(λ) = (1 − λ)Ĥ + λĤ =Ĥ − λ∆Ĥ is assumed to have a non-degenerate ground state |G(λ) and is uniformly gapped ∆(λ) ≥ ∆ > 0, for all λ ∈ [0, 1]. The basic idea is that, as the parameter λ is varied from 0 to 1, we keep track of how fast |G(λ) changes, as well as Ŝ l λ ≡ G(λ)|Ŝ l |G(λ) . Suppose |G(λ) is properly normalized such that G(λ)|G(λ) = 1 and phase chosen such that G(λ)| d dλ |G(λ) = 0 for any λ ∈ [0, 1]. Then first order non-degenerate perturbation theory gives d dλ where E 0 (λ) is the ground state energy ofĤ(λ) andP G (λ) ≡1 − |G(λ) G(λ)| is the projection operator to the space of excited states. Therefore, we have the integral formula for the FSE δS l : The integrand looks similar to the ground state correlator betweenŜ l and ∆Ĥ, so it's natural to guess that it should decay exponentially in a gapped system. This is verified by the following theorem (which is similar to the exponential clustering theorem in non-degenerate gapped ground states [67]): Theorem 1. LetÂ X ,B Y be arbitrary local observables with unit norm, supported on non-overlapping regions X, Y , respectively. Let |G be the unique ground state of locally-interacting HamiltonianĤ with spectral gap ∆. Then the quantity f XY ≡ G|Â XP G E0−ĤB Y |G + c.c. is upper bounded by where v is the LR velocity, W (x) is the Lambert-W function, C is a constant depending onÂ X andB Y , and d XY is the distance between X and Y .
Proof. Using the identity where α > 0 is a parameter to be specified later, we have where we use the notation (x) ≤1 = min{x, 1}, (x) ≥0 = max{x, 0}, C is a constant coming from the LR bound that does not depend on r, α, C 2 , C 3 are coefficients that only weakly depend on r, α, and in the last line we substituted T = τ r and α = 1/(λr). Notice that Eq. (S6) holds for arbitrary positive λ, since the parameter α > 0 introduced in Eq. (S4) can be chosen arbitrarily. If we choose λ to be any positive value, we can immediately prove the exponential decay of |f XY | in d XY which leads to Eq. (7) in the main text, since all the three terms in the rhs of Eq. (S6) decay exponentially in r. If we want a better bound, we can choose λ to maximize the smallest decay coefficient min{λ∆ 2 /4, min 1≤τ ≤1/v [τ 2 /λ − ln(vτ )], 1/(λv 2 )}, to make the rhs of Eq. (S6) decay in r as fast as possible. For ∆ ≤ v we choose λ = 2/(∆v) and the maximum is at τ = 1/v, while for ∆ > v, we choose λ to be the solution to the equation λ∆ 2 /2 + ln λv 2 /2 = 1, and the maximum occurs at τ = λ/2. In the end we arrive at Eq. (S3). This finishes the proof of our theorem.
Inserting Eq. (S3) into Eq. (S2), we get where ξ is chosen to be slightly larger than the ξ in Eq. (S3) to compensate the √ d XY prefactor, and the constant C is adjusted accordingly. This proves Eq. (7) in the main text.
Proof of Eq. (9): bound for dynamics initiated from a correlated initial state We assume 2vt < L, since these are the only times we will apply the FSE bounds; at longer times the error bounds become too large to be useful. We have whereÔ l is a unit norm operator that only acts nontrivially in X l , and x , x are the floor and ceiling of x, respectively. The first term in the third line comes from the trivial bound Ŝ l (t) −Ŝ l−1 (t) ≤ 2 Ŝ , while the second term comes from the LR bound Ŝ l (t) −Ŝ l−1 (t) ≤ C(2vt/l) ηl/2 . In the fourth line we use the inequality (2vt/l) ηl/2 ≤ e η(vt−l/2) to facilitate the calculation, and we sum the geometric series in the last two lines, with constants C 1 , C 2 depending at most weakly on t, L. Combined with the second term of Eq. (8), we obtain Eq. (9). One can explicitly compute as needed the constants appearing in this bound for a given model and initial stateρ.

Comparison of FSE bounds to perturbation theory at early times
We can gain insights into the accuracy of our error bounds by comparing them to the true FSE at lowest order in time. We can do this analytically using the Taylor series expansion of the true FSE. Recall from Eq. (2) that FSE is defined as whereĤ L is the L-site finite Hamiltonian and ψ L is the initial product state restricted to this finite chain. The Taylor series expansion of the time evolved operators can be expressed as sums of Lie clusters (nested commutators) involvingŜ and terms in the Hamiltonian [see, e.g. Eq. (S11)]. Many small clusters appear in both e iĤ L tŜ e −iĤ L t ψ L and e iĤtŜ e −iĤt ψ and therefore cancel each other. Only those clusters whose spatial length is at least half of the system size may contribute to FSE. In the following we treat OBC and PBC separately, starting with OBC.
In OBC, FSE is due to those nested commutators in e −iĤtŜ e −iĤt in which boundary terms appear at least once, since only such clusters are not canceled by any cluster in e −iĤ L tŜ e −iĤ L t . The lowest order cluster containing at least a boundary term is proportional to Ĉ (Ŝ,V j ) ψ t D(Ŝ,Vj ) /D(Ŝ,V j )!, whereV j the boundary term closest toŜ, and C(Ŝ,V j ) is the shortest nested commutator joiningŜ andV j . Assuming that the expectation value Ĉ (Ŝ,V j ) ψ does not vanish (which is true for most initial product states |ψ ), the true FSE has the same t dependence as the bound Eq. (5).
In PBC there is a qualitative difference. Any term in the Taylor expansion of e iĤtŜ e −iĤt ψ in Eq. (S9) whose spatial span is smaller than L do not contribute to FSE, because they cancel the corresponding terms in e iĤ L tŜ e −iĤ L t ψ L due to translation invariance and the product nature of the initial state, as shown in Fig. S1. Only those terms in e iĤtŜ e −iĤt ψ which are too long to be embeddable into an L-site periodic system can contribute. More precisely, in PBC, the rhs of Eq. (S9) is contributed by the following two classes of terms: (1) terms in e iĤtŜ e −iĤt ψ whose spatial span is larger than L, i.e. terms that are too long to be embeddable in the L-site PBC chain; (2) terms in e iĤ L tŜ e −iĤ L t ψ L which wrap around the whole periodic system (wraps around the torus, or the circle in d = 1). Such terms do not generally cancel any terms in e iĤtŜ e −iĤt ψ .
The leading order term in both classes are proportional to t L , where L is the number of Hamiltonian terms in the smallest L-site-unembeddable cluster starting fromŜ. This t L behavior is the same as the improved PBC bound in Eq. (6), and the exponent L is roughly twice as large as the OBC exponent D(Ŝ,V j ). In the next section, we give a rigorous bound for the sum of all terms in each class.
FIG. S1. An example of canceling terms in the Taylor expansions of e iĤtσy i e −iĤt and e iĤ L tσy i e −iĤ L t , in the case of 1D TFIM. Vertices represent termsĥj of the HamiltonianĤ = jĥ j , as well as the observableŜ =σ y i , and edges are drawn between terms that do not commute. The spatial span of the two Lie clusters shown in this figure (cluster of open circles, squares, and triangles) are equal and shorter than the system size. The expectation values of these two terms exactly cancel due to translation invariance and the initial state being a product state.

IMPROVED BOUND FOR PBC
For simplicity, we first focus on the 1D case, and later we show that a bound in higher dimension can be obtained by repeatedly using the 1D bound. Consider a translation invariant quantum system on a 1D periodic lattice with L unit cells, described by the HamiltonianĤ L , and letĤ =Ĥ L→∞ be the Hamiltonian in the thermodynamic limit. We will derive an upper bound for the sum of all terms in e −iĤtŜ e −iĤt that may contribute to the FSE (i.e. whose spatial span is larger than L), using ideas motivated by Ref. [63], and then do the same for e −iĤ L tŜ e −iĤ L t .
Let us focus on the first class of terms, i.e. terms in e iĤtŜ e −iĤt whose spatial span is larger than L, since the second class can be treated in an identical way. We write the Hamiltonian in the thermodynamic limit aŝ whereĥ j denotes a local term inĤ. It is convenient to introduce the notion of the commutativity graph G, defined in Ref. [64]. This is a graph whose vertices j are associated withĥ j and which has edges from j to j if and only if h j andĥ j do not commute. The observableŜ is represented as an external vertex s on G, linked to all the vertices j whoseĥ j do not commute withŜ, as shown in Fig. S1. Now we write down the Taylor expansion of e iĤtŜ e −iĤt . We use bold letters h j = adĥ j to denote the adjoint of the corresponding operatorĥ j , e.g.
(it) n n! j1,...,jn∈G T (s,j1,j2,...,jn)∈Ts where T (s, j 1 , j 2 , . . . , j n ) denotes the causal forest of the sequence (s, j 1 , j 2 , . . . , j n ), as defined in Ref. [63], and T s denotes the set of causal trees starting from the vertex s. In simple terms, we are summing over all the non-vanishing connected Lie clusters on G starting from the vertex s.
Our goal is to upper bound the sum over all the terms in Eq. (S11) whose spatial span is larger than L. For such a Lie cluster, let j i be the first term in the sequence j 1 , . . . , j n such that the spatial span of the subsequence h ji . . . h j1Ŝ is larger than L. This means that the spatial span of h ji−1 . . . h j1Ŝ is less than or equal to L. We call j i the earliest unembeddable vertex (EUV) of the sequence (s, j 1 , j 2 , . . . , j n ). Notice thatĥ ji must act nontrivially on at least two unit cells, because otherwise j i can never be the EUV of any sequence.
The basic idea is to classify the Lie clusters in Eq. (S11) into different families, with each family having the same EUV, and derive an upper bound for the sum over all terms within each family. To this end, let us denote by [e iHt (Ŝ)] k the sum of all the L-site unembeddable terms in the rhs of Eq. (S11) whose EUV is k, i.e.
[e iHt (Ŝ)] k ≡ n≥0, T (s,j1,j2,...,jn)∈Ts, EUV(s,j1,j2,...,jn)=k (it) n n! h jn . . . h j2 h j1Ŝ . (S12) We limit our explicit proof to the nearest neighbor interacting case in which every termĥ j in the Hamiltonian acts non-trivially on at most two neighboring unit cells. The proof for the general case is the same as for the nearestneighbor interacting case, but involves keeping track of more complicated notation. Therefore for the general case we omit the proof and present only the final result. Returning to the nearest-neighbor interacting case, for an arbitrary operatorÂ, let x L A , x R A denote the x-coordinates of the leftmost and rightmost unit cells on which the operatorÂ acts, and let n A ≡ x R A − x L A + 1 denote the spatial span ofÂ (we write x L,R k , n k for the operatorĥ k for simplicity). In the nearest neighboring interacting case, n k = 2 if k is an EUV. Without loss of generality, let us suppose that x L k < x L S (the other case x R k > x R S can be treated similarly). In this case we need to have x R S ≤ x R k + L − 1, because otherwise the cluster is already unembeddable beforeĥ k is attached, contradicting the assumption that k is the EUV. Then we have the following theorem: (2) can be proved by Taylor expanding both sides and explicitly comparing terms, using the same spirit as the proof of Lemma 5 and Lemma 6 in Ref. [63]. Intuitively, e iH k t Ŝ gives the sum of all terms in the Taylor expansion of e iHt Ŝ that act inside the region [x L k + 1, . . . , x L k + L], so (e iH k t Ŝ − e iH k t Ŝ ) gives the sum of all terms in e iHt Ŝ that act inside the region [x L k + 1, . . . , x L k + L] and act non-trivially on x L k + L, since those terms which do not act on x L k + L are canceled by −e iH k t Ŝ . Further, only those terms in (e iH k t Ŝ − e iH k t Ŝ ) that act non-trivially on x R k = x L k + 1 can survive the commutation withĥ k in the rhs of Eq. (S13). In short, the surviving terms in (e iH k t Ŝ − e iH k t Ŝ ) are those terms that span exactly L unit cells [x R k , . . . , x L k + L], which are the terms in e iHt Ŝ that are L-site embeddable but become unembeddable immediately after attaching h k . Therefore the rhs of Eq. (S13) gives the sum of all Lie clusters in e iHt (Ŝ) with k being the EUV, since the action of e iH(t−t ) happens at a time later than t and can not change the earliestness ofĥ k . [The last sentence can be better understood by noticing that i t 0 dt e iH(t−t ) h k e i(H−h k )t Ŝ simply gives the sum of all terms in e iHtŜ in whichĥ k appears at least once, with theĥ k at time t being the earliest appearance. Therefore it is natural to expect that i t 0 dt e iH(t−t ) h k (e iH k t Ŝ − e iH k t Ŝ ) gives a subset of the terms in i t 0 dt e iH(t−t ) h k e i(H−h k )t Ŝ , whose EUV is k.] For more general locally interacting Hamiltonians which involve terms that act non-trivially on more than 2 neighboring unit cells, Thm. (2) is generalized to We can get an upper bound for the operator norm of the rhs of Eq. (S15) using the triangle inequality and unitary invariance of operator norm. The result is (S17) The second class of terms [those arising from e iH L t (Ŝ)] can be bounded in a similar way, and it turns out that the resulting upper bound for the second class is identical to the first class, Eq. (S17). Adding up the two classes and taking n k = 2, we get In the following we give two different approaches to upper bound the commutator in the rhs of Eq. (S17). The first one is based on a combination of the methods in Ref. [63,64]. This method leads to the tightest bound but has a high computational cost (potentially exponential in the system size) except in a few simple cases. The second one is based on the differential equation method in Ref. [64], which is slightly looser than the first one, but is much easier to compute, and can also lead to a simple analytic upper bound such as Eq. (6).

CHEN-LUCAS BOUND FOR LR COMMUTATORS
In this section we present the tightest LR bounds for locally interacting systems, obtained by applying the method in Ref. [63] to the commutativity graph introduced in Ref. [64]. The goal is to upper bound [ĥ j ,Ŝ(t)] which appears in the simple bound in Eq. (4) and [ĥ k , e iH k t (Ŝ) − e iH k t (Ŝ)] which appears in the improved PBC bound in Eq. (S17). We begin with the first one. The bound for the second one is a generalization of the first one and is qualitatively smaller at early times. For simplicity we focus on 1D in this section.
The bound for [ĥ j ,Ŝ(t)] is obtained by generalizing Thm. 4 of Ref. [63] to the commutativity graph G. The result is (we present the time integrated version for simplicity) where P jS is the set of all irreducible paths on G from S to j, n(P ) is the number of vertices in P , and h i = ĥ i . An irreducible path P on graph G is a simple path (a path without repeated vertices) in which any two nonconsecutive vertices in P are not adjacent in G. That is, let P = (i n = j, i n−1 , i n−2 , . . . , i 2 , i 1 = S), then we have Fig. S2 for examples of irreducible paths. The bound for [ĥ j , e iH j t (Ŝ) − e iHj t (Ŝ)] is obtained in a similar way: where the first sum is over all vertices q that appears inĤ j −Ĥ j , Y s,jq is the set of all irreducible Y -shapes on G with root s and end points j, q, n j (Y ) is the number of vertices on the j-branch of Y and similarly for n q (Y ), so that n(Y ) = n j (Y ) + n q (Y ) + n s (Y ) + 1. The definition of an irreducible Y -shape with root s and endpoints j, q generalizes that of an irreducible path: it is a three-branch tree with a branch point y (y may coincide with one of s, j, q) such that the three branches P sy , P jy , P qy are irreducible paths, and any vertex in P sy \{y} is not linked (with respect to G) to any vertex in (P jy ∪ P qy )\{y}. The binomial coefficient in Eq. (S20) arises due to the different ways of relative time ordering the operators on the j-branch and the q-branch of Y .
We can see that at early times, the improved PBC bound in Eq. (S20) grows like t min n(Ys)−1 , where the minimum is taken over all the L-cell-unembeddable Y -shapes with root s (notice that the set of smallest L-cell-unembeddable Y -shapes always contain an irreducible one). Put it another way, the small-t exponent of the improved PBC bound is equal to the minimum number of Hamiltonian terms needed to be attached toŜ to make the cluster unembeddable (or, the number of Hamiltonian terms in the smallest unembeddable Lie cluster containing S), which agrees with perturbation theory (provided that the expectation value in |ψ of the leading term doesn't vanish).
In general, the rhs of Eqs. (S19,S20) can only be calculated numerically. Yet there are a few special cases in which we can obtain simple analytic expressions due to the simple structure of the commutativity graph. In the following we show the bound for 1D TFIM as an example. Similar bounds apply to any 1D model whose commutativity graph is a single chain, another example is the FHM in the large-U limit [64].
Example: 1D TFIM in PBC We write the Hamiltonian aŝ whereẐ j,j+1 ≡σ z jσ z j+1 andX j ≡σ x j . For illustrative purpose, takeŜ =σ x i . The commutativity graph G is simply a 1D ring, as shown in Fig. S1. In this case, there are only two irreducible paths between any two points in G. Inserting Eq. (S19) evaluated for the PBC TFIM into Eq. (4) we obtain where we assume for simplicity that L is an odd integer. The improved PBC error bound for δ σ x i (t) ψ is given by Eq. (S18), which in the current case becomes When |j − i| = L, we have e iHj t (X i ) =X i , so we can simply use Eq. (S19). The result is similar to the first term in Eq. (S22) with substitution L → 2L − 1. Inserting Eq. (S24) along with the |j − i| = L case into Eq. (S23), we get

BOUNDING THE PBC ERROR BOUND BY SOLVING A LINEAR DIFFERENTIAL EQUATION
Apart from a few special cases, the computational complexity of the Chen-Lucas bound, Eq. (S20), grows exponentially with system size, since the number of irreducible paths on G grows exponentially in general. For some models with very complicated G, the computation of the Chen-Lucas bound may take even longer than the quantum dynamics simulation itself. For this reason, in this section we give an alternative method based on Ref. [64], which is slightly looser than the Chen-Lucas bound but whose computational time complexity grows only quadratically with system size. In addition, with some further simplifications this method leads to the simple analytic expression in Eq. (6).
The goal here is to bound the FSE of a dynamics simulation on a periodic cluster of size L 1 ×L 2 ×. . .×L d . The first step is to extend the 1D bound to d dimensions. To this end, denote by C j a cluster of size L 1 ×L 2 ×. . .×L j ×∞×. . .×∞, that is, for i ≤ j, the i-th direction is periodic with size L i , while for i > j the i-th direction is infinite. Then we have TakingÂ =Ŝ in Eq. (S36) and inserting Eqs. (S35, S40), we get where t1,2 ≡ 0≤t1≤t2≤t3 dt 1 dt 2 and S l ≡ 2 √ h l ( ls ∈ G) Ŝ . Notice that M is symmetric, and so are e M t , M e M t , etc. Inserting Eq. (S41) into Eq. (S39), we obtain where The restriction on the summation over (i, j) follows from the discussion above Thm. 2 and the definition ofĤ i ,Ĥ i in Eq. (S14), which is essentially the requirement that i is the EUV of some Lie cluster starting from the vertex s, and thatĥ j is a term inĤ 2 −Ĥ 1 ≡Ĥ i −Ĥ i . Notice that the set S is translation invariant in directions L 2 , . . . , L d . Namely, if (i, j) ≡ (( r i , α i ), ( r j , α j )) ∈ S, then for any r 1⊥ ⊥x, r 2⊥ ⊥x, we have (i , j ) ≡ (( r i + r 1⊥ , α i ), ( r j + r 2⊥ , α j )) ∈ S. Let T ⊥ denote the group of all lattice translation vectors in directions L 2 , . . . , L d , such that S ∼ = S/T ⊥ × T ⊥ . We can therefore decompose the sum over i, j as (i,j)∈S = (xi,xj )∈S/T ⊥ ( r i⊥ , r j⊥ )∈T ⊥ . We have (we abbreviate We now need to sum over x i , x j satisfying restrictions defined in Eq. (S44). Notice that all such (x i , x j ) satisfy |x i −x j | = L−1. It turns out to be more convenient to extend the summation to all (x i , x j ) satisfying |x i −x j | = L−1. This still gives an upper bound since the rhs of Eq. (S42) is always non-negative. We let x j = x i ± (L − 1) and sum over x i from −∞ to ∞: This is the bound for the j = 1 term in Eq. (S26). Those j > 1 terms in Eq. (S26) can be treated in an almost identical way; the only difference is that since the directions L 1 , . . . , L j−1 are now periodic, the integrals over k in is the LR speed in the p-th direction, and κ p0 is the position of the minimum. On the other hand, using the Taylor expansion of G ij (t) in Eq. (S31), one can show that the upper bound for f (t) given in Eqs. (S39) and (S43) has a Taylor series expansion with the same leading term as the Chen-Lucas bound in Eq. (S20) and with non-negative coefficients, i.e. f (t) = n≥m f n t n where f n ≥ 0, and m = min n p (Y S ) − 1 is the number of Hamiltonian terms of the smallest Y -shape starting from S that is L p -cell unembeddable in the p-th direction. Therefore, Proposition 1 says In summary, we have where C p ≡ C κp0 e κp0 and L p ≡ min n p (Y S ) − 1. In general, the exponent L p is linearly related to L p , i.e. L p = η p L p − µ p,S . This finishes the proof of Eq. (6).