Tightening the Lieb-Robinson bound in locally-interacting systems

The Lieb-Robinson (LR) bound rigorously shows that in quantum systems with short-range interactions, the maximum amount of information that travels beyond an effective"light cone"decays exponentially with distance from the light-cone front, which expands at finite velocity. Despite being a fundamental result, existing bounds are often extremely loose, limiting their applications when quantitative accuracy is necessary. We derive improved LR bounds in systems with finite-range interactions that qualitatively strengthen previous bounds in three ways: (1) At large distance $d_{XY}$, the bound decays like $e^{-d_{XY}\ln d_{XY}}$, faster than the previous exponential decay. (2) At short times, the bound is $O(t^{d_{XY}})$, which is much smaller than previous bounds, which grow proportional to $t$, and is often tight. (3) The LR velocity vastly improves previous bounds. As an example, when applied to the two-dimensional transverse field Ising model~(TFIM) with interaction strength $J$ and transverse field $h$, our bound for the LR velocity is $\min\{4.27\sqrt{Jh},15.1J,12.1h\}$. Compared to the best previously known bound on the LR speed, $16 e J\approx 43.5J$, the new bound is more than 10 times tighter for $J=h$, at least 2.8 times tighter for every value of $J$ and $h$, and remains finite as $J\rightarrow\infty$, while previous bounds diverge. These tighter LR bounds have immediate consequences. For example, we obtain qualitatively tighter bounds for the correlations in gapped ground states. We also expect the bounds to enable new applications where quantitative accuracy is necessary.


I. INTRODUCTION
The Lieb-Robinson (LR) bound [1,2] has revealed how locality shapes and constrains quantum matter. It says that in short-range (exponentially decaying) interacting quantum systems on a lattice, the influence of any local perturbation is restricted to an effective light cone expanding at a finite speed, apart from an exponentially decaying tail outside the light cone. It has direct implications for many body quantum dynamics, as it has been used to understand the timescales necessary to generate correlation and entanglement [3][4][5][6], to transfer information in a quantum chanel [7], and to prethermalize a system [8][9][10][11]. It has also found applications to equilibrium properties, for example in proving that correlations decay exponentially in a gapped ground state (the exponential clustering theorem [12][13][14][15][16]), that systems with halfinteger spin per unit cell are necessarily gapless (the Lieb-Schultz-Mattis-Hastings theorem [17][18][19]), the area law for entanglement entropy in ground states of gapped systems [20,21], the stability of topological order [22], and the formal proof for the quantization of the Hall conductance as integer multiples of e 2 /h [23]. It is also the theoretical basis for analyzing the efficiency of some numerical algorithms [24][25][26][27][28][29]. It has been generalized to systems with long-range power-law decaying interaction [13,[29][30][31][32][33][34][35][36][37][38], to n-partite connected correlation functions [39], and to dynamics with Markovian dissipation [40][41][42] as well. Furthermore, recent developments in experimental technology of atomic, molecular and optical physics have enabled experimentalists to observe the effective light-cone and measure the LR speed [31,43,44].
Despite the LR bounds' fundamental consequences and efforts to tighten them [13,14,41,[45][46][47], the bounds are often extremely loose (except some special bounds on free particle systems [48][49][50][51][52]), severely limiting their applicability to quantitative problems. There are three qualitative ways in which these bounds are loose. To see these, note that in locally interacting systems, current LR bounds are typically of the form whereÂ X ,B Y are local operators supported on regions X, Y , respectively,Â X (t) = e iĤtÂ X e −iĤt , s XY is the minimal graph-theoretical distance between points in X and Y , v and ξ are constants that only depend on the HamiltonianĤ (and the lattice structure). The first way that this bound is loose is that at short time t → 0, perturbative arguments show that the left-hand side (LHS) of Eq. (1) grows like O(t d XY ), where the exponent d XY grows linearly with s XY when s XY is large, much slower than the O(t) in the RHS of Eq. (1). Second, at large distance s XY → ∞, the LHS actually decays faster than exponential, as already observed by Ref. [2] (although incorrectly claimed to decay like a Gaussian there). A simple example is non-interacting systems, where the LHS decays at large distance like O(e −s XY ln s XY ). Third, the LR speeds v appearing in existing LR bounds are typically much larger than the actual speed of information propagation. For example, in 1D TFIM at critical point, the previous best LR speed is 8eJ ≈ 21.7J while the actual speed determined from exact solution is 2J. Often, the bounds on v are even looser: again in the TFIM, for large J, previous bounds grow ∝ J, while the true LR speed remains finite.
In this paper, we will prove new LR bounds in locally interacting systems (i.e. systems with finite-range interactions) that replace the RHS of Eq. (1) by a new analytic expression, which qualitatively tightens the bounds in all three aforementioned aspects, in some cases to the tightest possible. To achieve this, we introduce and employ a method different from the approaches in prior LR bounds. Typically, the LR bound is proved by first applying the Heisenberg equation and triangle inequality to obtain an integral inequality for the norm of unequal time commutator, and then using the integral inequality iteratively [13,14,45]. One usually get a summation over products of operators along various paths, and then upper bounds the summation by a suitable analytic function. Our approach in this paper is closer to the original method of Lieb and Robinson [1]. The key idea is that in locally interacting systems, we can obtain a system of linear integral inequalities which only involve the unequal time commutator of the Hamiltonian terms. We can then find a system of differential equations (with suitable initial conditions) whose solution gives an upper bound for any solution to the integral inequalities.
This approach leads directly to the tightest bound we obtain, which is given by the solution to a system of coupled linear differential equation with a number of variables proportional to the size of the system, Eq. (10). If the quantitative tightness of the bound is the sole criterion for applying the LR bound, one should solve these either analytically, or numerically if necessary, which can be done quickly and straightforwardly for systems with many thousands of sites. In translationally invariant systems, the solution to the system of linear differential equations can be reduced to an integral by Fourier transformation. For example, in the d-dimensional TFIM, our where ω(k) corresponds to eigenvalue of the linear operator encoding the differential equations, which can be found by diagonalizing a (d + 1) × (d + 1) matrix for each k. The bound for other operators takes a similar form. Based on the Fourier integral representation, we provide a general method to analyze the asymptotic behaviors of the solutions-which shows that the small t and large r behavior is qualitatively tighter than prior bounds, and are often tightest possible-and derive an explicit formula to extract the LR speed. Our bound for the LR speed vastly improves previous ones, as summarized in Table I [ 53]. For example, they are parametrically tighter for the TFIM at large J/h or in large dimension, and even for J/h = 1 and d = 2, they represent more than a 10-fold improvement. Although expressions like Eq. (2) provide the tightest bounds, we can find simple analytic bounds while compromising the tightness only marginally. We will show that We note that this bound applies even to systems lacking translational invariance (under mild realistic requirements on the Hamiltonian to be introduced later). Here d XY is the distance between operatorsÂ X (t),B Y (0) on the commutativity graph to be introduced later (d XY grows linearly with the distance between the two operators in real space), C is a constant, and u is related to the LR speed given in Table I (their precise relation is  given in Sec. IV and Sec. V). This bound also tightens previous ones in all the three aspects mentioned above, and is tightest possible in most of these aspects. For example, the small-t exponent is the same as the exact one obtained by perturbative arguments and the large-x exponent is saturated by some free particle systems.  The tighter LR bounds also improve several important results that rely on it. We mention two here. The first, which we explicitly derive in this paper, is that our LR bounds give rise to a tighter bound for the ground state correlation length in systems with a spectral gap. For example, when applied to the 1D TFIM, our new bound gives ξ −1 ≥ ln J h − 1.3 in the limit J/h → ∞, in agreement with the exact solution, while the previous best bound approaches a constant. As a second example, the improvement of the bound also enables new applications where quantitative accuracy is the key to get helpful results, for example in upper bounding the error of several numerical algorithms [24-26, 28, 29]. Previous LR bounds would give extremely loose numerical error bounds, which renders the error bound practically useless in typical situations where only modest system sizes can be simulated with current computational capabilities. We expect the tighter LR bounds to give practically useful numerical error bounds for reasonably small system sizes.
Our paper is organized as follows. In Sec. II we outline the general method to upper bound the unequal time commutator by the solution of a set of first order linear differential equations. In Sec. III we specialize to trans-lation invariant systems, where one can obtain Fourier integral solutions to the differential equations and obtain the bound in Eq. (2). We also derive an explicit formula for the LR speed by studying the analytic properties of the Fourier integrals. In Sec. IV we derive a power series solution to the aforementioned differential equations in an arbitrary graph, from which we prove the LR bound in Eq. (3). In Sec. V we give examples of our new LR bounds in some classic models (TFIM, FH). In Sec. VI we introduce a special treatment that can be used to obtain a tighter bound when some of the parameters in the Hamiltonian become large. In Sec. VII we show how our LR bounds give a better bound on the ground state correlation length in systems with a spectral gap. Finally in Sec. VIII we summarize the essential tools we introduced, the key results we obtained, and discuss other possible applications of the new bound.

II. GENERAL METHOD: INFORMATION PROPAGATION ON COMMUTATIVITY GRAPHS
There are two major reasons why previous LR bounds are loose. First, the previous methods do not make use of the details of the Hamiltonian. For example, the derivations in Ref. [13] only makes use of how the operator norms of the interaction terms decay with interaction range, completely ignoring all other properties such as commutativity. Therefore, the previous results actually bound the "worst case" Hamiltonian that satisfies the finite range condition. Second, previous derivations use the integral inequalities iteratively, leading to a sum of products of terms along different paths. Since this sum over paths is difficult to calculate analytically, one has to use the triangle inequality many times to get a simple analytic expression, which makes the resulting bound looser.
In the following, we derive our LR bounds using a method that avoids these two difficulties. Motivated by the LR bound of a special model studied in Ref. [45] in the context of topological matter, we make use of the details of the Hamiltonian by carefully exploiting the commutation relations of the terms in the Hamiltonian. To facilitate this procedure, we introduce a graphical tool which we call the commutativity graph that helps visualize the commutation relations of Hamiltonian terms, and will be introduced in Sec. II A. In Sec. II B we will use Heisenberg equation and triangle inequality to upper bound the unequal time commutator and will arrive at a similar integral inequality as previous LR bounds. However, instead of using it iteratively, we will find a set of differential equations whose solution naturally provides an upper bound for any quantities satisfying the integral inequalities. A simple analytic bound for solutions to these differential equations will be discussed in the next section, using a simpler method than the typical sum over paths in the previous approaches, and which produces much tighter bounds.

A. The commutativity graph
We begin by introducing a useful graphical tool, namely the commutativity graph, to represent the Hamiltonian of locally interacting systems. Consider a locallyinteracting quantum system in arbitrary dimension with an arbitrary lattice structure, or on an arbitrary graph. The HamiltonianĤ can in general be written aŝ whereγ j are local Hermitian operators with unit norm γ j = 1, and h j are constant parameters. The commutativity graph G of the HamiltonianĤ is constructed as follows. Each local operatorγ i is represented by a vertex i, the parameter h i is attached to vertex i, and we link two vertices i, j by an edge i, j if and only if the corresponding terms do not commute [γ i ,γ j ] = 0. The resulting graph necessarily reflects the locality, since local operators acting on non-overlapping spatial regions must commute, so there is no edge between their representative vertices. Some simple examples of commutativity graphs are shown in Fig. 1. Notice that the same Hamiltonian may have different decompositions in the form of Eq. (4) and therefore different commutativity graphs, due to the freedom in how to partition terms ofĤ. Nevertheless, for convenience we will simply speak of "the commutativity graph ofĤ" when a certain decomposition is implicitly assumed. We will discuss how to choose the decomposition at the end of Sec. II B. In the following we will derive a Lieb-Robinson bound based on a differential equation on the commmutivity graph of H.
where the dot onγ i (t) means time derivative, and the graph geometry provides a natural way to label the summation.
The general task of this paper is to find an upper bound for unequal time commutators [Â(t),B(0)] of local oper-atorsÂ,B, i.e. operators with finite support. Let us first consider the special case whenÂ =γ i is a term of the Hamiltonian. In this case, we are interested in the unequal time commutatorγ B i (t) = [γ i (t),B(0)], for an arbitrary local operatorB. Using Eq. (5) and the Jacobi identity, we have the evolution equation Now we bound the time dependence of the operator norm ofγ B i (t), the fundamental object controlled by LR bounds.
. Using this and applying the basic inequalities Â +B ≤ Â + B and Â ·B ≤ Â · B , we obtain where in the last line we used γ i (t ) = γ i = 1. Although Eq. (8) bounds the operator norm that is of interest, the right hand side depends on this unknown operator norm. Nevertheless, one can use the generalized Grönwall's inequality [54] is the solution to Eq. (8) with the inequality replaced by equality: whereγ B i (0) = γ B i (0) . Differentiating Eq. (9), one obtains the first order linear differential equatioṅ with initial condition We define the support of operator B on the commutativity graph as S(B) = {i| [γ i ,B] = 0}. Then the initial condition can be simplified asγ B i (0) = 2 B (i ∈ S(B)), where the notation (P ) is defined as (P ) = 1 if P is true and (P ) = 0 if P is false, for an arbitrary statement P .
The differential equation Eq. (10) gives an upper bound for the unequal time commutator of a Hamiltonian term with an arbitrary local operator. To obtain the upper bound ofÂ B (t) = [Â(t),B(0)], whereÂ is an operator not in the Hamiltonian, we can perform similar calculations in Eqs. (5)(6)(7)(8) to get For later convenience, we introduce the Green's function as a useful tool to discuss solutions to Eq. (10). First notice that Eq. (10) can be rewritten in a symmetric forṁ with initial condition G ij (0) = δ ij , so any initial value problem can be obtained by linear combinations of the Green's functions: and therefore the LR bound of [γ i (t),B Y (0)] can be expressed as where Y = S(B) is the support of operatorB on the commutativity graph (from now on we useÂ X to mean that operatorÂ has support X on the commutativity graph). The LR bound for the unequal time commutator between arbitrary local operators [Â X (t),B Y (0)] can be obtained by inserting Eq. (16) into Eq. (12): Eqs. (14) and (17) [or Eqs. (10) and (12)] are the main results of this section, which upper bound the unequal time commutators by the solution to a set of linear differential equations. Note that the bound Eq. (17) can in general be efficiently computed. One can obtain G ij (t) by solving Eq. (14), a set of N H linear coupled ordinary differential equations, where N H is the number of terms in the Hamiltonian, which grows linearly with system size in locally-interacting systems. This is a substantial reduction from solving the original many-body problem, whose cost in general grows exponentially with the system size. Furthermore, in translation invariant systems, we can solve these equations by a Fourier transform, which allows us to derive further simplified analytic upper bounds, as shown in the next section.
We end this section by commenting on the issue noted at the end of Sec. II A. As we mentioned there, we can write the same Hamiltonian in the form of Eq. (4) in different ways, leading to different commutativity graphs. The differential equations and the resulting LR bound will be different as well. One may ask how to choose a decomposition of the HamiltonianĤ so that the resulting LR bound is tightest. While there is not a general statement about which decomposition is best, for spin-1/2 and fermionic systems we know a particularly good choice of decomposition that in many cases give us the tightest bound. For such systems, it is always possible to write the Hamiltonian in the form of Eq. (4) [55,56] if we require thatγ 2 i = 1 andγ i ,γ j for i = j either commute or anti-commute. More precisely, the algebraic relations of {γ i } are related to their commutativity graph G as follows: For example, in spin-1/2 systemsγ i can be chosen as products of Pauli matrices, while in fermionic systems they are products of Majorana fermion operators.
More specific examples are given in Sec. V. If the additional constraints Eq. (18) still fail to specify the decomposition uniquely, we would choose [among those satisfying Eq. (18)] a decomposition such that the total number ofγ i terms is minimal. We call it a minimal Clifford decomposition. The LR bound derived from this decomposition is typically tightest since the triangle inequality [A, B] ≤ 2 A B (which is used many times in our derivation) is saturated when {A, B} = 0. One can also show that the resulting LR bound has the tightest small-t exponent in typical models, as discussed in Sec. IV.

C. An alternative treatment for interacting fermions
Our previous analysis applies to the fermionic case as well, since our derivations did not rely on any special property of the Hamiltonians beyond being locally-interacting. However, for the case of interacting fermions in an arbitrary d-dimensional lattice, we can use a slightly different method which sometimes gives us a tighter LR bound.
In the following we consider an arbitrary locallyinteracting fermionic system with up to quartic terms, but the method generalizes straightforwardly to systems with higher order interacting terms. We will use the Majorana representation for convenience, where the Hamiltonian can in general be written aŝ where t ij and U ijkl are totally antisymmetric, and the Majorana operators satisfŷ We first upper bound the norm of the unequal time anticommutatorĉ The first term of Eq. (22) can be expanded into the sum of three terms using {abc, d} = ab{c, d}+a{b, d}c+{a, d}bc. Therefore Eq. (23) is the analog of Eq. (8). Therefore, using the same argument as we used in Sec. II B, we can conclude with initial condition c im (0) = ĉ im (0) = 2δ im . Once we know the LR bound for {ĉ i (t),ĉ j (0)}, we can calculate LR bounds for arbitrary local operators using identities like [ab, cd] = a{b, c}d−ac{b, d}+{a, c}db−c{a, d}b, since an arbitrary local operator can be expressed as products of basic Majorana operators c j . Eq. (24) takes a similar form as Eq. (10) and can be treated similarly. For example, the Green's function method in Sec. II B can also be applied, and the analyses in the following two sections apply to both equations equally well. We currently do not have a general criteria to determine which equation will lead to a tighter LR bound for a specific fermionic system, and the results have to be compared case by case. But compared to LR bounds obtained prior to the current paper, both methods give improved LR speeds, correct power law growth of correlations at short time, and the superexponentially decaying tail. See the application to the FH model in Sec. V B for a detailed comparison.

III. ANALYTIC BOUNDS IN TRANSLATION INVARIANT SYSTEMS
In the last section we upper bounded the unequal time commutator by the solution to an integral inequality [Eq. (17)] with the integrand given by the solution to a set of linear differential equations [Eq. (14)]. In translation invariant systems, we can obtain the solution to these equations in the form of a Fourier integral, and derive a simpler, fully analytic upper bound by methods of complex analysis.
We will start our discussion by considering Eq. (14) in a general periodic lattice structure and obtain a formal solution in terms of a Fourier integral. Suppose there are l vertices in a primitive unit cell, labeled by index α = 1, . . . , l. We use capital I to label unit cells with coordinate labeled by r I [for example, if {α 1 , . . . , α d } is a set of primitive lattice translation vectors, then coordinate r I = (n 1 , n 2 , . . . , n d ), n a ∈ Z, a = 1, . . . , d labels the primitive cell at physical position R I = d a=1 n a α a ], so that every vertex i can be labeled by a pair (I, α). Taking advantage of translation invariance, we can use the Fourier integral representation of the Green's func- [Notice that other discrete symmetries of the lattice (such as inversion, reflection, discrete rotation, etc.) do not play a significant role in this context, so we are using integer coordinates of unit cells and consequently the Fourier momentum simply lies in (−π, π] d rather than the first Brillouin zone in typical solid state contexts.] Inserting Eq. (25) into Eq. (14) we geṫ The formal solution to Eq.
We now derive a simple analytic upper bound for the Fourier integral of Eq. (25) using methods of complex analysis. Since the integrand of Eq. (25) is a complex analytic function of k except at infinity, and since the integrand is periodic in shifts of the real axis by 2π, we can send k → k + i κ without changing the result of the integral, as illustrated in Fig. 2. Therefore we have where r = r I − r J , in the second line we used the triangle inequality for the integral, in the third line we used where ω m (i κ) is the largest eigenvalue of H (i κ) (must be positive since H (i κ) is nonnegative). Therefore we arrive at the important result: While this expression looks quite similar to previous LR bounds with exponential decaying tail, described by Eq. (1), the point here is that we can choose κ = κ( r, t) to depend on r, t to minimize the RHS at each point in space-time. This makes the bound decaying faster than exponential at large distance, as we prove in Sec. IV, and see in examples in Sec. V. We end this subsection by giving an explicit formula for the LR speed that emerges from this bound. The LR speed of the Green's function G ij (t) (considered as a function of r, t) in Eq. (25) is, by definition, the speed at which r must change as a function of t in order for G ij (t)'s magnitude to stay constant. While the Eq. (25) is not simple enough to exactly extract the LR speed, we can get an upper bound for the LR speed from the upper bound of G ij (t) in Eq. (29). If we choose κ = κ 0 ≡ sgn( r)κ 0 [where sgn( r) ≡ (sgn(r 1 ), sgn(r 2 ), . . . , sgn(r d ))], we get the old form of LR bound with exponential decaying tail where | r| = d j=1 r j is the Manhattan distance between the corresponding unit cells. Therefore the LR speed of |G( r, t)| must be upper bounded by an optimal choice of Notice that while κ 0 depends on sgn( r), in typical models where Equation (31) is one of the main results of our paper, for upper bounding the LR speed of a locally interacting many body Hamiltonian. Compared to previous methods for upper bounding v LR , for example the methods given in Ref. [2] or Ref. [3], our formula Eq. (31) is not only simpler but also gives much tighter bounds, as will be shown in specific examples in Sec. V.

IV. A SIMPLE LR BOUND FOR AN ARBITRARY GRAPH
In this section we will prove the simple LR bound in Eq. (3) for an arbitrary locally interacting system. The only assumption we make about the system is that the commutativity graph G has an upper bound on the degree of vertices and the weight H ij of the edges. We will first express the Green's function in an arbitrary graph as a Taylor series of t, illustrate the graph-theoretical meaning of the expansion coefficients, and then use this expansion to prove the LR bound.
The formal solution to Eq. (14) is where G Notice that H ij can be considered as the adjacency matrix of the weighted commutativity graph G with the weight of the edge (i, j) being H ij . In this way, the meaning of G (n) ij is the sum of the weights of all possible paths in G connecting vertices i and j, where the weight of a path is the product of the weights of all its edges. [Note: there can be duplicate edges in a path, and the weight of path p is w p = e∈p w de e where d e is the multiplicity of edge e.] Therefore, we have G (n) ij = 0 for n < d ij where d ij is the graph theoretical distance (length of the shortest path) between i and j. In summary we have We now prove the LR bound in Eq. (3). First let us prove that there exists a velocity u such that for all pairs of vertices i, j we have Inserting Eq. (34), the above inequality is equivalent to Since the H ij are non-negative, the coefficients G (n) ij are non-negative, so the LHS of Eq. (36) is a non-decreasing function of t while the RHS is independent of t. Therefore we only need to prove that the inequality holds for t = d ij /u: which can be simplified to This is exactly the statement that the Green's function has a finite LR velocity u. One can show that Eq. (38) holds extremely generally, requiring only locality and bounded magnitude of the H ij . For translation invariant systems, we have already proved this in the last section [see Eq. (30)], except that distance is measured differently in Eq. (30) and Eq. (38): in the former case we used the Manhattan distance between unit cells while in the latter case we use graphtheoretical distance in the commutativity graph, so u is in general different from v LR . But since c 1 ≤ d ij /r ij ≤ c 2 for some constants c 1 , c 2 (this follows straightforwardly from the fact that each real space unit cell is associated with a finite number of terms inĤ), existence of a finite v LR guarantees existence of u. In the next section we will show in some specific models how to find a tight u. For the general case without translation invariance, existence of a finite u can be proved under the assumption that the weighted graph G has an upper bound on the degree of vertices and the weights of the edges, that is, every vertex has at most λ edges and the weights satisfy H ij ≤ h for all pairs of neighboring vertices i, j. With these assump- for t ≤ d XY /u, where h XY is a geometric factor defined as In the third line of Eq. (40) we used the inequality(1 + 1/d ij ) dij ≤ e, and in the last line we used (ut/x) x ≤ (ut/y) y e y−x for t ≤ y/u ≤ x/u [this inequality holds since the function (ut/x) x e x is decreasing in x for x ≥ ut].
. In summary we have where C = 2e max{1, 2 h XY u } A B for X ∩ Y = ∅ and ∀t ≥ 0. This completes the proof of the LR bound in Eq. (3). [Eq. (42) is actually a slightly tighter version].
The new form of LR bound in Eq. (42) substantially and qualitatively improves previous LR bounds, not only because it has the superexponential decaying tail e −r ln r , but also because it has much tighter small-t exponent, and a much tighter LR speed, as will be seen in examples in the next section. The small-t exponent η( r) of a bound f ( r, t) is defined by its limiting behavior f ( r, t) ∝ t η( r) as t → 0. We now argue that if a minimal Clifford decomposition [as defined below Eq. (18) where adĤ is the adjoint map ofĤ acting on the space of operators defined as adĤ (B) ≡ [Ĥ,B], P ij is the set of shortest paths connecting i and j (note: by our convention a path p ∈ P ij does not contain i and j), the string operatorŜ p along a path p ∈ P ij is defined aŝ S p = k∈p h kγk (with a suitable ordering). Therefore the small-t exponent is η ij = d ij − 1 provided that the sum on the second line of Eq. (43) is nonzero. A sufficient condition for this to be true is if no two string operatorŝ S p andŜ p along different paths p = p are proportional to each other [Notice that each individual termγ iŜpγj is a product of invertible operators (due to the relation γ 2 j = 1) and is therefore always nonzero]. This condition is satisfied by typical models such as the TFIM, FH model, the Heisenberg XYZ model, and we believe that it is automatically satisfied by the minimal Clifford decomposition of the Hamiltonian. Therefore the small-t exponent d XY + 1 of the bound in Eq. (42) is saturated by exact results, since when X = S(γ i ), Y = S(γ j ), we have d ij − 1 = d XY + 1. This has never been achieved in previous LR bounds, where the small-t exponent is typically just η( r) = 1, as in Eq. (1) for example [57].

V. EXAMPLES
In this section, we apply the techniques developed above to calculate the LR velocity for the bound of Eq. (3) explicitly for two models as examples: the d-dimensional TFIM (Sec. V A) and the FH model (Sec. V B). The techniques are general and can be straightforwardly applied to arbitrary models. We will introduce techniques to further tighten the LR velocity in Sec. VI.

A. The d-dimensional TFIM
The Hamiltonian for the d-dimensional hypercubic lattice Ising model with a transverse field iŝ withγ r,j =σ z rσ z r+êj whereê j is the unit vector in the j th direction, andγ r,0 =σ x r . We will assume J ≥ 0, h ≥ 0 (the resulting LR bound only depends on |J| and |h|). We have written the Hamiltonian in this form so that Eq. (18) is satisfied. It is easily verified that the square of each term is equal to unityγ 2 r,α = 1, α = 0, 1, 2, . . . , d, and any two terms either commute or anticommute. Fig. 3 shows the commutativity graph of this model for the d = 2 case.
We can use the general Fourier integral representation of the Green's function in Eq. (25) where H ( k) , defined in Eq. (27), is in this case a (d + 1) × (d + 1) matrix where 1 is the identity matrix and with cos k ≡ d j=1 cos k j . The eigenvalues of H ( k) are 0 and ω ± ( k) = ±ω( k) = ± 8Jh(d + cos k), from which we can compute the upper bound for the LR speed given by Eq. (31): where the constant X y is defined as the solution to the equation x arcsinh(x) = √ x 2 + 1 + y, and X 0 ≈ 1.50888. The previous best bound for v Ising is 8edJ ≈ 21.7dJ obtained from the method in Ref. [2], so our bound in Eq. (49) is a major, parametric improvement when J or d is large. Even for J = h, d = 2, Eq. (49) represents a 10fold improvement over the previous bound. See Table I for a summary of the comparison. For the case when h is large, an extension of our method will be used to derive a bound which also improves the previous results there as well, as will be shown in Sec. VI.
The LR bound for arbitrary local operators can be obtained from Eqs. (16) and (12). For example, for Finally, [σ x r (t),σ z 0 (0)] also has a simple LR bound in the form of Eq. (42). Notice that in this case d XY is the distance between points X = r and Y = 0 on the commutativity graph, which is exactly 2| r|. Therefore, the parameter u needed to satisfy Eq.
As another example, consider the LR bound for [σ x r (t),σ x 0 (0)] . To get the tightest coefficient we can use Eq. (35) and Eq. (16), where in this case i is the red triangle at r while Y includes the four blue circles near 0, see Fig. 3. We obtain B. The FH model We next take the FH model as an example of interacting fermions, using the method in Sec. II C. The Hamiltonian in a general lattice iŝ (53) In the following, we limit our discussion to a bipartite lattice. We first split the fermion creation and annihilation operators to Majorana operatorŝ where E is the set of even sites, and the Majorana operators satisfyĉ † j,α =ĉ j,α , {ĉ j,α ,ĉ l,β } = 2δ jl δ αβ .
The The norm of the operatorĉ αγ jk (t) ≡ {ĉ j,α (t),ĉ k,γ (0)} is upper bounded by the solution to Eq. (24), which in the current case becomeṡ with initial condition c αγ jk (0) = ĉ αγ jk (0) = 2δ jk δ αγ . Eq. (57) can be easily solved by a Fourier transform. In d-dimensional hypercubic lattice, we obtain the LR bound where I α (t) is the modified Bessel function of the first kind, and I α (t) ≡ dJ. At small U/J, this is 3.02dJ, significantly improving the previous best bound 16edJ ≈ 43.5dJ [2] (improved bounds at larger U are given at the end of this section and in Sec. VI). See Table I for a summary of the comparison.
In the derivations above we did not use the commutativity graph method, but the method in Sec. IV can still be applied to derive the simple LR bound, with d jl = | x j − x l |, and u = v FH,1 , and we have for ∀j, l, α, β, t.
We now compare the bound in Eq. (59) with the exact unequal time anti-commutator at the non-interacting point U = 0 in a d-dimensional square lattice, to investigate how tight the bound is. The Hamiltonian is diagonalized asĤ where ω( k) = 2J cos k. Thereforê where J n (t) is the Bessel function. Comparing this with our bound in Eq. (58), we see that at the non-interacting point U = 0, our LR bound just replaces the Bessel function J xj − x l (2Jt) of the exact solution by the modified Bessel function I xj − x l (2Jt). Since both J x (t) and I x (t) have the asymptotic form 1 x! ( t 2 ) x when x/t is large, we conclude that our LR bound Eq. (58) has the tightest large-x and small-t exponent at the non-interacting point, where the large-x exponent ζ of a bound f (x, t) is defined as sup{ζ| lim x→∞ e ζx ln x f (x, t) exists, ∀t > 0}, [Our bounds Eqs. (58) and (59)  As a comparison, we can also use the commutativity graph method to derive a different LR bound for the FH model which is tighter than Eq. (58) when U/J > 3.135. The commutativity graph of the Hamiltonian in Eq. (56) is shown in Fig. 4, from which we can calculate the eigenfrequencies ω ± (k) = J cos k ± J 2 cos 2 k + 4U J(1 + cos k). (63) and the LR speed is again computed from Eq. The comparison between the LR speeds calculated from different methods is shown in Fig. 5. Both methods introduced so far in this section substantially improve the previous bound for small U , providing better results up to U/J = 21.07. Results that also improve the previous bound at large U are derived in Sec. VI.

VI. IMPROVING THE LIEB-ROBINSON VELOCITY BY ELIMINATING UNPHYSICAL DIVERGENCES
Examining the large-J or large-h limits of v Ising in Eq. (49) reveals that the LR velocity diverges. In reality, however, the velocity of information propagation is expected to be finite. The divergence of the LR velocity is therefore not physical. For example, in the large-J limit of the Ising model, the velocity is expected to be proportional to h rather than the √ Jh of Eq. (49). Such discrepancies often appear when a Hamiltonian is a sum of terms, and one set of the terms is multiplied by coefficients that become very large. Unfortunately, this dis- FIG. 6. The graph G F is constructed from the commutativity graph G of the 2D TFIM in Fig. 3 by removing large parameters, by the construction in Sec. VI. Left: G F for F being the set of vertices representingσ z iσ z j terms, i.e., removing parameter J; Right: G F for F being the set of vertices representingσ x i terms, i.e., removing parameter h.
crepancy renders the LR bound infinitely weak in these limits, a limitation that has plagued prior LR bounds. We introduce a method, extending the results above, that leads to a finite and physically reasonable LR velocity. The basic idea is to derive LR bounds from coupled differential equations on a new subset of the commutativity graph that has removed the operators associated with large coefficients. This is accomplished by an appropriate unitary transform of the remaining operators. Although similar methods were employed in previous works [13,33,36], they were only used to remove single site terms, while our method can be used to remove arbitrary mutually commuting terms in the Hamiltonian.
Consider the commutativity graph G = (V G , E G ) of a locally-interacting spin Hamiltonian defined in Eq. (4), where V G is the set of vertices and E G is the set of edges. Let F ⊂ V G denote a subset of disjoint vertices of G, such that {h j |j ∈ F } is the set of large parameters we want to get rid of. The fact that elements of F are disjoint vertices in G means that the corresponding operators {γ j |j ∈ F } mutually commute. We now construct a new graph G F = (V , E ) from G in the following way: (1) the vertices V = V G \F of G F are obtained from V G by removing elements of F ; (2) E is defined such that for ∀i, j ∈ V , ij ∈ E if and only if either ij ∈ E G or there exists an l ∈ F such that il ∈ E G , lj ∈ E G . In the left (right) panel of Fig. 6 we show the graph G F after removing parameter J (h), respectively.
Let us denote byÎ the sum of all Hamiltonian terms corresponding to vertices in F (i.e. the unwanted terms): The important step here is to consider the evolution of the operatorγ i (for i ∈ G F ) in a way that the unwanted termÎ is "rotated away" by a unitary transformation: Notice thatγ i (t, t) =γ i (t), the operator in Heisenberg picture. Taking derivative with respect to s, we have where the sum is over all terms inĤ −Î that do not commute withγ i (0, t − s). Eq. (67) should be considered as the analog of Eq. (5). Now the same derivations as in Eqs. (6)(7)(8) In this case the generalized Grönwall's inequality [54] indicates that γ B i (t) is upper bounded by the solution γ B i (t) to the differential equatioṅ with initial value γ B i (0) = 2 B . Eq. (69) is the same as Eq. (10) with G replaced by G F .
In summary, when the parameters of a set of disjoint vertices of G becomes very large, we can get a better bound for the speed of information propagation by solving the same differential equations on the new graph G F , which has the large parameters removed, and links added when two vertices of G F are both connected to a (removed) vertex in the original graph. Then all the methods we introduced in Sec. III and Sec. IV apply equally well. The simple form of the bound in Eq. (3) is still true provided that d XY and u are understood as the distance and LR speed on the new graph G F , respectively [58].
As an example, in the case of the TFIM, the left (right) panel of Fig. 6 shows the graph G F constructed from G after eliminating the parameter J (the parameter h) for the d = 2 case, and the resulting upper bound for the LR speed is 4X 0 dh (4X d−1 d dJ). Therefore, combined with Eq. (49), we have the upper bound for the LR speed for the d-dimensional TFIM. This is a major improvement compared to the previous best bound v Ising ≤ 8edJ (calculated from the method in Ref. [2]), as we not only removed unphysical divergence of v Ising in the large-J limit, but also tightened it by a factor of at least 2.4 for all values of d, J and h.
As another example, we consider the FH model. Here, the tightest large-U bound is calculated from the commutativity graph G F which is constructed from G in Fig. 4 by first merging the four blue vertices on every column (i.e., grouping α=1,2,3,4 iĉ j,αĉl,α on each link jl into a singleγ i term) and then eliminating the parameter U according to the prescription of this section. The result for the LR speed is v FH ≤ 8X 0 J. Combined with the results of Sec. V B, we have the upper bound for the 1D FH model, which is also a significant improvement of the previous best bound v FH ≤ 16eJ, see Fig. 5 for the comparison.

VII. CONSEQUENCES OF THE NEW LR BOUNDS ON GROUND STATE CORRELATION DECAY
Our tighter LR bounds can immediately translate into improved bounds in all areas that LR bounds are used, such as bounding the timescale for dynamical processes [3,6,11], studying equilibrium properties like exponential decay of connected correlation functions [12][13][14][15][16], quantifying the entanglement area law [20], and providing error bounds on numerical algorithms [24-26, 28, 29]. In this section we demonstrate this with one of the most important consequences of LR bounds, establishing the decay of ground state correlations in gapped, locally interacting systems [12][13][14][15].
We will achieve this goal by directly replacing the previous LR bound by our new bound Eq. (3) in the proof of Ref. [13], and obtain a new upper bound on the correlation length ξ. The new bound is not only a quantitative improvement, but when applied to TFIM, it is qualitatively tighter in the limit J/h → ∞, where the bound of Ref. [13] approaches a constant while our new bound scales as ξ −1 ≈ ln J h . See Table II for comparison. Let us start with the following well-established inequality which relates the ground state correlation to the unequal time commutator [2,13] where c 0 is a constant independent of α and t, and we assumed G|Â X |G = G|B Y |G = 0 for simplicity. We will split the integral over t into two regions which we will treat differently: for t < r/u we apply our new LR bound in Eq. (3), while for t ≥ r/u we use the trivial where r = d XY , C is a constant defined in Eq. (42), and in the second step we applied the inequality e −αt 2 ≤ e −α(r/u) 2 +2α(r/u)t for t ≥ r/u. If we choose α −1 = λr, then all three terms above decay exponentially with distance r, and the first term is upper bounded by c πu max 0≤τ ≤1/u e −rτ 2 /λ u r τ r−1 . Choosing λ to maximize the smallest decay coefficient of the three terms, the lower bound for ξ −1 is where W (x) is the Lambert W function (the product logarithm) defined as the solution to the equation W e W = x.
We applied this bound to the 1D TFIM, and in Table II we give a comparison between our bound with the previous best bound given in Ref. [13], and also with the exact solution given in Ref. [59] [60]. As we see, in this case our new bound is always tighter than the old one. This is an especially strong improvement in the limit J h where the old bound gives a constant, while the bound in Eq. (74) gives ξ ∝ log(J/h), in agreement with the exact solution.

VIII. CONCLUSIONS
We introduced a method to obtain tighter LR bounds by taking into account the details of the Hamiltonian, and used this method to obtain new analytic forms of the LR bounds which tighten previous ones in three significant ways. The bounds established in this paper have superexponential decaying tails e −d XY ln d XY , where d XY is the separation between the observables on the commutivity graph, instead of previous bounds' exponentials e −µd XY . They have a much tighter, and often the tightest possible, scaling at short times for all distances, namely O(t d XY ), in contrast to previous bounds' O(t) scaling. They also have much tighter LR speeds, often parametrically so. For example, in the TFIM for J h, the bounds of this paper result in v LR ∝ h instead of v LR ∝ J.
The method introduced in this paper follows a simple prescription. To apply this method, one has to (1) Write the Hamiltonian in a suitable decomposition, ideally in a form where any two terms either commute or anticommute. (2) Draw the commutativity graph introduced in Sec. II A and write down the differential equation Eq. (10) whose solution gives LR bounds for the corresponding unequal time commutators; (3) Solve the differential equation. Depending on one's specific needs, we have discussed different methods to do the last step. For applications of the LR bounds in quantitative estimates, one can solve the linear differential equations numerically. In an infinite system with translation invariance, one can write the solution in a Fourier integral form, which can also be computed efficiently, and in this case we derived a simple formula that upper bounds the LR speed. For general systems lacking translation invariance, we have derived a bound of the general form , which is slightly looser than direct numerical solution to the differential equations, but is simple, much tighter than previous bounds, and retains all of the qualitative improvements offered here, namely the small-t behavior, the superexponential decaying tail, and tighter LR speeds. We have illustrated our methods in the examples of TFIM and FH model in Sec. V, where the resulting LR bounds demonstrate significant improvements over the previous best bounds.
Our tightening of the LR bounds may have profound consequences in the study of both equilibrium and dynamical properties of quantum many body systems where the LR bounds have been applied. As a first demonstration, we showed how our new bound leads to a tighter bound on the ground state correlation length in gapped systems, which qualitatively improves previous best bounds in the limit J/h → ∞ of the 1D TFIM, and the asymptotic behavior of our bound agrees with the exact solution. Our method to obtain much tighter LR speeds will also be useful in the quantum information context, where efforts have been made to design protocols for faster quantum state transfer [7], or schemes that create certain entangled states [61][62][63]. Since the LR bounds can be used to bound the time needed for these dynamical processes [3], a sufficiently tight bound not only serves as a criterion to evaluate the performance of a protocol, but also sets a theoretical upper limit on any protocols based on a given physical platform. While previous applications of the LR bounds are mostly limited to analytic proofs, the significant quantitative tightening of the bound shown in this paper will enable numerical applications as well. One such case is to use the LR bound to upper bound the error of a local observable simulated with a finite system. Furthermore, although we only focused on locally-interacting systems in this paper, we expect that our methods can be generalized to tighten LR bounds in systems with power-law decaying few-body interactions.