Strictly linear light cones in long-range interacting systems of arbitrary dimensions

In locally interacting quantum many-body systems, a velocity of the information propagation is finitely bounded and the linear light cone can be defined. Outside the light cone, amount of the information propagation rapidly decays with the distance. When systems have long-range interactions, it is highly nontrivial whether such a linear light cone exists or not. We herein consider generic long-range interacting systems with decaying interactions as $R^{-\alpha}$, where $R$ is the distance. We rigorously prove the existence of the linear light cone in the regime of $\alpha>2D+1$ with the spatial dimension $D$. Our result is expressed as the Lieb-Robinson bound in the form of $\|[O_i(t), O_j]\| \lesssim t^{2D+1} (R-\bar{v} t)^{-\alpha}$ with $\bar{v}=\mathcal{O}(1)$ for arbitrary two operators $O_i$ and $O_j$ which are separated by a distance $R$.

Introduction.-In deep understanding of many-body physics, we necessarily encounter the question on how fast information propagates in the dynamics. In relativistic systems, the information propagation is absolutely prohibited outside the light cone, while in nonrelativistic systems the rigorous light cone cannot be defined. Even in the case, Lieb and Robinson proved in 1972 [1] that an effective light cone can be defined, outside which the amount of the information propagation exponentially decreases with the distance. Here, the effective light cone is characterized by the so-called Lieb-Robinson velocity.
In case of short-range interacting spin systems, the effective light cone is characterized by a finite velocity, and the information propagation is restricted inside of the linear light cone. However, when we consider longrange interacting systems, existence of a linear light cone is quite subtle, since long-range integration induces a fast propagation of information. We here mean by the long-range interaction that interaction strength between separated sites shows a power-law decay as R −α (α > 0) with R the distance. The recent experiments have realized long-range interacting systems with various values of α [71][72][73][74][75][76][77][78][79][80], and exploring universal aspects of the long-range interacting systems attract more and more attention [81][82][83][84][85][86][87]. Depending on the exponent α, both of the linear and the non-linear light cones can appear. From this background, it has been one of the most important and intriguing open problems to clarify whether linear right cone can exist in long-range interacting systems, and the general criterion for it.
In this letter, we rigorously prove the linear light cone in generic long-range interacting systems under the condition of α > 2D + 1. As a related work, in one-dimensional two-body interacting system, the longrange Lieb-Robinson bound has been proved very recently in the form of [ t/R for α > 3 [108], which implies a non-trivial upper bound up to the distance R = O(t). In our analyses, the Hamiltonian is not restricted to few-body interactions, and applicable to arbitrary spatial dimensions. Our Lieb-Robinson bound is given in the stronger form of [O i (t), O j ] t 2D+1 (R −vt) −α . However, only the above commutator relation is not sufficient to upper-bound the whole amount of information propagation outside the light cone. In order to give the linear light cone in strict sense (see (5) and Fig. 1), we also prove that error of the local approximation of O i (t) decays as t D+1 R −α+D outside the light cone (see (7) below). Our result can improve various existing analyses which depend on the polynomial light cone [49, [105][106][107].
Setup and main results.-We consider a quantum many-body system on n sites, where each of the sites sits on a vertex of the D-dimensional graph with Λ the total set (|Λ| = n). We assume that a finite dimen-  . 1. (color online) When we consider a time-evolution Oi(t) of a local operator Oi, the quasi-locality of the interaction ensures that Oi(t) is well-approximated by an operator which is defined on a ball region i [r] having distance at most r from i. We say that the dynamics defines the linear light cone if we can achieve an arbitrary approximation error for r = O(t) as in (5).
sional Hilbert space is assigned to each of the sites. For arbitrary subsystems X, Y ⊂ Λ, we define d X,Y as the shortest path length on the graph that connects X and Y . If X ∩ Y = ∅, d X,Y = 0. For a subset X ⊆ Λ, we define diam(X) := max i,j∈X d i,j , the cardinality |X| as the number of vertices contained in X, and define the complementary subset of X as X c := Λ \ X.
We consider a Hamiltonian as where each of the interaction terms {h Z } Z⊆Λ acts on the sites in Z ⊆ Λ. Our analysis can be also applied to a time-dependent Hamiltonian, but for the simplicity we consider the time-independent Hamiltonians. Notably, we here do not assume the few-body interaction; that is, |Z| can be arbitrarily large in Eq. (1). The only assumption throughout the discussion is the following power-law decay of the interactions: for an arbitrary site pair of {i, j} ⊂ Λ, where Z:Z {i,j} denotes the summation which picks up all the interaction terms {h Z } Z⊆Λ including the sites i and j, · · · is the operator norm, and g 0 is an O(1) constant with g 0 ≥ 1 We are now interested in the time-evolution by the Hamiltonian H. For the simplicity, we consider an operator O i which is locally defined on the site i and analyze Mainly, we focus on the following two values: and where i[r] denotes the set of sites having distance at most r from the site i. The value (4) characterizes the local approximation of a time-evolved operator O i (t) in the region i[r] (see Fig. 1). We note that the decay of (4) is a stronger notion than that of (3). We here define the linear light cone in the following sense. We say that the Hamiltonian dynamics e −iHt defines a linear light cone iff for an arbitrary error δ ∈ R and a fixed t where v t,δ is a constant non-increasing with respect to |t| (i.e., v ∞,δ = const.). From the definition, the amount of information propagation is smaller than δ outside the region separated by the distance v t,δ |t|.
We here show our main theorems: for ∀j ∈ Λ with d i,j = R. Also, the Lieb-Robinson bound for (4) is given by where R >v|t| is considered and C H , C H andv are constants which only depend on the parameters {D, g 0 , α} and a geometric constant which is defined by the lattice structure. We emphasize that the same upper bound is obtained for generic operators O X and O Y (see the supplementary materials [109]).
Here, the coefficient log 2D (R + 1) appears for a technical reason due to the macroscopic interactions (i.e., h Z with |Z| 1) in Eq. (1). The first inequality (6) is stronger than the second one (7) in the sense that asymptotic decay is as small as O(R −α ). From the inequality (7), the Lieb-Robinson velocity v δ,t is given by Intuitive explanation of α > 2D + 1 and strategy of the proof.-Before showing the proof outline of the main theorem, we intuitively explain why the condition α > 2D + 1 is necessary. To concentrate on the physics, let us fix the time t, and focus on the simplest case where the Hamiltonian is given by H = H 1 +H , where H 1 has only nearest-neighbor interactions and H 2 has a length scale from to 2 as follows: When we consider the unitary operator e −iH1t , the standard Lieb-Robinson bound [1, 14, 15] clearly gives the linear light cone. As for e −iH t , the Lieb-Robinson bound gives e −c(x/ −v t) and v is proportional to the one-site energy: where j∈Λ h i,j ≤ g 0 j: ≤di,j ≤2 d −α i,j is a summation of all the interaction terms which act on a site i. Hence, the linear light cone is obtained for α > D + 1 and t, which gives a necessary condition; indeed, for α ≤ D + 1, there exists a counterexample for the linear light cone [94,95]. Therefore, the product of unitary operator e −iH1t e −iH t retains the linear light cone as long as α > D + 1.
On the other hand, in considering the unitary operator e −i(H1+H )t , we have to decompose the contributions of H 1 and H as where H (H 1 , τ ) := e iH1τ H e −iH1τ and T denotes the time-ordering operator. Because one-site operator spreads up to a distance O(τ ) due to the time-evolution e −iH1τ (Fig. 2), the one-site energy is now given by where we formally denote H (H 1 , τ ) = Z h τ,Z . Therefore, the time-evolution T e . When the length scale is given by O(t), the linear-light cone retains only for α > 2D + 1. In summary, the operator spreading changes the effective one-site energy by t D times, which yields the condition of α > 2D + 1 for the linear light cone.
In our proof, we decompose all the length scales into pieces and consider the multi-unitary decomposition by generalizing Eq. (8). We then obtain the Lieb-Robinson bound for each of the decomposed unitary operators and connect them into a single Lieb-Robinson bound. The technical difficulties lie in that we need to connect infinitely many Lieb-Robinson bounds; in the step by step connections, a simple estimation makes the Lieb-Robinson velocity diverge rapidly, and hence highly refined analyses are required to obtain a finite velocity.
Sketch of the Proof.-We herein show the essential ideas and give the details in the supplementary materials [109]. For the proof, we first decompose the length scale into ≤ t and > t for a fixed t, and consider the following decomposition of the total Hamiltonian into where H ≤ ( ∈ N) is defined as the operator which picks up all the interaction terms with length scale smaller than : In the case where the length is short range or = O(1), the Hamiltonian H ≤ gives the Lieb-Robinson bound with a finite velocity. However, we here consider the case of = t with t depending on the time t. When the length scale is in the middle range between = O(t 0 ) and = ∞, it is no longer trivial whether or not the light-cone for the dynamics by H ≤ t is retained.
In order to obtain the Lieb-Robinson bound for H ≤ with generic , we further decompose it as follows (Fig. 3): We define a set of the length scales In this choice, the q increases by double exponential function with respect to q.
In deriving the Lieb-Robinson bound for e −iH ≤ t , we iteratively take a longer length scales into account. We begin with the unitary operator e −iH1t . As long as 1 is independent of t, the Lieb-Robinson bound for e −iH1t is short range with a velocity of v 1 = O(1). Then, by using this Lieb-Robinson bound, we derive a new Lieb-Robinson bound for e −i(H1+H2)t = e −iH1:2t , where we define H 1:q := q s=1 H s with 1 ≤ q ≤ q * . We repeat this process, extending the length scales as 1 → 2 → · · · → q * = . In each of the steps, based on the Lieb-Robinson bound for e −iH1:q−1t , we update it to the Lieb-Robinson bound for e −iH1: qt We start from the unitary operator e −iH1:2t : The operator spreading by e −iH1t is order of O(v 1 t), and hence as long as ∆t 2 /v 1 , H 2 (H 1 , τ ) (τ ≤ ∆t) has the same length scale as the original one, namely 2 . Thus, we can obtain the Lieb-Robinson bound with a linear light cone for T e − ∆t 0 H2(H1,τ )dτ . For generic time t, we consider U 2,t = U m 2,∆t2 U 2,∆t 2 , with t = m∆t 2 + ∆t 2 , where ∆t 2 ∝ 2 /v 1 and ∆t 2 < ∆t 2 . In this way, the unitary operator e −i(H1+H2)t is decomposed to pieces of short-time unitary operators which have a length scale at most O( 2 ). By appropriately connecting all the Lieb-Robinson bound, this results in a Lieb-Robinson bound with the following form: For general q −1, we define v q−1 as the Lieb-Robinson velocity for e −iH1:q−1t , and analyze the unitary operator U q,t := e −iH1:qt . We decompose it as U q,t = U m q,∆tq U q,∆t q with 3. (color online) Decomposition of the length scale. We first decompose the total Hamiltonian into two regimes, namely H ≤ t and H > t . The Hamiltonian H ≤ t includes the interactions up to the length scale t which depends on the time, and H > t includes the other interactions. In order to obtain the middle-range Lieb-Robinson bound, we further decompose the range [1, t] into q * pieces. We start from the Hamiltonian only including the length scale 1, and we iteratively take into account the longer length scales.
and t = m∆t q + ∆t q , where ∆t q ∝ q /v q−1 and ∆t q < ∆t q . From the choice of ∆t q , we can ensure that H q (H 1:q−1 , τ ) has the same length scale as that of H q , namely q . We then obtain the Lieb-Robinson bound for H 1:q as follows: where v q depends on the v q−1 . In this way, we iteratively estimate the Lieb-Robinson velocity v q by using v q−1 . We can derive the following recursion relation: We give more explicit form in [109]. The length scale q is now lower-bounded by a double exponential function with respect to q. Therefore, lim q→∞ v q converges to a constant v * , we obtain the following "middle-range Lieb-Robinson bound" for e −iH ≤ t : In this way, we can ensure that H ≤ retains the linear light cone for t, while for t we cannot. We choose as t = |t|η withη := 1 − α−2D−1 2(α−D) < 1 and decompose the total time evolution as where We apply the middle-range Lieb-Robinson bound (10) to e −iH ≤ t t and utilize the standard recursion approach [14, 15] to the long-range contribution by U > t . Because of the middle-range Lieb-Robinson bound, we can ensure that . By taking into account the effect, after intricate calculations, we obtain the Lieb-Robinson bound for U > t as We then connect the two Lieb-Robison bound for e −iH ≤ t t (10) and U > t (11) to prove the inequality (3). In order to give the linear light cone in strict sense (see (5) and Fig. 1), we also prove that error of the local approximation of O i (t) decays as t D+1 R −α+D outside the light cone (see (7) below).
Summary and discussion.-In this letter, we rigorously proved the existence of the linear light cone (see (5) for the definition) in generic long-range interacting systems, where the interaction decays as R −α (α > 2D + 1) with respect to a distance R. Our Lieb-Robinson bound in Theorem 3 roughly gives the commutator relation as t 2D+1 (R −vt) −α and rapidly decays beyond r vt. Moreover, the error of the local approximation of O i (t) was estimated Our result is quite general in that only the condition (2) is imposed to the Hamiltonian and no other assumptions such as the few-body interactions are required. Our work may improve many of the previous works which depends on polynomial light-cone [49, [104][105][106][107]: for example, the speed of the scrambling time [51, 52], Haah-Hastings-Kothari-Low algorithm [48,49] in simulating the long-range interacting systems, and so on. Moreover, although we consider the Hamiltonian dynamics throughout the letter, we expect that the same analysis is applicable to generic Markovian dynamics by following Refs. [18,110].
We finally mention open questions. The most important problem is whether or not we can improve the condition α > 2D + 1 for the linear light cone. The improvement seems to be highly challenging since we utilize this condition throughout the analyses. We believe that only the restriction (2) is too loose and additional conditions such as few-body interaction of Hamiltonian are necessary for the improvement of the condition. In tackling the problem, the simplest case that includes only two length scale as in Eq. (8) may be a good starting point. So far, even if we restrict ourselves to two-body-interacting Hamiltonians with only two-length scales, we have not succeeded to obtain the linear light cone for α ≤ 2D + 1.
As another problem, it is also intriguing to improve the polynomial (or superlinear) light cone in the regimes of α ≤ 2D+1. The state of the art analysis [49] gives the polynomial light cone in the form of r = t (α−D)/(α−2D) . We expect that our present analysis should lead to a better light cone in high-dimensional systems.

ACKNOWLEDGMENTS
The work of T. K. was supported by the RIKEN Center for AIP and JSPS KAKENHI Grant No. 18K13475. K.S. was supported by JSPS Grants-in-Aid for Scientific Research (JP16H02211 and JP19H05603).  We here recall the setup. We consider a quantum spin system with n spins, where each of the spin sits on a vertex of the D-dimensional graph with Λ the total spin set (|Λ| = n). We assume that a finite dimensional Hilbert space is assigned to each of the spins. For a partial set X ⊆ Λ, we denote the cardinality, that is, the number of vertices contained in X, by |X| (e.g. X = {i 1 , i 2 , . . . , i |X| }). We also denote the complementary subset of X by X c := Λ \ X.

B. Definition of the long-range interacting Hamiltonian
We consider a Hamiltonian as where each of the interaction terms {h Z } Z⊆Λ acts on the sites in Z ⊆ Λ. Note that we here do not assume the few-body interaction (i.e., |Z| = O(1)). Throughout the paper, for arbitrary operators A and O, we denote for the simplicity of the notation. In order to characterize the long-range interaction of the Hamiltonian, we impose the following assumption for the Hamiltonian: Assumption 1 (Power-law decaying interactions). We assume the power-law decay of the interaction in the following sense: for an arbitrary site pair of {i, j} ⊆ Λ, where · · · denotes the operator norm and g 0 is an O(1) constant with g 0 ≥ 1.
The above condition includes the following assumption: Assumption 2 (Power-law decaying interactions (alternative)). When the Hamiltonian (S.1) satisfies the condition (S.4), the following inequality also holds where g is bounded from above by Our purpose is to estimate the upper bound of where we use the assumption (S.4) in the last inequality. By using the inequality (S.15) with ξ = 1, we obtain Also, for an arbitrary By combining the inequalities (S.8) and (S.9) with (S.7), we have the inequality of We thus obtain Eq. (S.6).

C. Coarse grained set
For a subset X ⊆ Λ, we define X[r] as We also define a coarse grained total set Λ (r) as the minimum subset such that Λ (r) [r] = Λ, namely where Λ (0) = Λ. Similarly, for arbitrary subset X ⊆ Λ, we define X (r) ⊆ Λ (r) as follows: where X (0) = X. From the definition, we notice that We introduce the geometric parameter γ which is determined by the lattice structure. We define γ as a constant of O(1) which satisfies for arbitrary ξ ∈ N, for r ≥ ξ. By using the parameter γ, we obtain for an arbitrary subsystem X ⊆ Λ, where diam(X) = max i,j∈X d i,j . Furthermore, we compare two subset X (ξ) and , and hence where in the second inequality we utilize the third inequality in (S.16).
is defined by extending the original subset X by a distance r. (b) The subset Λ (r) ⊆ Λ is the coarse grained lattice (orange sites) which is defined as the minimum subset such that Λ (r) [r] = Λ. For arbitrary subset X ⊂ Λ, X (r) ⊂ Λ (r) is the coarse grained subset (pink sites) which is defined as the minimum subset such that X (r) [r] ⊇ X. The cardinality of the subset X (r) is roughly (1/r) D times as the original one, namely |X (r) | ≈ (1/r) D |X| as shown in (S.16).

A. Definitions
We first define G(x, t, X , Y)-Lieb-Robinson bound and G(x, t, X )-Lieb-Robinson bound as follows: The definition clearly implies G(x, t, X , Y) ≤ 1 and we assume that the cut off is automatically included in the definition. For example, when we give G(x, t, X , Y) in the form of However, we omit the max(1, · · · ) for the simplicity of the notation. The function G(x, t, X ) is connected to an error in approximating time-evolved operator on a local region. We show the following lemma which has been given in Ref.
where O X = 1 and X[r] was define in Eq. (S.11).

B. Main theorem
We here show our main theorems: Theorem 3 (Lieb-Robinson bound for long-range interacting systems). Let us consider the long-range interacting Hamiltonian H satisfying the assumption 1. For |t| ≥ 1, the Hamiltonian H satisfies the G H (x, t, X , Y)-Lieb-Robinson bound as where C H , C H andv are constants which only depend on the parameters {D, g 0 , α, γ}. Note that we use the notation in Eq. (S.13) for X (v|t|) and Y (v|t|) .
Remark. From the theorem, outside the light-cone of R ≥v|t|, the commutator is as small as O(R −α ). The first inequality (S.23) is stronger than the second one (S.24) when subset Y is small. However, in the case where Y is infinitely large, the inequality (S.23) is useless. For example, in order to obtain the local approximation (S.22) for time-evolved operators, we need the second inequality (S.24). From the inequality (S.24), the light-cone of v δ,t is given by where we use the condition α > 2D + 1.

C. Proof of Theorem 3
Before discussing the long-range interaction, we consider the Lieb-Robinson bound which is obtained from the following Hamiltonian with an interaction truncation: As long as = O(1), the Hamiltonian gives the Lieb-Robinson bound with a finite velocity. However, when the length depends on the time t, it is highly nontrivial whether the unitary operator e −iH ≤ t satisfies the Lieb-Robinson bound with a finite velocity or not. This kind of "middle-range" Lieb-Robinson bound is critical in discussing the long-range Lieb-Robinson bound. We here show the following theorem which we will prove in Sec. IV: Theorem 4 (Middle-range Lieb-Robinson bound). We consider a Hamiltonian H which satisfy Assumption (S.5) with α > 2D + 1. Then, for arbitrary , the Hamiltonian H ≤ satisfies the G ≤ (x, t, X )-Lieb-Robinson bound with where v * is a constant which only depend on the parameters {D, g 0 , α, γ} and From Theorem 4, we ensure that as long as ≤ O(t), the Lieb-Robinson bound retains the linear light cone. We prove the long-range Lieb-Robinson bound by using this middle-range Lieb-Robinson bound. For the purpose, we consider the following decomposition of the total Hamiltonian into where we define H ≤ t by Eq. (S.25). Note that under the condition of α > 2D + 1 we haveη < 1 and the dynamics by H ≤ t retains the linear light cone. We then decompose the unitary operator e −iHt into the form of We have already obtained the Lieb-Robinson bound for e −iH ≤ t t from Theorem 4. Hence, we need to consider the Lieb-Robinson bound for the unitary operator of We can prove the following theorem on the Lieb-Robinson bound for the unitary operator U > t (see Sec. VI for the proof): Theorem 5 (Contribution by the long-range interacting terms). Under the choice of t by Eq. (S.29), the unitary operator (S.31) satisfies the G > t (x, t, X , Y)-Lieb-Robinson bound as whereJ 0 and r 0 are constants which only depend on the parameters {D, g 0 , α, γ}.
By using Theorems 4 and 5, we prove the main theorem as follows. We here estimate the norm of the commutator where we set ξ j * = R/2 and appropriately determine {ξ j } j * −1 j=1 afterwards. From the inequality (S.22) in Lemma 2, we have We thus obtain where we use Theorems 4 and 5 with = t = tη. We, in the following, determine {ξ j } j * −1 j=1 so that where we use the inequality in (S.16). Thus, we need to choose ξ j as with c 1 and c 2 constants which depend on v * ,C 0 andη. By applying the inequality (S.37) to (S.36), we obtain On the other hand, we obtain for ξ j ≤ R/2 By applying the inequality (S.41) to (S.40), we obtain We here define c 3 as where the inequality comes from the condition of |t| ≥ 1. We here chooser 0 t such that for R ≥r 0 v * |t| > 2r 0 v * |t| the following inequality holds Now,r 0 is a constant which only depends onC 0 , r 0 and v * . The above inequality reduces the inequality (S.43) to We note that the subset Y is now restricted by the condition diam(Y ) ≤ 2v * |t|.

A. Relation between local Lieb-Robinson bound and global Lieb-Robinson bound
We first consider the Lieb-Robinson function G(x, t, X , Y) for local subsets X , Y with diam(X), diam(Y ) ≤ ξ (ξ: constant). In our analysis, it is crucial to connect the Lieb-Robinson function for local operators to that for generic operators.
To make the motivation clearer, we consider an connection of the Lieb-Robinson bounds from different timeevolutions e −iH1t and e −iH2t . We are now interested in the commutator of [e iH2t e iH1t O X e −iH1t e −iH2t , O Y ] . The obstacle in the connection comes from the subset dependence of the Lieb-Robinson bound. We often obtain the Lieb-Robinson bounds for U 1 and U 2 in the form of for j = 1, 2 and d X,Y = x, where F j (x, t) reflects the details of the locality of U j . In the connection, we consider the decomposition Then, for arbitrary subsets L ⊆ Λ and X ⊆ Λ with diam(X) ≤ 2ξ, we obtain the Lieb-Robinson function G(x, t, X, L) which is bounded from above by for ∀X ⊆ Λ and ∀Y ⊆ Λ such that diam(Y) ≤ 2ξ, where we assume F(t, X , Y) = F(−t, X , Y) and the function c(x) monotonically increasing with respect to x. Then, for arbitrary two subsets L, L ⊆ Λ, we obtain the Lieb-Robinson function G(x, t, L, L ) which is bounded from above by (S.62)

Proof of Theorem 6
We first prove the inequality (S.54). For the proof, we consider the norm of , the partial trace with respect to an arbitrary subset X 0 ⊆ Λ is given by where we use L (ξ) [ξ] ⊃ L in the last equation. Also, we notice that We thus obtain We obtain similar relations to (S.67) and (S.68) as follows: where we define U L[R] c := ñ s=1 U Xs . By using the inequality (S.54), we have and hence We therefore obtaiñ By using the inequality (S.15), the summation with respect to i is bounded from above by By using the form of where we use ξ ≥ ξ 0 ≥ 1 and the inequality

Proof of Corollary 7
First of all, because the derivation of (S.54) does not rely on diam(X ) and the form of L(x), we obtain the inequality (S.59) in exactly the same way.
Our task is to prove the inequality (S.60). We start from the inequality (S.72): For an arbitrary integer R, we obtain the same inequality as (S.77) as follows: We now calculate the summation with respect to s for L(x) = x −α0 . Assuming R + 1 ≥ ξ(α 0 − D + 4 + r), we obtain By combining the inequalities (S.91) and (S.92), we have which reduces the inequality (S.90) tõ

B. Reformation of the Lieb-Robinson bound
In our analysis, we often utilize the G(x, t, X )-Lieb-Robinson bound with a fixed form of G(x, t, X ). We here prove the following lemma: with C ≥ 1. Then, for |t| ≥ ∆t > 0, the function G(x, t, X ) is bounded from above by (S.97)

Proof of Lemma 8
We first note that as long as |t| ≥ ∆t the inequality G(x, t, X ) ≤ 1 can hold only for x ≥ v∆t. Then, for x ≥ v∆t, we have Second, from the definition (S.13), we have X ⊆ (X (r) )[r] for an arbitrary subset X ⊆ Λ and hence For |t| ≥ ∆t, we obtain which gives the inequality (S.96). By combining the inequalities (S.98) and (S.100), we have This completes the proof.

A. statement
Theorem 4 (Middle-range Lieb-Robinson bound) We consider the Hamiltonian H which satisfies Assumption (S.5) with α > 2D + 1. Then, for arbitrary , the Hamiltonian H ≤ satisfies the G ≤ (x, t, X )-Lieb-Robinson bound with where v * is a constant which only depend on the parameters {D, g 0 , α, γ} and

B. Decomposition of the length scale
We define a set of the length scales { q } q * q=1 as where we choose {η q } q * q=1 appropriately such that { q } q * q=1 become integer and there exists an integer q * ∈ N satisfying q * = (S.111) We define 1 as a constant which only depends on the parameters {D, g 0 , α, γ}. We show the condition for 1 in Sec. IV C. We then define Hamiltonian H q which picks up all the interaction terms with length from q−1 to q : 1, 2, . . . , q * − 1), where we set 0 = 1. From the assumption 2, we obtain From the definition of η (S.106), we have D+1 q

C. Conditions for 1
In the proof, we adopt various kinds of condition for 1 . We choose 1 so that the following conditions for q (q ≥ 1) are satisfied. We note that all the conditions depend only on the parameters {D, g 0 , α, γ}. We summarize all the conditions in the following: 1. Used in the inequality (S.138): whereC 0 is defined in Eq. (S.104).

D. Proof of Theorem 4
We first consider the case of ≤ 1+η 1 . In this case, we cannot define q * such that q * = in Eq. (S.106), but the Hamiltonian H ≤ has a finite length scale which is independent of t. We here consider the Lieb-Robinson bound by a Hamiltonian which has a length scale at most ξ. As long as ξ = O(1), we trivially obtain the Lieb-Robinson velocity of order O(1). However, it is possible that the velocity increases polynomially with ξ. Our purpose here is to prove that the Lieb-Robinson velocity is at most linear to ξ as long as α > 2D + 1. We prove the following lemma: We note that the same inequality holds for time-dependent Hamiltonian.
From Lemma 9 or the inequality (S.126), we obtain whereη is defined in Eq. (S.106). Because 1+η 1 v 0 depends only on the parameters {D, g 0 , α, γ}, we obtain the inequality (S.103) for the case of ≤ 1+η 1 . We next focus on the case of > 1+η 1 , where we can define q * ≥ 2 such that q * = . We here define the Hamiltonian H 1:q as follows: For the proof, we need to estimate the Lieb-Robinson bound for the Hamiltonian H 1:q * . We define that {H 1:q } q * q=1 satisfy the G q (x, t, X )-Lieb-Robinson bound for q = 1, 2, . . . , q * , namely, for arbitrary operators O L and O L . We aim to prove the Lieb-Robinson bound for H 1:q in the form of whereC 0 is defined in Eq. (S.104), and v q is given by as defined in Eq. (S.127). By remembering that the q is lower-bounded by a double exponential function with respect to q as in (S.111), we can ensure and v ∞ depends only on the parameters {D, g 0 , α, γ}. From the inequalities (S.128) and Eq. (S.133), we can choose v * = max( 1+η 1 v 0 , v ∞ ). We thus prove the theorem. In the following, we prove the inequality (S.131) by induction method. We first consider the case of q = 1. From the inequality (S.126), we obtain which clearly reduces to the form of (S.131). We then adopt the assumption for H 1:q−1 (q ≥ 2) with the Lieb-Robinson bound of and derive the Lieb-Robinson bound for H 1:q in the form of (S.131). We first restrict ourselves to and reduce the inequality (S.135) to the form of (S.139) Note that v q−1 ≤ 2v q−1 is a consequence from the condition (S.116). The derivation of the inequality (S.137) is given by as follows: By using Lemma 8 with r = q , p = D − 1, ∆t = ∆t q , v = v q−1 and ξ = q−1 /2, we obtain where we use η ≤ η q ≤η. From the condition (S.116), we have δ q−1 ≤ 1/2, which yields 2(1 − δ q−1 ) ≤ 1, and hence the inequality (S.140) reduce to (S.137). We then consider the unitary operator e −iH1:qt with t ≤ ∆t q . We start from e −iH1:qt = e −iH1:q−1t T e Hq(H1:q−1,τ )dτ by the functionG q (x, t, X ). In order to estimatẽ G q (x, t, X ), we first estimate the quasi-locality of H q (H 1:q−1 , τ ). We prove the following proposition: Proposition 10. Let˜ q be an length scale such that for fixed τ ≤ ∆t q . Then, under the assumption of (S.137), we obtain a decomposition the time-evolved Hamiltonian H q (H 1:q−1 , τ ) in the form of for an arbitrary positive integer s 0 , where h τ,Z[s˜ q ] is an operator which acts on the subset Z[s˜ q ].
Also, we can prove the following proposition for quasi-local operator: Proposition 11. LetH be a Hamiltonian such that For an arbitrary subset X s.t. diam(X) ≤ s 0 ξ, we assume Then, the HamiltonianH satisfies GH (x, t, X )-Lieb-Robinson bound with We note that the same inequality holds for time-dependent Hamiltonian.
From the inequalities (S.138) and τ ≤ ∆t q = q /(12v q−1 ), by choosing˜ q = 5 q /12, we have 2γ 2 e v q−1 τ / q−1 e −˜ q /(2 q−1) ≤ 2γ 2 e 2vq−1( q /(12vq−1))/ q−1 e −(5 q /2)/(12 q−1) = e − ηq q−1 /24+log(2γ 2 ) < 1, (S.151) where in the first inequality we use v q−1 ≤ 2v q−1 which is given in (S.138), and in the second inequality ηq q−1 ≥ η q−1 ≥ log(2γ 2 ) is derived from the condition (S.117). From the choice of˜ q = 5 q /12, we have Hence, in order to apply Proposition 11 to H q (H 1:q−1 , τ ) (τ ≤ ∆t q ), we choose {ξ, µ,g} as We then obtain the Lieb-Robinson functionG q (x, t, X ) for dynamics by the time-evolved Hamiltonian H q (H 1:q−1 , τ ) (0 ≤ τ ≤ t with t ≤ ∆t q ): where we use the condition (S.151) to obtain 1 − 4 q−1 / q ≥ 11/12, or q / q−1 = ηq q−1 ≥ 48. We further simplify the bound as x +ṽ q |t| + (e 2ṽ q |t|) 1/2 · exp log(e 2ṽ q |t|) 8 q x . (S.156) where we use v q−1 ≥ v 1 = 1 v 0 and the inequality (S.115) for g q in the first inequality and the condition (S.118) in the second inequality. This reduces the inequality (S.156) tõ for t ≤ ∆t q , where we use −η q−1 /e 2 exp( −η q−1 /e 2 ) ≤ −η q−1 in deriving the first term. We then obtain the function G q (x, t, X , Y) for the Hamiltonian H 1:q . For the purpose, we first estimate We then obtain the general form by using Theorem 6. For arbitrary ∆t q < ∆t q , we consider for diam(X ), diam(Y) ≤ νξ (ν ≥ 1). Also, let H 0 be a Hamiltonian which satisfy the Lieb-Robinson bound as We then obtain the Lieb- where r * is defined as an integer which satisfies In the following, we prove the form of For the parameters ξ q and R q , we obtain where in the first inequality we use the condition (S.121), and in the second inequality we use e q−1 log(ζ1)/ξq ≤ e which is a direct consequence from the condition (S.119). We prove the inequality (S.167) in the subsection IV E. From the inequality (S.167), we can write G q (x, t m , X , Y) as where we use (m − 1)∆t q ≤ t m and ξ q ≤ q /2 in (S.170) for the inequality and define v q as where we have defined ∆t q = q /(12v q−1 ) in Eq. (S.136). In this way, we obtain the G q (x, t, X , Y)-Lieb-Robinson bound only for local subsets X , Y (diam(X ), diam(Y) ≤ q ). Finally, by using Theorem 6 with ξ = ξ 0 = q /2, we obtain the Lieb-Robinson bound for generic operators as follows: We thus prove the inequality (S.131). This completes the proof of Theorem 4.

E. Proof of the inequality (S.167)
We prove the inequality (S.167) by induction. For m = 0, we have G q (x, t 0 , X , Y) = 0 and the inequality is trivially obtained.
Then, we assume the form of (S.167) for a fixed m and consider the case of m + 1. We here consider the Lieb-Robinson bound for e −iH1:qtm+1 = e −iH1:qtm e −iH1:q−1∆tq T e For the first, we consider the Lieb-Robinson bound for e −iH1:qtm e −iH1:q−1∆tq , which we characterize by From the form of G q−1 (x, t, X ) given in (S.137), we can decompose We notice that for 0 ≤ x 0 ≤ x, where the inequality 1/ q−1 ≥ 2/ξ q comes from the condition (S.120). From Proposition 12, we have the following bound: where ν q and r * q,1 are defined by and We note that ν q ≥ 2 comes from the condition (S.121). Now, r * q,1 is given by We also obtain 1 + 2r * q,1 with ∆t q = q /(12v q−1 ), where we use the condition (S.122), ξ q ν q = q , v q−1 ≤ 2v q−1 and ν q ≥ 2. The inequalities (S.182) and (S.183) reduces the inequality (S.179) to where we define which gives the Lieb-Robinson bound for e −iH1:qtm+1 , namely G q (x, t m+1 , X , Y). In order to apply Proposition 12 to this case, we first simplify the inequality (S.158) tõ where we use the condition (S.120) and the definition of ξ q in the second inequality. Thus, by using Proposition 12 with (S.184) and (S.187), we have the upper bound of where r * q,2 is chosen such that it satisfies From the condition (S.123), the condition (S.190) is satisfied by choosing r q,2 = 1. Hence, we have where we use ν q ≥ 2. Because of R q,1 R q,2 ≤ R q from the definition (S.169), we prove the inequality (S.167).

V. DETAILED PROOF OF THE COMPLEMENTARY STATEMENTS FOR THE MIDDLE-RANGE LIEB-ROBINSON BOUND
A. Proof of Proposition 10: Quasi-locality of operator after time-evolution

Statement
Proposition 10. Let˜ q be an length scale such that for fixed τ ≤ ∆t q . Then, under the assumption of (S.137), we obtain a decomposition the time-evolved Hamiltonian H q (H 1:q−1 , τ ) in the form of for an arbitrary positive integer s 0 , where h τ,Z[s˜ q ] is an operator which acts on the subset Z[s˜ q ].

Proof of Lemma 13
First, we obtain from the Lieb-Robinson bound of Eq. (S.198): On the summation with respect to Z, we obtain where we use the inequalities (S.197), (S.14) and (S.16) in the second, fourth and fifth inequalities, respectively. For an arbitrary subset X s.t. diam(X) ≤ s 0 ξ, we assume Then, the HamiltonianH satisfies GH (x, t, X )-Lieb-Robinson bound with We note that the same inequality holds for time-dependent Hamiltonian.

Proof of Proposition 11
We We notice that for L m * the summation with respect to Z s m * is not restricted to {Z s m * |Z s m * ∩ L = ∅}. We, in the following, aim to upper-bound the mth term L m in Eq. (S.224). By using the condition (S.215), we obtain the upper bound of Because the necessary condition for Z sm ∩ L = ∅ is given by the mth term L m is bounded from above by x z e −x+1 dx ≤ e −(µ−1)S0+1 z!, (S.232) where in the first inequality we use the condition µ > 1. By using the above inequality, for m ≤ d L,L /ξ, the inequality (S.231) reduces to s1≥1,s2≥1,··· ,sm≥1 where in the last inequality we use mD+m−1 m ≤ 2 mD+m−1 . Therefore, by combining the inequalities (S.230) and (S.233), we finally arrive at the inequality for L m as where we use the definition ofṽ in Eq. (S.217). In the same way, we can derive the similar upper bound for L m * . For L m * , only the difference is that the constraint of S m ξ ≥ d L,L in (S.230) is removed, and hence we obtain Then, by using the upper-bounds for {L m } m * m=1 , we have We thus prove the inequality (S.216). This completes the proof.
We note that the same inequality holds for time-dependent Hamiltonian.
Remark. When we consider a Hamiltonian with k-body interactions (k = O(1)) (i.e., k-local Hamiltonian), the inequality (S.240) is trivially obtained for α > D; for example, please see Sec. 2.3 in Ref. [111] or Appendix C in Ref.
[24]. However, the k-dependence of the Lieb-Robinson bound is roughly given by The Hamiltonian (S.238) have k = O(ξ D ), and hence the Lieb-Robinson velocity is not given by O(ξ), but given by O(ξ D+1 ). Therefore, in order to obtain the Lieb-Robinson velocity of O(ξ), we need to utilize the power-law decay as in (S.239) with α > 2D + 1.

Proof of Lemma 9
The proof is almost the same as Proposition 11. We obtain the same inequality as (S.223): and Thus, by using |Z| ≤ γr D for diam(Z) = r which comes from (S.10), we obtain By applying the inequalities (S.250) and (S.251), we obtain the upper bound of L m * by where we use m * ! ≥ (m * /e) m * and v 0 was defined in Eq. (S.241). We thus obtain the first inequality in (S.240) from the choice of m * as in Eq. (S.246). Finally, we need to prove For the purpose, we use the following upper bound for x ≥ ξv 0 |t|, which yields the inequality (S.253) and This completes the proof. for diam(X ), diam(Y) ≤ νξ (ν ≥ 1). Also, let H 0 be a Hamiltonian which satisfy the Lieb-Robinson bound as where the function G v ,t (x) is assumed to satisfy We then obtain the Lieb-Robinson function G (x, t, t , X , Y) for e −iH0t e −iH 0 t as where r * is defined as an integer which satisfies

Proof of Proposition 12
We consider the norm of the commutator where we determine r * afterward. From the Lieb-Robinson bound (S.163) and Lemma 2, the norm of O Y [sξ] is bounded from above by By using the decomposition (S.263), we obtain We then calculate an upper bound of ∞ s=r * +1 By using the inequality (S.258) with x 0 = r * , we obtain where in the last inequality we use (S.86). By applying the inequality (S.272) to (S.273), we finally obtain ≤2Cγe −(d X,Y −r * ξ−v|t|)/ξ (1 + 2r * /ν) D + 2e 2 γ 2 D!(2r * /ν + 1 + 4/ν) D F (r * ξ, t ) . (S.273) By choosing r * in the above inequality such that we obtain the inequality (S.259). Because of r * ≥ 1 and ν ≥ 1, we have and hence the condition (S.274) is simplified as This completes the proof.

VI. CONTRIBUTION OF THE LONG-RANGE INTERACTIONS
whereJ 0 and r 0 are constants which only depend on the parameters {D, g 0 , α, γ}.

B. Proof of Theorem 5
In the following, we defineα asα which givesη in Eq. (S.29) asη From the assumption of α > 2D + 1,α > 0. Also, the assumption 2 implies From Theorem 4, H ≤ t satisfies the G ≤ t (x, t, X )-Lieb-Robinson bound as By using Lemma 8, we first rewrite the function G ≤ t (x, t, X ) in the form of Because of v * ≥ v 0 1 from Sec. IV and the condition of |t| ≥ 1, we have where v 0 ≥ 1 is derived from the definitions (S.127) and (S.6) with g 0 ≥ 1. Thus, the condition (S.116) for 1 gives 1 − δ t ≥ 1/2, (S.285) which reduces (S.282) to For the derivation of (S.277), we first estimate the quasi-locality of H > t (H ≤ t , τ ) for τ < t. We prove the following proposition: Proposition 14. Let us consider a Hamiltonian A as follows: where A satisfies Assumptions 1 and 2 with g 0 = g a,0 and g = g a . We also consider a Hamiltonian A 0 which satisfying G 0 (x, t, X )-Lieb-Robinson bound as follows: We then obtain the decomposition of where J is a constant which depends on {γ, α, D, C}.