Tightening the Lieb-Robinson Bound in Locally-Interacting Systems

TL;DR

This paper introduces a new method that significantly tightens Lieb-Robinson bounds in locally-interacting quantum systems, especially at large local Hilbert space dimensions, improving understanding of information propagation limits.

Contribution

The authors develop a novel approach that yields much sharper Lieb-Robinson bounds, reducing the growth of the velocity with system parameters compared to previous bounds.

Findings

LR velocity grows slower with local Hilbert space dimension D

In Heisenberg model, v remains constant at large S

In multiorbital Hubbard models, v scales as √N instead of N

Abstract

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. We introduce a method that dramatically and qualitatively improves LR bounds in models with finite-range interactions. Most prominently, in systems with a large local Hilbert space dimension $D$, our method gives an LR velocity that grows much slower than previous bounds with $D$ as $D\to \infty$. For example, in the Heisenberg model with spin $S$, we find $v\leq$ const. compared to the previous $v\propto S$ which diverges at large $S$, and in multiorbital Hubbard models with $N$ orbitals, we find $v\propto \sqrt{N}$ instead of previous $v\propto N$, and similarly in the $N$-state truncated Bose-Hubbard model and Wen's quantum rotor model. Our bounds also scale qualitatively better in some systems when the spatial dimension or certain model parameters become large, for example in the $d$-dimensional quantum Ising model and perturbed toric code models. Even in spin-1/2 Ising and Fermi-Hubbard models, our method improves the LR velocity by an order of magnitude with typical model parameters, and significantly improves the LR bound at large distance and early time.

Tables

Table 1. Table 1: Comparison between previous LR speeds and those introduced in this paper, in the case of 2D spin-1/2 TFIM, 1D SU(2) FH model, and perturbed toric code (PTC) model  Kitaev ( 2003 ) with on site h ​ ∑ j σ ^ j x ℎ subscript 𝑗 subscript superscript ^ 𝜎 𝑥 𝑗 h\sum_{j}\hat{\sigma}^{x}_{j} perturbation. The LR speed v LR subscript 𝑣 LR v_{\mathrm{LR}} of this paper is upper bounded by the minimum of all the expressions in the different rows of each model. For each column divided into two subcolumns, the first subcolumn shows analytic expressions, and the second subcolumn shows their approximate numerical values. (For those bounds of the FH model depending on U / J 𝑈 𝐽 U/J we have chosen a representative point). The constant X y subscript 𝑋 𝑦 X_{y} is defined as the solution to the equation x ​ arcsinh ⁡ ( x ) = x 2 + 1 + y 𝑥 arcsinh 𝑥 superscript 𝑥 2 1 𝑦 x\operatorname{arcsinh}(x)=\sqrt{x^{2}+1}+y . Z U / J subscript 𝑍 𝑈 𝐽 Z_{U/J} is defined in Eq. ( 80 ) and plotted in Fig.  5 .

model	$v_{\text{LR}}$ (this paper)		$v_{\text{LR}}$ (Ref. Hastings (2010))
2D TFIM	$2X_0\sqrt{2Jh}$	$4.27\sqrt{Jh}$	$16eJ$	$43.5J$
	$8X_{\frac{1}{2}}J$	$15.1J$
	$8X_{0}h$	$12.1h$
1D FH	$2X_{\frac{3U}{4J}}J$	$4.14J_{\left(U=J\right)}$	$16eJ$	$43.5J$
	$8X_{0}J$	$12.1J$
	$Z_{U/J}J$	$7.05J_{\left(U=5J\right)}$
PTC	$8X_{\frac{1}{2}}h$	$15.1h$	$32e$	$87.0$
	$2X_0\sqrt{2h}$	$4.27\sqrt{h}$
	$8X_0$	$12.1$

Table 2. Table 2: Summary of qualitative improvements of v LR subscript 𝑣 LR v_{\mathrm{LR}} in this paper for models with large on-site Hilbert spaces or large spatial dimensions: S 𝑆 S indicates the size of spins in spin models, N 𝑁 N the maximal number of particles on a single-site in Hubbard models, and d 𝑑 d the spatial dimension. In Wen’s rotor model, we only list the result in the most interesting limit S → ∞ , J ≪ g formulae-sequence → 𝑆 much-less-than 𝐽 𝑔 S\to\infty,J\ll g , where the previous best methods from Refs.  Hastings ( 2010 ); Hamma  et al. ( 2009 ) are infinitely looser than our result.

model	$v_{\text{LR}}$ (this paper)	$v_{\text{LR}}$ (Ref.Hastings (2010))
large-$d$ TFIM	$\propto \sqrt{d}$	$\propto d$
large-$S$ (spin) TFIM	$\propto \sqrt{S}$	$\propto S^2$
large-$N$ SU($N$) FH	$\propto \sqrt{N}$	$\propto N$
Spin-$S$ Heisenberg	const.	$\propto S$
Wen’s rotor model,	$\propto \sqrt{dgJ}$	$\propto \sqrt{Jg}dS$ (Ref.Hamma et al. (2009))
$S\to \infty$ and $J\ll g$		$\propto dg$ (Ref.Hastings (2010))
$N$-state truncated BH	$\propto \sqrt{N}$	$\propto N$

Equations

$$\Vert [{\hat{A}}_X(t),{\hat{B}}_Y(0)]\Vert \le C{e}^{-s_{XY}/\xi }(e^{v|t|/\xi }-1),$$

$$\Vert [{\hat{\sigma }}_{\vec{r}}^x(t),{\hat{\sigma }}_{\vec{0}}^z(0)]\Vert \le {\int }_{-\pi }^{\pi }\cosh [\omega (\vec{k})t]e^{i\vec{k}\cdot \vec{r}}\frac{d^{d}k}{(2\pi {)}^d},$$

$$\Vert [{\hat{A}}_X(t),{\hat{B}}_Y(0)]\Vert \le C{\left(\frac{u|t|}{d_{XY}}\right)}^{d_{XY}}.$$

$$\hat{H}=\sum _j{h}_j{\hat{\gamma }}_j,$$

$$i{\dot{\hat{\gamma }}}_i(t)=[{\hat{\gamma }}_i(t),\hat{H}]=\sum _{j:⟨ij⟩\in G}{h}_j[{\hat{\gamma }}_i(t),{\hat{\gamma }}_j(t)],$$

$$i{\dot{\hat{\gamma }}}_i^B(t)=\sum _{j:⟨ij⟩\in G}{h}_j\{[{\hat{\gamma }}_i^B(t),{\hat{\gamma }}_j(t)]+[{\hat{\gamma }}_i(t),[{\hat{\gamma }}_j(t),\hat{B}(0)]]\}.$$

$$i{\dot{\hat{\tau }}}_i^B(t)={\hat{U}}^{†}(t)\sum _{j:⟨ij⟩\in G}{h}_j[{\hat{\gamma }}_i(t),[{\hat{\gamma }}_j(t),\hat{B}(0)]]\hat{U}(t).$$

$$\Vert {\hat{\gamma }}_i^B(t)\Vert -\Vert {\hat{\gamma }}_i^B(0)\Vert$$

$$\le {\int }_0^t\Vert {\dot{\hat{\tau }}}_i^B(t^{\prime })\Vert d{t}^{\prime }$$

$$\le \sum _{j:⟨ij⟩\in G}|h_j|{\int }_0^t\Vert [[{\hat{\gamma }}_j(t^{\prime }),\hat{B}(0)],{\hat{\gamma }}_i(t^{\prime })]\Vert d{t}^{\prime }$$

$$\le 2\sum _{j:⟨ij⟩\in G}|h_j|{\int }_0^t\Vert {\hat{\gamma }}_j^B(t^{\prime })\Vert d{t}^{\prime }.$$

$${\bar{\gamma }}_i^B(t)={\bar{\gamma }}_i^B(0)+2\sum _{j:⟨ij⟩\in G}|h_j|{\int }_0^t{\bar{\gamma }}_j^B(t^{\prime })d{t}^{\prime },$$

$${\dot{\bar{\gamma }}}_i^B(t)$$

$${\bar{\gamma }}_i^B(0)=\Vert [{\hat{\gamma }}_i,\hat{B}]\Vert .$$

$$\Vert {\hat{A}}^B(t)\Vert -\Vert {\hat{A}}^B(0)\Vert \le {\int }_0^t\sum _{i\in S(\hat{A})}2\Vert \hat{A}\Vert |h_i|{\bar{\gamma }}_i^B(t^{\prime })d{t}^{\prime }.$$

$${\dot{\bar{\Gamma }}}_i^B(t)$$

$${\dot{G}}_{ij}(t)=\sum _{k:⟨ik⟩\in G}{H}_{ik}{G}_{kj}(t),$$

$${\bar{\Gamma }}_i^B(t)=\sum _j{G}_{ij}(t){\bar{\Gamma }}_j^B(0),$$

$$\Vert [{\hat{\gamma }}_i(t),{\hat{B}}_Y(0)]\Vert \le \sum _{j\in Y}{G}_{ij}(t)\sqrt{\frac{h_j}{h_i}}2\Vert \hat{B}\Vert ,$$

$$\Vert [{\hat{A}}_X(t),{\hat{B}}_Y(0)]\Vert -\Vert [{\hat{A}}_X(0),{\hat{B}}_Y(0)]\Vert$$

$$\le 4\Vert \hat{A}\Vert \Vert \hat{B}\Vert {\int }_0^t\sum _{i\in X,j\in Y}\sqrt{|h_i{h}_j|}{G}_{ij}(t^{\prime })d{t}^{\prime }.$$

$${\hat{\gamma }}_i^2$$

$$\{{\hat{\gamma }}_i,{\hat{\gamma }}_j\}$$

$$\left[{\hat{\gamma }}_i,{\hat{\gamma }}_j\right]$$

$$\hat{H}=\sum _{i<j}{t}_{ij}i{\hat{c}}_i{\hat{c}}_j+\sum _{i<j<k<l}{U}_{ijkl}{\hat{c}}_i{\hat{c}}_j{\hat{c}}_k{\hat{c}}_l.$$

$${\hat{c}}_i^{†}$$

$$\{{\hat{c}}_i,{\hat{c}}_j\}$$

$$i{\dot{\hat{c}}}_i(t)=2\sum _j{t}_{ij}i{\hat{c}}_j(t)+2\sum _{j<k<l}{U}_{ijkl}{\hat{c}}_j(t){\hat{c}}_k(t){\hat{c}}_l(t).$$

$$i{\dot{\hat{c}}}_{im}(t)$$

$$\Vert {\hat{c}}_{im}(t)\Vert -\Vert {\hat{c}}_{im}(0)\Vert$$

$$\le 2{\int }_0^t\sum _{j<k<l}|U_{ijkl}|\Vert {\hat{c}}_{jm}(t^{\prime })+{\hat{c}}_{km}(t^{\prime })+{\hat{c}}_{lm}(t^{\prime })\Vert d{t}^{\prime }$$

$$+2{\int }_0^t\sum _j|t_{ij}|\Vert {\hat{c}}_{jm}(t^{\prime })\Vert d{t}^{\prime }.$$

$${\dot{c}}_{im}(t)$$

$$G_{ij}(t)\equiv {\int }_{-\pi }^{\pi }\frac{d^{d}k}{(2\pi {)}^d}{G}_{\alpha \beta }^{(\vec{k})}(t)e^{i\vec{k}\cdot ({\vec{r}}_I-{\vec{r}}_J)}.$$

$${\dot{G}}_{\alpha \beta }^{(\vec{k})}(t)=\sum _{\gamma =1}^l{H}_{\alpha \gamma }^{(\vec{k})}{G}_{\gamma \beta }^{(\vec{k})}(t),$$

Full text

I Introduction

The Lieb-Robinson (LR) bound Lieb and Robinson (1972); Hastings (2010) 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 Bravyi et al. (2006); Eisert and Osborne (2006); Hazzard et al. (2014); Gogolin and Eisert (2016), to transfer information in a quantum chanel Bose (2007); Epstein and Whaley (2017), and to prethermalize a system Calabrese and Cardy (2006); Rigol et al. (2007); Gong and Duan (2013); Eisert et al. (2015). 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 Hastings (2004a); Hastings and Koma (2006); Nachtergaele and Sims (2006); Cramer and Eisert (2006); Hastings (2006); Haegeman et al. (2013); Kliesch et al. (2014)), that systems with half-integer spin per unit cell are necessarily gapless (the Lieb-Schultz-Mattis-Hastings theorem Lieb et al. (1961); Hastings (2004b); Nachtergaele and Sims (2007)), the area law for entanglement entropy in ground states of gapped systems Hastings (2007a, b); Eisert et al. (2010), the stability of topological order Hastings and Wen (2005); Bravyi et al. (2010); Bravyi and Hastings (2011); Michalakis and Zwolak (2013), and the formal proof for the quantization of the Hall conductance as integer multiples of $e^{2}/h$ Hastings and Michalakis (2015). It is also the theoretical basis for analyzing the accuracy of some numerical algorithms Osborne (2006a, 2007a, 2007b); Woods et al. (2015); Woods and Plenio (2016); Tran et al. (2019). It has been generalized to systems with long-range power-law decaying interaction Hastings and Koma (2006); Eisert et al. (2013); Richerme et al. (2014); Gong et al. (2014); Foss-Feig et al. (2015); Maghrebi et al. (2016); Matsuta et al. (2017); Else et al. (2020); Guo et al. (2020); Tran et al. (2019); Chen and Lucas (2019a); Kuwahara and Saito (2020), to $n$-partite connected correlation functions Tran et al. (2017), and to dynamics with Markovian dissipation Hastings (2004c); Poulin (2010); Sweke et al. (2019) 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 Cheneau et al. (2012); Jurcevic et al. (2014); Richerme et al. (2014).

Despite the LR bounds’ fundamental consequences and efforts to tighten them Nachtergaele and Sims (2006); Hastings and Koma (2006); Hamma et al. (2009); Poulin (2010); Schuch et al. (2011); Roberts and Swingle (2016), the bounds are often extremely loose (except some special bounds on free particle systems Burrell and Osborne (2007); Cramer et al. (2008); Burrell et al. (2009); Damanik et al. (2014); Gong et al. (2019)), severely limiting their applicability. To see this, note that in locally interacting systems, previous LR bounds are typically of the form

[TABLE]

where $\hat{A}_{X},\hat{B}_{Y}$ are local operators supported on regions $X,Y$, respectively, $\hat{A}_{X}(t)=e^{i\hat{H}t}\hat{A}_{X}e^{-i\hat{H}t}$, $s_{XY}$ is the minimal graph-theoretical distance between points in $X$ and $Y$, and $v$ and $\xi$ are constants that only depend on the Hamiltonian $\hat{H}$ (and the lattice structure). There are two qualitative ways in which these bounds are loose. First, and arguably most importantly, the LR speeds $v$ appearing in existing bounds are typically much larger than the actual speed of information propagation. For example, in the one-dimensional (1D) transverse-field Ising model (TFIM) at critical point, the previous best LR speed is $8eJ\approx 21.7J$ while the actual speed determined from exact solution is $2J$. Previous bounds on $v$ can even become infinitely loose in certain limits: for example in the TFIM at large $J$, previous bounds grow $\propto J$; in the spin-$S$ Heisenberg model, previous bounds grow $\propto S$; in the classical limit of the SU$(N)$ Fermi-Hubbard (FH) model, previous bounds grow $\propto N$, while in all these cases the true LR speed remains finite. Previous bounds are also qualitatively loose at short time and large distance: as $t\to 0$, perturbative arguments show that the left-hand side (LHS) of Eq. (1) grows like $O(t^{cs_{XY}})$, where $c$ is a constant, much slower than the $O(t)$ in the RHS of Eq. (1); similarly, as $s_{XY}\to\infty$, the LHS actually decays faster than exponential. A simple example is non-interacting systems, where the LHS decays at large distance like $O(e^{-s_{XY}\ln s_{XY}})$.

In this paper, we focus on locally interacting systems (i.e. systems with finite-range interactions), and introduce a method that qualitatively tightens the LR bounds in all the aforementioned aspects. The key new ingredient in our method is the commutativity graph, which is a graphical tool that helps us taking advantage of the commutativity structure of the Hamiltonian terms, a feature that has been ignored by previous methods (except in a special model in Ref. Hamma et al. (2009)). By applying the Heisenberg equation and triangle inequality, we obtain a system of linear integral inequalities which only involve the unequal time commutator of the Hamiltonian terms. We then find a system of linear differential equations whose solution gives an upper bound for any solution to the integral inequalities. These differential equations can be naturally understood as the information flow equation on the commutativity graph.

If quantitative tightness is the sole criterion for applying the LR bound, one should solve these equations either analytically, or numerically if necessary, which can be done quickly and straightforwardly for systems with many thousands of sites [since the number of variables is proportional to the size of the system, see Eq. (10)]. 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 LR bound for $[\hat{\sigma}^{x}_{\vec{r}}(t),\hat{\sigma}^{z}_{\vec{0}}(0)]$ is

[TABLE]

where $\omega(k)$ corresponds to eigenvalue of the linear operator encoding the differential equations, which can be found by diagonalizing a $(d+1)\times(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. The LR speed obtained by this method vastly improves previous ones, as summarized in Tables 1 and 2.

As we will see, there are three broad scenarios where the LR speed appearing in our bounds is qualitatively tighter: (1) when the number of local degrees of freedom becomes large, e.g. the large-$S$ limit in TFIM, Heisenberg XYZ models, and Wen’s quantum rotor model; the Bose-Hubbard (BH) model truncated to allow at most $N$ particles per site, and the large-$N$ limit in SU($N$) FH model, (2) when the parameters of a set of commuting operators becomes large (e.g. large-$J$ limit in TFIM, perturbed toric code model), and (3) TFIM and Wen’s quantum rotor model in large spatial dimension. The bounds we present shed light onto these important limits which have previously resisted bounds with the qualitatively correct scaling. They also have important physical implications, for example in exactly solvable models with topological order (e.g., toric code Kitaev (2003) and string net models Levin and Wen (2005)), perturbed by an arbitrary bounded, local perturbation. Our results show that their LR velocity vanishes as the strength of the perturbation approaches zero, in contrast to prior bounds which retained a finite LR velocity in this limit. Even in situations with low spatial dimension and small local degrees of freedom, e.g. at $J=h$ in 2D TFIM or $U=J$ in 1D FH, our bound represents more than a 10-fold improvement.

Although expressions like Eq. (2) provide the tightest bounds in this paper, we can find simple analytic bounds while compromising the tightness only marginally. We will show that

[TABLE]

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 $\hat{A}_{X}(t),\hat{B}_{Y}(0)$ on the commutativity graph ($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 Tables 1, 2 (their precise relation is given in Sec. IV and Sec. V). This bound also tightens previous ones in the two aspects mentioned above, and the asymptotic behavior at small time and large distance is often tightest possible. 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 bound gives $\xi\propto(\ln\frac{J}{h})^{-1}$ in the limit $J/h\to\infty$, in agreement with the exact solution, while the previous best bound approaches a constant; in the spin-$S$ Heisenberg model, when $S\to\infty$, our bound gives $\xi\leq c_{1}$ (a constant), while the previous bound $\xi\leq c_{2}(S+1)$ diverges linearly in $S$. 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 Osborne (2006a, 2007a, 2007b); Woods and Plenio (2016); Tran et al. (2019). 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 on the commutativity graph. In Sec. III we specialize to translation 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 LR bounds in some paradigmatic models (TFIM, spin-$S$ Heisenberg, $N$-state truncated BH, SU($N$) FH, and Wen’s quantum rotor model). 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. In Sec. VIII we summarize our methods and results, and discuss their potential applications along with some future directions. In Appen. A we briefly introduce a more general method that unifies some of the seemingly different approaches in the main text.

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. Hastings and Koma (2006) 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. Hamma et al. (2009) 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.1. In Sec. II.2 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.

II.1 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 locally-interacting quantum system in arbitrary dimension with an arbitrary lattice structure, or on an arbitrary graph. The Hamiltonian $\hat{H}$ can in general be written as

[TABLE]

where $\hat{\gamma}_{j}$ are local Hermitian operators with unit norm $\|\hat{\gamma}_{j}\|=1$, and $h_{j}$ are constant parameters.

The commutativity graph $G$ of the Hamiltonian $\hat{H}$ is constructed as follows. Each local operator $\hat{\gamma}_{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 $\langle i,j\rangle$ if and only if the corresponding terms do not commute $[\hat{\gamma}_{i},\hat{\gamma}_{j}]\neq 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 $\hat{H}$. Nevertheless, for convenience we will simply speak of “the commutativity graph of $\hat{H}$” when a certain decomposition is implicitly assumed. We will discuss how to choose the decomposition at the end of Sec. II.2. In the following we will derive an LR bound based on a differential equation on the commmutivity graph of $\hat{H}$.

II.2 Upper bound for $\|[\hat{\gamma}_{i}(t),\hat{B}(0)]\|$

The Heisenberg equation for the operator $\hat{\gamma}_{i}(t)=e^{i\hat{H}t}\hat{\gamma}_{i}e^{-i\hat{H}t}$ is

[TABLE]

where the dot on $\hat{\gamma}_{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 $[\hat{A}(t),\hat{B}(0)]$ of local operators $\hat{A},\hat{B}$, i.e. operators with finite support. Let us first consider the special case when $\hat{A}=\hat{\gamma}_{i}$ is a term of the Hamiltonian. In this case, we are interested in the unequal time commutator $\hat{\gamma}^{B}_{i}(t)=[\hat{\gamma}_{i}(t),\hat{B}(0)]$, for an arbitrary local operator $\hat{B}$. Using Eq. (5) and the Jacobi identity, we have the evolution equation for $\hat{\gamma}^{B}_{i}(t)$

[TABLE]

Substituting $\hat{\gamma}^{B}_{i}(t)=\hat{U}(t)\hat{\tau}^{B}_{i}(t)\hat{U}^{\dagger}(t)$ where $\hat{U}(t)$ is the unitary operator satisfying $\hat{U}(0)=1$ and $i\dot{\hat{U}}(t)=-\sum_{j:\langle ij\rangle\in G}h_{j}\hat{\gamma}_{j}(t)\hat{U}(t)$ into Eq. (6), we have

[TABLE]

Now we bound the time dependence of the operator norm of ${\hat{\gamma}}_{i}^{B}(t)$, the fundamental object controlled by LR bounds. Since ${\hat{\gamma}}_{i}^{B}$ and ${\hat{\tau}}_{i}^{B}$ are related by unitary transforms, $\|\hat{\gamma}^{B}_{i}(t)\|=\|\hat{\tau}^{B}_{i}(t)\|$ (recall that $\|\hat{U}\hat{A}\hat{V}^{\dagger}\|=\|\hat{A}\|$ for unitary operators $\hat{U},\hat{V}$). Using this and applying the basic inequalities $\|\hat{A}+\hat{B}\|\leq\|\hat{A}\|+\|\hat{B}\|$ and $\|\hat{A}\cdot\hat{B}\|\leq\|\hat{A}\|\cdot\|\hat{B}\|$, we obtain

[TABLE]

where in the last line we used $\|\hat{\gamma}_{i}(t^{\prime})\|=\|\hat{\gamma}_{i}\|=1$.

Although Eq. (II.2) 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 Chandra and Davis (1976) to show that $\|\hat{\gamma}^{B}_{i}(t)\|\leq\bar{\gamma}^{B}_{i}(t)$, where ${\bar{\gamma}}_{i}^{B}(t)$ is the solution to Eq. (II.2) with the inequality replaced by equality:

[TABLE]

where $\bar{\gamma}^{B}_{i}(0)=\|\hat{\gamma}^{B}_{i}(0)\|$. Differentiating Eq. (9), one obtains the first order linear differential equation

[TABLE]

with initial condition

[TABLE]

We define the support of operator $\hat{B}$ on the commutativity graph as $S(\hat{B})=\{i|~{}[\hat{\gamma}_{i},\hat{B}]\neq 0\}$. Then the initial condition can be simplified as $\bar{\gamma}^{B}_{i}(0)=2\|\hat{B}\|(i\in S(\hat{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 $\hat{A}^{B}(t)=[\hat{A}(t),\hat{B}(0)]$, where $\hat{A}$ is an operator not in the Hamiltonian, we can perform similar calculations in Eqs. (5-II.2) to get

[TABLE]

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 form

[TABLE]

where $\bar{\Gamma}^{B}_{i}(t)=\sqrt{h_{i}}\bar{\gamma}^{B}_{i}(t)$. We denote the real symmetric coefficient matrix as $H_{ij}=2\sqrt{|h_{i}||h_{j}|}(\langle ij\rangle\in G)$. The Green’s function $G_{ij}(t)$ is defined as the solution to the differential equation

[TABLE]

with initial condition $G_{ij}(0)=\delta_{ij}$, so any initial value problem can be obtained by linear combinations of the Green’s functions:

[TABLE]

and therefore the LR bound of $\|[\hat{\gamma}_{i}(t),\hat{B}_{Y}(0)]\|$ can be expressed as

[TABLE]

where $Y=S(\hat{B})$ is the support of operator $\hat{B}$ on the commutativity graph (from now on we use $\hat{A}_{X}$ to mean that operator $\hat{A}$ has support $X$ on the commutativity graph). The LR bound for the unequal time commutator between arbitrary local operators $\|[\hat{A}_{X}(t),\hat{B}_{Y}(0)]\|$ can be obtained by inserting Eq. (16) into Eq. (12):

[TABLE]

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.1. 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 $\hat{H}$ 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) Nussinov and Ortiz (2009); Cobanera et al. (2011) if we require that $\hat{\gamma}^{2}_{i}=1$ and $\hat{\gamma}_{i},\hat{\gamma}_{j}$ for $i\neq j$ either commute or anti-commute. More precisely, the algebraic relations of $\{\hat{\gamma}_{i}\}$ are related to their commutativity graph $G$ as follows:

[TABLE]

For example, in spin-$1/2$ systems $\hat{\gamma}_{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. (II.2) still fail to specify the decomposition uniquely, we would choose [among those satisfying Eq. (II.2)] a decomposition such that the total number of $\hat{\gamma}_{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]\|\leq 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.

II.3 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 locally-interacting 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 as

[TABLE]

where $t_{ij}$ and $U_{ijkl}$ are totally antisymmetric, and the Majorana operators satisfy

[TABLE]

We first upper bound the norm of the unequal time anticommutator $\hat{c}_{ij}(t)=\{\hat{c}_{i}(t),\hat{c}_{j}(0)\}$. The Heisenberg equation for the operator $\hat{c}_{i}(t)=e^{i\hat{H}t}\hat{c}_{i}e^{-i\hat{H}t}$ is

[TABLE]

Using Eq. (21), we have the evolution equation for $\hat{c}_{im}(t)$

[TABLE]

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

[TABLE]

Eq. (II.3) is the analog of Eq. (II.2). Therefore, using the same argument as we used in Sec. II.2, we can conclude that $\|\hat{c}_{im}(t)\|\leq c_{im}(t)$ where $c_{im}(t)$ is the solution to the differential equation

[TABLE]

with initial condition $c_{im}(0)=\|\hat{c}_{im}(0)\|=2\delta_{im}$. Once we know the LR bound for $\{\hat{c}_{i}(t),\hat{c}_{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.2 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.4 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 $\alpha=1,\ldots,l$. We use capital $I$ to label unit cells with coordinate labeled by $\vec{r}_{I}$ [for example, if $\{\boldsymbol{\alpha}_{1},\ldots,\boldsymbol{\alpha}_{d}\}$ is a set of primitive lattice translation vectors, then coordinate $\vec{r}_{I}=(n_{1},n_{2},\ldots,n_{d}),~{}n_{a}\in\mathbb{Z},~{}a=1,\ldots,d$ labels the primitive cell at physical position $\mathbf{R}_{I}=\sum^{d}_{a=1}n_{a}\boldsymbol{\alpha}_{a}$], so that every vertex $i$ can be labeled by a pair $(I,\alpha)$. Taking advantage of translation invariance, we can use the Fourier integral representation of the Green’s function $G_{ij}(t)\equiv G_{I\alpha;J\beta}(t)$,

[TABLE]

[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 $(-\pi,\pi]^{d}$ rather than the first Brillouin zone in typical solid state contexts.] Inserting Eq. (25) into Eq. (14) we get

[TABLE]

where the $l\times l$ coefficient matrix $H^{(\vec{k})}_{\alpha\beta}$ is the Fourier transform of $H_{I\alpha;J\beta}$ (notice that due to translation invariance $H_{I\alpha;J\beta}$ only depends on $\vec{r}_{I}-\vec{r}_{J}$):

[TABLE]

The formal solution to Eq. (26) is $G^{(\vec{k})}_{\alpha\beta}(t)=[e^{H^{(\vec{k})}t}]_{\alpha\beta}$, where $H^{(\vec{k})}$ is considered as an $l\times l$ matrix.

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 $\vec{k}$ except at infinity, and since the integrand is periodic in shifts of the real axis by $2\pi$, we can send $\vec{k}\to\vec{k}+i\vec{\kappa}$ without changing the result of the integral, as illustrated in Fig. 2. Therefore we have

[TABLE]

where $\vec{r}=\vec{r}_{I}-\vec{r}_{J}$, in the second line we used the triangle inequality for the integral, and in the third line we used the inequality $|[e^{H^{(\vec{k}+i\vec{\kappa})}t}]_{\alpha\beta}|\leq|[e^{H^{(i\vec{\kappa})}t}]_{\alpha\beta}|$ [which follows from the triangle inequality since the coefficients $H_{I\alpha;J\beta}$ in Eq. (27) are non-negative]. In the last line we used the inequality $|[e^{H^{(i\vec{\kappa})}t}]_{\alpha\beta}|\leq\|e^{H^{(i\vec{\kappa})}t}\|\leq c_{\vec{\kappa}}e^{\omega_{m}(i\vec{\kappa})t}$, where $\omega_{m}(i\vec{\kappa})$ is the eigenvalue of $H^{(i\vec{\kappa})}$ with largest magnitude [$\omega_{m}(i\vec{\kappa})$ must be real positive according to the Perron-Frobenius theorem], $c_{\vec{\kappa}}$ is a function of $\vec{\kappa}$ independent of $t$, and the last inequality can be proved by diagonalizing $H^{(i\vec{\kappa})}$ into Jordan canonical form. Therefore we arrive at the important result:

[TABLE]

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 $\vec{\kappa}=\vec{\kappa}(\vec{r},t)$ to depend on $\vec{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 $\vec{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 $\vec{\kappa}=\vec{\kappa}_{0}\equiv\mathrm{sgn}(\vec{r})\kappa_{0}$ [where $\mathrm{sgn}(\vec{r})\equiv(\mathrm{sgn}(r_{1}),\mathrm{sgn}(r_{2}),\ldots,\mathrm{sgn}(r_{d}))$], we get the old form of LR bound with exponential decaying tail

[TABLE]

where $|\vec{r}|=\sum^{d}_{j=1}r_{j}$ is the Manhattan distance between the corresponding unit cells. Therefore the LR speed of $|G(\vec{r},t)|$ must be upper bounded by an optimal choice of $\kappa_{0}$:

[TABLE]

Notice that while $\vec{\kappa}_{0}$ depends on $\mathrm{sgn}(\vec{r})$, in typical models $\omega_{m}(i\vec{\kappa}_{0})$ does not depend on $\mathrm{sgn}(\vec{r})$ [see Sec. V]. In the case that $\omega_{m}(i\vec{\kappa}_{0})$ does depend on $\mathrm{sgn}(\vec{r})$, Eq. (31) would represent a directional LR speed, and to get the overall LR speed we should take the maximum over all possible $2^{d}$ values of $\mathrm{sgn}(\vec{r})$.

It is immediately clear from Eq. (17) that the LR speed of the Green’s function gives an upper bound for the LR speed of the unequal time commutator between any local operators $\|[\hat{A}_{X}(t),\hat{B}_{Y}(0)]\|$, since, by inserting Eq. (30) into Eq. (17), we get

[TABLE]

where $C=4\|\hat{A}\|\|\hat{B}\|\sqrt{h_{X}h_{Y}}c_{\vec{\kappa}_{0}}/\omega_{m}(i\vec{\kappa}_{0})$, $h_{X}\equiv\max_{i\in X}h_{i}$, $r_{XY}=\min_{i\in X,j\in Y}r_{ij}$ and we assumed $X\cap Y=\emptyset$. This proves that $\|[\hat{A}_{X}(t),\hat{B}_{Y}(0)]\|$ has the same LR speed $v_{\text{LR}}$.

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_{\text{LR}}$, for example the methods given in Ref. Hastings (2010) or Ref. Bravyi et al. (2006), 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

[TABLE]

where $G^{(n)}_{ij}=[H^{n}]_{ij}$. 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}=\prod_{e\in p}w^{d_{e}}_{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

[TABLE]

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

[TABLE]

for some constant $c$. Inserting Eq. (34), the above inequality is equivalent to

[TABLE]

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$:

[TABLE]

which can be simplified to

[TABLE]

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 graph-theoretical distance in the commutativity graph, so $u$ is in general different from $v_{\text{LR}}$. But since $c_{1}\leq d_{ij}/r_{ij}\leq 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 $\hat{H}$), existence of a finite $v_{\text{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 $\lambda$ edges and the weights satisfy $H_{ij}\leq h$ for all pairs of neighboring vertices $i,j$. With these assumptions we have $G^{(n)}_{ij}\leq h^{n}\lambda^{n},~{}\forall i,j\in G,~{}n\geq 0$. Inserting into Eq. (34), we have

[TABLE]

for $t<d_{ij}/\lambda he$. This proves Eq. (38) if we choose, for example, $u=2\lambda he$.

Now we use the bound Eq. (35) on $G_{ij}(t)$ to bound a general unequal-time commutator $\|[\hat{A}_{X}(t),\hat{B}_{Y}(0)]\|$. Inserting Eq. (35) into Eq. (17), we have

[TABLE]

for $t\leq d_{XY}/u$, where $h_{XY}$ is a geometric factor defined as

[TABLE]

In the third line of Eq. (IV) we used the inequality$(1+1/d_{ij})^{d_{ij}}\leq e$, and in the last line we used $(ut/x)^{x}\leq(ut/y)^{y}e^{y-x}$ for $t\leq y/u\leq x/u$  [this inequality holds since the function $(ut/x)^{x}e^{x}$ is decreasing in $x$ for $x\geq ut$]. For $t>d_{XY}/u$ we have the trivial bound $\|[\hat{A}_{X}(t),\hat{B}_{Y}(0)]\|\leq 2\|\hat{A}\|\|\hat{B}\|\leq 2e\|\hat{A}\|\|\hat{B}\|\left(\frac{ut}{d_{XY}+1}\right)^{d_{XY}+1}$. In summary we have

[TABLE]

where $C=2e\max\{1,2c_{\vec{\kappa}_{0}}\frac{h_{XY}}{u}\}\|A\|\|B\|$ for $X\cap Y=\emptyset$ and $\forall t\geq 0$. This completes the proof of the LR bound in Eq. (3). [Eq. (42) is actually a slightly tighter version].

The 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 $\eta(\vec{r})$ of a bound $f(\vec{r},t)$ is defined by its limiting behavior $f(\vec{r},t)\propto t^{\eta(\vec{r})}$ as $t\to 0$. We now argue that if a minimal Clifford decomposition [as defined below Eq. (II.2)] of the Hamiltonian is used in Eq. (4), then the bound in Eq. (42) has a small-$t$ exponent that generically agree with the exact result, and thus are the tightest possible. For simplicity, consider the LR bound for $[\hat{\gamma}_{i}(t),\hat{\gamma}_{j}(0)]$. The Baker-Campbell-Hausdorff formula gives

[TABLE]

where $\mathrm{ad}_{\hat{H}}$ is the adjoint map of $\hat{H}$ acting on the space of operators defined as $\mathrm{ad}_{\hat{H}}(\hat{B})\equiv[\hat{H},\hat{B}]$, $P_{ij}$ is the set of shortest paths connecting $i$ and $j$ (note: by our convention a path $p\in P_{ij}$ does not contain $i$ and $j$), the string operator $\hat{S}_{p}$ along a path $p\in P_{ij}$ is defined as $\hat{S}_{p}=\prod_{k\in p}h_{k}\hat{\gamma}_{k}$ (with a suitable ordering). Therefore the small-$t$ exponent is $\eta_{ij}=d_{ij}-1$ provided that the sum on the second line of Eq. (IV) is nonzero. A sufficient condition for this to be true is if no two string operators $\hat{S}_{p}$ and $\hat{S}_{p^{\prime}}$ along different paths $p\neq p^{\prime}$ are proportional to each other [Notice that each individual term $\hat{\gamma}_{i}\hat{S}_{p}\hat{\gamma}_{j}$ is a product of invertible operators (due to the relation $\hat{\gamma}^{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(\hat{\gamma}_{i}),Y=S(\hat{\gamma}_{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 $\eta(\vec{r})=1$, as in Eq. (1) for example 111Refs. Bravyi et al. (2006); Osborne (2006b); Chen and Lucas (2019b); Chessa and Giovannetti (2019) also derived LR bounds which decay superexponentially in distance, and have small-$t$ exponents proportional to the Manhattan distance in real space. But these bounds are still looser than ours in most cases of TFIM. For example, their small-$t$ exponent is $|\vec{r}|$ while our bound in Eq. (51) has small-$t$ exponent $2|\vec{r}|$, and at large $r\equiv|\vec{r}|$, their bound decays like $e^{-r\ln r}$ while our bound decays like $e^{-2r\ln r}$, due to the fact that the distance on the commutativity graph of TFIM is twice the distance on real space. Also, Ref. Osborne (2006b) only gives the 1D result and the method is hard to generalize to higher dimensions, while the higher dimensional cases of Refs. Bravyi et al. (2006); Chen and Lucas (2019b) requires solving a hard combinatorical problem of summing over paths, the computational difficulty of which may grow exponentially with system size..

V Examples

In this section, we apply the techniques developed above to calculate the LR velocity for the bound of Eq. (3) explicitly for five models as examples: the $d$-dimensional TFIM (Sec. V.1), the spin-$S$ Heisenberg model (Sec. V.2), the $N$-state truncated BH model (Sec. V.3), the SU($N$) FH model (Sec. V.4), and Wen’s quantum rotor model (Sec. V.5). The techniques are general and can be straightforwardly applied to arbitrary models. All these examples demonstrate significant quantitative improvements over the previous LR speeds, and in a number of limits (e.g. large-$d$, large-$S$ or large-$N$) our bounds are qualitatively tighter for having better scalings. We will introduce techniques to further tighten the LR velocity in Sec. VI.

V.1 The $d$-dimensional TFIM

In the following we will first derive a bound in the familiar spin-1/2 case, and then extend this result to arbitrary spin in Sec. V.1.2.

V.1.1 Spin-1/2 case

The Hamiltonian for the $d$-dimensional hypercubic lattice Ising model with a transverse field is

[TABLE]

with $\hat{\gamma}_{\vec{r},j}=\hat{\sigma}^{z}_{\vec{r}}\hat{\sigma}^{z}_{\vec{r}+\hat{e}_{j}}$ where $\hat{e}_{j}$ is the unit vector in the $j$th direction, and $\hat{\gamma}_{\vec{r},0}=\hat{\sigma}^{x}_{\vec{r}}$. We will assume $J\geq 0$, $h\geq 0$ (the resulting LR bound only depends on $|J|$ and $|h|$). We have written the Hamiltonian in this form so that Eq. (II.2) is satisfied. It is easily verified that the square of each term is equal to unity $\hat{\gamma}^{2}_{\vec{r},\alpha}=1$, $\alpha=0,1,2,\ldots,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)

[TABLE]

where $H^{(\vec{k})}$, defined in Eq. (27), is in this case a $(d+1)\times(d+1)$ matrix

[TABLE]

so that

[TABLE]

where $\mathds{1}$ is the identity matrix and

[TABLE]

with $\cos\vec{k}\equiv\sum^{d}_{j=1}\cos k_{j}$. The eigenvalues of $H^{(\vec{k})}$ are [math] and $\omega^{\pm}(\vec{k})=\pm\omega(\vec{k})=\pm\sqrt{8Jh(d+\cos\vec{k})}$, from which we can compute the upper bound for the LR speed given by Eq. (31):

[TABLE]

where the constant $X_{y}$ is defined as the solution to the equation $x\operatorname{arcsinh}(x)=\sqrt{x^{2}+1}+y$, and $X_{0}\approx 1.50888$. The previous best bound for $v_{\mathrm{Ising}}$ is $8edJ\approx 21.7dJ$ obtained from the method in Ref. Hastings (2010), so our bound in Eq. (49) is a major, parametric improvement when $J$ or $d$ is large. Here we emphasize that to our knowledge this is the first LR bound whose speed scales sublinearly with spatial dimension $v\propto\sqrt{d}$–all previously known LR bounds have velocities that scale linearly with $d$ 222The reason for this is that most of the previous LR bounds are very generic to include a large family of models. As we will see at the end of Sec. V.4, the speed of $d$-dimensional free fermions grows linearly at large $d$, so $v\leq(\text{const.})\times d$ is the best thing a generic bound can do.. Even for $J=h,d=2$, Eq. (49) represents a 10-fold improvement over the previous bound. See Table 1 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.2.

The LR bound for arbitrary local operators can be obtained from Eqs. (16) and (12). For example, for $\|[\hat{\sigma}^{x}_{\vec{r}}(t),\hat{\sigma}^{z}_{\vec{0}}(0)]\|$ we have

[TABLE]

Finally, $\|[\hat{\sigma}^{x}_{\vec{r}}(t),\hat{\sigma}^{z}_{\vec{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=\vec{r}$ and $Y=\vec{0}$ on the commutativity graph, which is exactly $2|\vec{r}|$. Therefore, the parameter $u$ needed to satisfy Eq. (38) can be chosen as $u=2v_{\mathrm{Ising}}$. Inserting Eq. (35) into Eq. (50) we get the simple LR bound

[TABLE]

where $C$ is a constant.

V.1.2 Extending to arbitrary spin

The $d$-dimensional spin-$S$ TFIM is defined by the Hamiltonian

[TABLE]

where $\hat{S}^{\alpha}$ is the spin-$S$ operator satisfying $[\hat{S}^{\alpha},\hat{S}^{\beta}]=i\sum_{\gamma}\epsilon^{\alpha\beta\gamma}\hat{S}^{\gamma}$, $\alpha,\beta,\gamma\in\{x,y,z\}$. Eq. (52) is a direct generalization of the spin-1/2 case Eq. (44), so in this case we can simply set $\hat{\gamma}_{\vec{r},j}=\frac{1}{S^{2}}\hat{S}^{z}_{\vec{r}}\hat{S}^{z}_{\vec{r}+\hat{e}_{j}}$ and $\hat{\gamma}_{\vec{r},0}=\frac{1}{S}\hat{S}^{x}_{\vec{r}}$, so that the commutativity graph for $\hat{\gamma}_{\vec{r},j},\hat{\gamma}_{\vec{r},0}$ is the same as in Fig. 3, with $J,h$ rescaled to $4JS^{2},2hS$, respectively. In this way we can borrow the result from the spin-1/2 case Eq. (49) to get $v_{\mathrm{LR}}\leq 2X_{0}\sqrt{dJh}(2S)^{3/2}$. This is already qualitatively tighter than the previous best bound $v_{\mathrm{LR}}\leq 8edJ(2S)^{2}$.

But we can even do much better than this. Recall from the discussion at the end of Sec. II.2 that the resulting LR bound would typically be tightest if we decompose $\hat{H}$ into a form where any two terms either commute or anti-commute. Now the difficulty is that in the spin-$S$ case the spin operators $\hat{S}^{x},\hat{S}^{z}$ neither commute nor anti-commute. To overcome this, we can decompose each spin-$S$ operator as a sum of spin-1/2 Pauli operators:

[TABLE]

by standard addition of angular momentum. The physical Hilbert space on each site is spanned by all the eigenstates of the operator $(\hat{S}^{x})^{2}+(\hat{S}^{y})^{2}+(\hat{S}^{z})^{2}$ with eigenvalue $S(S+1)$, which are simply the states that are fully symmetric under permutation of indices $a$. In the following we will derive an LR bound for operators in the enlarged Hilbert space, which automatically gives a bound on operators acting on the physical Hilbert space, since the physical operators (e.g. $\hat{H},\hat{S}^{\alpha}_{j}$) do not couple physical states to unphysical states.

Inserting Eq. (53) into Eq. (52), the Hamiltonian in the enlarged space reads

[TABLE]

where $\hat{Z}^{ab}_{ij}\equiv\hat{\sigma}^{z}_{i,a}\hat{\sigma}^{z}_{j,b}$, $\hat{X}^{a}_{j}\equiv\hat{\sigma}^{x}_{j,a}$. The commutativity graph of the enlarged Hamiltonian in Eq. (54) can be regarded as a decorated version of the spin-1/2 case shown in Fig. 3, where each triangle is replaced by a set of $2S$ triangles and each circle is replace by a set of $\binom{2S}{2}$ circles, each of them linked to a distinct pair of two neighboring triangles, as shown in Fig. 4. In this case Eq. (10) becomes

[TABLE]

Let $X_{i}\equiv\sum_{1\leq a\leq 2S}\bar{X}^{a}_{i}$ and $Z_{ij}\equiv\sum_{1\leq a,b\leq 2S}\bar{Z}^{ab}_{ij}$, we then have

[TABLE]

This is exactly the same as the spin-1/2 case with substitution $h\to 2hS$. Therefore, borrowing the result in Eq. (49), we have

[TABLE]

growing as $\sqrt{S}$ at large $S$, a remarkable improvement over the previous bound which grows quadratically in $S$. Notice also that if we take the limit $S\to\infty$ at which $JS$ stays constant (the classical limit), then our LR speed stays finite while the previous bound diverges linearly in $S$.

V.2 The spin-$S$ Heisenberg XYZ model

The Hamiltonian for the spin-$S$ Heisenberg XYZ model is

[TABLE]

where the coefficient is normalized so that the excitation spectrum has a well-defined large-$S$ limit Haldane (1983). To get the tightest possible LR bound, we use the same trick as in Sec. V.1.2, decomposing each spin-$S$ operator as a sum of spin-1/2 Pauli operators as done in Eq. (53). The enlarged Hamiltonian reads

[TABLE]

where $\hat{X}^{ab}_{ij}\equiv\hat{\sigma}^{x}_{ia}\hat{\sigma}^{x}_{jb}$, and similarly for $\hat{Y}^{ab}_{ij},\hat{Z}^{ab}_{ij}$. The LR bounds for the norm of the commutators $[\hat{X}^{ab}_{ij}(t),\hat{B}],[\hat{Y}^{ab}_{ij}(t),\hat{B}]$, and $[\hat{Z}^{ab}_{ij}(t),\hat{B}]$ for an arbitrary local operator $\hat{B}$ are obtained from the solution to the differential equation (10), which in the current case becomes (we omit the superscript $B$ for notational simplicity)

[TABLE]

where it is assumed that $i,j$ are neighboring sites, and there are two other equations obtained by permuting $X,Y,Z$. The initial condition is similar to Eq. (11). Summing over indices $a,b$, we get

[TABLE]

where we use the notation

[TABLE]

The largest eigenfrequency in a $d$-dimensional hypercubic lattice is $\omega_{m}(\vec{k})=4J_{m}[\cos\vec{k}+d-\frac{1}{2S}]$, where $J_{m}=\|\mathbf{J}\|$, so

[TABLE]

We now make a comparison with previous best bound and the exact solution in the $S\to\infty$ limit. For simplicity we focus on the isotropic case $J_{x}=J_{y}=J_{z}=J$. We have $J_{m}=2J$, so our LR speed remains finite $v_{\mathrm{LR}}\leq 8dJX_{1}$ in the limit $S\to\infty$, while the previous bound $v_{\mathrm{LR}}\leq 8de\|\hat{h}_{ij}\|=16deJ(S+1)$ diverges linearly. Our bound is qualitatively tight in the sense that in the large-$S$ limit, it is only a finite factor of $2X_{1}\approx 4.47$ bigger than the exact speed $4dJ$ (the spin-wave group velocity), as calculated from the mean-field dispersion relation in Ref. Haldane (1983). We emphasize that this is one of the few known examples of a finite LR bound for a system with an infinite local Hilbert space Cramer et al. (2008); Nachtergaele et al. (2009); Schuch et al. (2011).

V.3 Truncated BH model

We now consider a spin-$S$ XY model with an additional $(S+S_{j}^{z})^{2}$ interaction at each site:

[TABLE]

where $\hat{S}^{\pm}=\hat{S}^{x}\pm i\hat{S}^{y}$. At each site $j$ there are $2S+1$ states $\{|m-S\rangle\}^{2S}_{m=0}$, and the action of operators $\hat{S}^{\pm}$ is given by

[TABLE]

therefore, in the limit $S\to\infty$ the Hamiltonian in Eq. (64) reproduces the BH Hamiltonian

[TABLE]

if we identify the mapping

[TABLE]

In the following we derive an LR bound for the truncated system in Eq. (64). Using the same trick as in Sec. V.2, the analog of Eq. (V.2) is (notice that the term linear in $\hat{S}^{z}$ commutes with $\hat{H}$ and can therefore be dropped)

[TABLE]

and another equation similar to the first one with $X$ and $Y$ exchanged. Summing over $a,b$ and doing a Fourier transform, we get

[TABLE]

where $X_{\vec{k}}$ is the Fourier transform of $X_{i}\equiv\sum_{j:\langle ji\rangle,a,b}X^{ab}_{ij}$, and $Z_{\vec{k}}$ is the Fourier transform of $Z_{i}\equiv\sum_{a\neq b}Z^{ab}_{i}$. At large $S$, we have $\omega_{m}(i\vec{\kappa})\approx 4\sqrt{2SUJd}\cosh\frac{\kappa}{2}$, so

[TABLE]

qualitatively tighter than the bound from previous methods which diverges linearly in $S$. It is known that the actual speed of information propagation in the $S\to\infty$ limit also diverges Eisert and Gross (2009), and mean-field approximation predicts that the BH model has $v\propto\sqrt{N}$ for an initial coherent state with average occupation number $N$ Menotti and Trivedi (2008), which suggests that the asymptotic behavior of our bound $v\propto\sqrt{S}$ is tight.

V.4 The SU($N$) FH model

We next take the SU($N$) FH model as an example of interacting fermions, using the method in Sec. II.3. The Hamiltonian in a general lattice is

[TABLE]

For $N=2$, this is the usual FH model. In the following, we limit our discussion to a bipartite lattice. We first split the fermion creation and annihilation operators to Majorana operators

[TABLE]

where $E$ is the set of even sites, $\bar{\sigma}=-\sigma$, and the Majorana operators satisfy Eq. (II.3). The Hamiltonian in the Majorana representation is

[TABLE]

The norm of the operator $\hat{c}_{j\sigma;l\sigma^{\prime}}(t)\equiv\{\hat{c}_{j\sigma}(t),\hat{c}_{l\sigma^{\prime}}(0)\}$ is upper bounded by the solution to Eq. (24), which in the current case become (for notational simplicity we omit the labels $l,\sigma^{\prime}$)

[TABLE]

with initial condition $c_{j\sigma}(0)=2\delta_{jl}\delta_{\sigma\sigma^{\prime}}$. Eq. (74) can be easily solved by a Fourier transform. In a $d$-dimensional hypercubic lattice, we obtain the LR bound

[TABLE]

where $I_{\alpha}(t)$ is the modified Bessel function of the first kind, and $I_{\vec{\alpha}}(t)\equiv\prod^{d}_{m=1}I_{\alpha_{m}}(t)$. The LR speed $v_{\mathrm{FH}}$ can be calculated from Eq. (31) with $\omega_{m}(\vec{k})=2J\cos\vec{k}+(N-1/2)U$, which gives

[TABLE]

At small $U/J$, this is $3.02dJ$, significantly improving the previous best bound $8edJN\approx 21.7dJN$ Hastings (2010).

At large $N$, the above LR speed scales like $N/\ln N$, which is already qualitatively tighter than the previous bound. But we can do even better than this, to prove an LR speed that grows like $\sqrt{N}$. This can either be achieved by using the commutativity graph method, or by using the fermion method with a slightly different treatment. In the following we take the latter approach for simplicity.

We go back to Eq. (73) and write down the Heisenberg equation for the operators $\hat{c}_{j\sigma}$ and $\hat{u}_{j\sigma\sigma^{\prime}}\equiv\hat{c}_{j\sigma}\hat{c}_{j\bar{\sigma}}\hat{c}_{j\sigma^{\prime}}\hat{c}_{j\bar{\sigma}^{\prime}}$

[TABLE]

where $p_{j}=(-1)^{(j\in E)}$. Following the derivations in Eqs. (6-11), one can show that the norm of the operators $\{\hat{c}_{j\sigma}(t),\hat{c}_{l\tau}(0)\},[\hat{u}_{j\sigma\sigma^{\prime}},\hat{c}_{l\tau}(0)]$ are upper bounded by the solution to the differential equations

[TABLE]

with initial condition $c_{j\sigma}(0)=2\delta_{jl}\delta_{\sigma\tau},u_{j\sigma\sigma^{\prime}}(0)=2\delta_{jl}(\tau\in\{\sigma,\bar{\sigma},\sigma^{\prime},\bar{\sigma}^{\prime}\})$. The largest eigenfrequency is

[TABLE]

and the LR speed is again computed from Eq. (31). The result is $v_{\mathrm{FH}}\leq Z_{(N-1)dU/J}J$, where

[TABLE]

At large $y$, $Z_{y}\approx X_{0}\sqrt{2y}$, so

[TABLE]

a significant qualitative improvement over the previous linear growth. In the classical limit $N\to\infty$ where $\langle\hat{n}_{j}\rangle/N$ and $UN$ stays constant, our LR speed remains finite while previous bound diverges linearly in $N$.

For the 1D SU(2) case, 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 $U/J\leq 80$. Results that also improve the previous bound at large $U$ are derived in Sec. VI.3.

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}=|\vec{x}_{j}-\vec{x}_{l}|$, and $u=\min\{2X_{\frac{(2N-1)U}{4dJ}}dJ,Z_{(N-1)dU/J}J\}$, and we have

[TABLE]

for $\forall j,l,\sigma,\sigma^{\prime},t$, where $C$ is a constant.

We now compare the above bounds with the exact unequal time anti-commutator at the non-interacting point $U=0$ in a $d$-dimensional square lattice, to investigate how tight these bounds are. The Hamiltonian is diagonalized as

[TABLE]

where $\omega(\vec{k})=2J\cos\vec{k}$. Therefore

[TABLE]

where $J_{n}(t)$ is the Bessel function, and we define $J_{\vec{n}}(t)\equiv\prod^{d}_{m=1}J_{n_{m}}(t)$. Inserting Eq. (84) into the LHS of Eq. (75), we get

[TABLE]

Comparing this with our bound in Eq. (75) [notice that the modified fermion method Eq. (78) also leads to the same result when $U=0$], we see that at the non-interacting point $U=0$, our LR bound just replaces the Bessel function $J_{\vec{x}_{j}-\vec{x}_{l}}(2Jt)$ of the exact solution by the modified Bessel function $I_{\vec{x}_{j}-\vec{x}_{l}}(2Jt)$. Since both $J_{x}(t)$ and $I_{x}(t)$ have the same asymptotic form $\frac{1}{x!}(\frac{t}{2})^{x}$ when $x/t$ is large, we conclude that our LR bound Eq. (75) has the tightest large-$x$ and small-$t$ exponent at the non-interacting point, where the large-$x$ exponent $\zeta$ of a bound $f(x,t)$ is defined as $\sup\{\zeta|\lim_{x\to\infty}e^{\zeta x\ln x}f(x,t)\text{ exists, }\forall t>0\}$, [Our bounds Eqs. (75) and (82) and the exact solution Eq. (85) all have $\zeta=1$].

V.5 Wen’s quantum rotor model

In this section we consider Wen’s quantum rotor model Wen (2003a, b); Levin and Wen (2006). The model is defined on a $d\geq 2$ dimensional hypercubic lattice whose sites are labeled by $\vec{r}$ and edges are labeled by $(\vec{r}\alpha)\equiv\langle\vec{r},\vec{r}+\vec{e}_{\alpha}\rangle$, $1\leq\alpha\leq d$. On each edge $(\vec{r}\alpha)$ there is a 2D quantum rotor, with angle variable $\hat{\theta}_{\vec{r}\alpha}$ and the dual angular momentum $\hat{S}_{\vec{r}\alpha}$, satisfying the canonical quantization condition $[\hat{\theta}_{\vec{r}\alpha},\hat{S}_{\vec{r}^{\prime}\beta}]=i\delta_{\vec{r},\vec{r}^{\prime}}\delta_{\alpha\beta}$. The Hamiltonian is

[TABLE]

where $\hat{W}_{\vec{r}\alpha\beta}=\exp(i\hat{\theta}_{\vec{r}\alpha}+i\hat{\theta}_{\vec{r}+\vec{e}_{\beta},\alpha}-i\hat{\theta}_{\vec{r}\beta}-i\hat{\theta}_{\vec{r}+\vec{e}_{\alpha},\beta})$. Notice that $\hat{W}_{\vec{r}\alpha\beta}=\hat{W}_{\vec{r}\beta\alpha}^{\dagger}$.

Based on some semiclassical treatments, Refs. Wen (2003a, b); Levin and Wen (2006) have shown that in the topologically-ordered phase $J\ll g$, this system hosts bosonic excitations propagating at speed $\sqrt{2gJ}$, behaving like artificial light waves. Ref. Hamma et al. (2009) studied a three-state truncated version of this model, and proved an LR bound with velocity $e\sqrt{2gJ}$. However, their method would fail in the untruncated version Eq. (86), since the term $(\hat{S}_{\vec{r}\alpha})^{2}$ has infinite operator norm $\|\hat{S}_{\vec{r}\alpha}^{2}\|=S^{2}$ (the number of states on a single site is $2S+1$), leading to an LR speed that diverges linear in $S$. Understanding what happens as $S\rightarrow\infty$ is important, since the argument that light emerges is valid only in this limit. In the following we use a slightly modified version of our method to overcome this difficulty, and prove that $v_{\mathrm{LR}}\leq 2X_{0}\sqrt{(d-1)gJ}$.

We begin by writing down the Heisenberg equations for $\hat{S}_{\vec{r}\alpha},\hat{W}_{\vec{r}\alpha\beta}$

[TABLE]

where we used $[\hat{S}_{\vec{r}\alpha},e^{\pm i\hat{\theta}_{\vec{r}^{\prime}\beta}}]=\pm e^{\pm i\hat{\theta}_{\vec{r}^{\prime}\beta}}\delta_{\vec{r},\vec{r}^{\prime}}\delta_{\alpha\beta}$. Taking the commutator with an arbitrary local operator $\hat{B}$, we have

[TABLE]

Now we get rid of the first term in the second equation by doing a change of variable $\hat{W}^{B}_{\vec{r}\alpha\beta}(t)=\hat{U}_{\vec{r}\alpha\beta}(t)\tilde{W}^{B}_{\vec{r}\alpha\beta}\hat{V}_{\vec{r}\alpha\beta}(t)$, where $\hat{U}_{\vec{r}\alpha\beta}(0)=\hat{V}_{\vec{r}\alpha\beta}(0)=\hat{I}$, $i\partial_{t}\hat{U}_{\vec{r}\alpha\beta}(t)=J[\hat{S}_{\vec{r}\beta}(t)+\hat{S}_{\vec{r}+\vec{e}_{\alpha},\beta}(t)-(\alpha\leftrightarrow\beta)]\hat{U}_{\vec{r}\alpha\beta}(t)$, and $i\partial_{t}\hat{V}_{\vec{r}\alpha\beta}(t)=\hat{V}_{\vec{r}\alpha\beta}(t)J[\hat{S}_{\vec{r}\beta}(t)+\hat{S}_{\vec{r}+\vec{e}_{\alpha},\beta}(t)-(\alpha\leftrightarrow\beta)]$. Then using the same derivations in Eqs. (7-10), we can prove that $\|\hat{S}^{B}_{\vec{r}\alpha}\|$ and $\|\hat{W}^{B}_{\vec{r}\alpha\beta}\|$ are upper bounded by the solution to the differential equation:

[TABLE]

with initial condition $S_{\vec{r}\alpha}(0)=\|\hat{S}^{B}_{\vec{r}\alpha}(0)\|,W_{\vec{r}\alpha\beta}(0)=\|\hat{W}^{B}_{\vec{r}\alpha\beta}(0)\|$. The maximal eigenfrequency at $\vec{k}=i\vec{\kappa}$ is $\omega_{m}(i\vec{\kappa})=4\sqrt{(d-1)gJ}\cosh(\kappa/2)$, therefore, using Eq. (31) we get

[TABLE]

This is our second example (after the large-$S$ Heisenberg XYZ model) of a finite LR speed in a system with an infinite local Hilbert space, and a second example (after TFIM) of an LR speed that grows sublinearly with spatial dimension $v\propto\sqrt{d}$. It is also a good bound quantitatively–in 2D, it is only larger than the semiclassical result $v=\sqrt{2gJ}$ by a factor of $\sqrt{2}X_{0}\approx 2.1$ 333Notice that the semiclassical result uses Cartesian distance instead of graph theoretical distance. If we use Cartesian distance as well, then the LR speed in all directions is upper bounded by $2\sqrt{gJ}\min_{\kappa>0}\frac{\sqrt{3+\cosh\kappa}}{\kappa}\approx 2.396\sqrt{gJ}$, which is only a factor of 1.69 larger than semiclassical result..

VI Improving the Lieb-Robinson velocity by eliminating large coupling constants

In several examples discussed in the previous section, our commutativity graph method dramatically improves the LR velocities when the local Hilbert space dimension or spatial dimension becomes large. Most prominently, in the large spin $S\to\infty$ limit of Heisenberg XYZ model and Wen’s quantum rotor model, our method gives finite LR velocities which removes the unphysical divergence in previous bounds. Nevertheless, there is still another type of unphysical divergence present in these bounds. This happens when a Hamiltonian is a sum of terms, and one set of mutually commuting terms is multiplied by coefficients that become very large. For example, in the large-$J$ or large-$h$ limits of TFIM in Eq. (49), the large-$U$ limit of SU($N$) FH model in Eqs. (76) and (81), and in the large-$J$ or large-$g$ limit of Wen’s quantum rotor model for any finite spin $S$, the bounds for the LR velocity diverge, while in reality, the actual velocities of information propagation in these limits are expected to be finite. This kind of divergence 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 removes this second type of unphysical divergence and gives a finite LR velocity in these limits. The basic idea is to derive LR bounds from coupled differential equations on a 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 Hastings and Koma (2006); Foss-Feig et al. (2015); Else et al. (2020), they were only used to remove single site terms, while our method can be used to remove arbitrary mutually commuting terms in the Hamiltonian. In the following we first present the general technique in Sec. VI.1, then in Sec. VI.2 to Sec. VI.4 we apply this technique to the specific examples mentioned above, and finally present a general result in Sec. VI.5.

VI.1 Eliminating large coupling constants on the commutativity graph

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\subset V_{G}$ denote a subset of disjoint vertices of $G$, such that $\{h_{j}|j\in 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 $\{\hat{\gamma}_{j}|j\in F\}$ mutually commute. We now construct a reduced graph $G^{\prime}_{F}=(V^{\prime},E^{\prime})$ from $G$ in the following way: (1) the vertices $V^{\prime}=V_{G}\backslash F$ of $G^{\prime}_{F}$ are obtained from $V_{G}$ by removing elements of $F$; (2) $E^{\prime}$ is defined such that for $\forall i,j\in V^{\prime}$, $\langle ij\rangle\in E^{\prime}$ if and only if either $\langle ij\rangle\in E_{G}$ or there exists an $l\in F$ such that $\langle il\rangle\in E_{G},\langle lj\rangle\in E_{G}$. In the left (right) panel of Fig. 6 we show the graph $G^{\prime}_{F}$ after removing parameter $J$ ($h$), respectively.

Let us denote by $\hat{I}$ the sum of all Hamiltonian terms corresponding to vertices in $F$ (i.e. the unwanted terms):

[TABLE]

The important step here is to consider the evolution of the operator $\hat{\gamma}_{i}$ (for $i\in G^{\prime}_{F}$) in a way that the unwanted term $\hat{I}$ is “rotated away” by a unitary transformation:

[TABLE]

Notice that $\hat{\gamma}_{i}(t,t)=\hat{\gamma}_{i}(t)$, the operator in Heisenberg picture. Taking derivative with respect to $s$, we have

[TABLE]

where the sum is over all terms in $\hat{H}-\hat{I}$ that do not commute with $\hat{\gamma}_{i}(0,t-s)$. Eq. (93) should be considered as the analog of Eq. (5). Now the same derivations as in Eqs. (6-II.2) for the commutator $\hat{\gamma}^{B}_{i}(s,t)=[\hat{\gamma}_{i}(s,t),\hat{B}(0)]$ result in

[TABLE]

In this case the generalized Grönwall’s inequality Chandra and Davis (1976) indicates that $\|\hat{\gamma}^{B}_{i}(t)\|$ is upper bounded by the solution $\gamma^{B}_{i}(t)$ to the differential equation

[TABLE]

with initial value $\gamma^{B}_{i}(0)=2\|\hat{B}\|$. Eq. (95) is the same as Eq. (10) with $G$ replaced by $G^{\prime}_{F}$.

In summary, when the parameters of a set of disjoint vertices of $G$ become very large, we can get a better bound for the speed of information propagation by solving the differential Eq. (95) on the reduced graph $G^{\prime}_{F}$, which has the large parameters removed, and links added when two vertices of $G_{F}^{\prime}$ 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 reduced graph $G^{\prime}_{F}$, respectively 444Notice however, that if the methods of this section are used, then the large-$x$ and small-$t$ exponents of the resulting bounds would not be tight anymore. For example, in the TFIM, both the large-$J$ bound and the large-$h$ bound gives small-$t$ exponent $\eta(\vec{r})=|\vec{r}|$, which is looser than the bound in Sec. V.1 and exact result by a factor of 2. Since a tighter LR speed is more important in typical applications, it is generally worthwhile to tighten the LR speed at the cost of compromising the large-$x$ and small-$t$ asymptotic behaviors..

VI.2 Example: TFIM at large $J$ or large $h$

As a first example, in the case of the TFIM, the left (right) panel of Fig. 6 shows the graph $G_{F}^{\prime}$ 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_{\frac{d-1}{d}}dJ$). Therefore, combined with Eq. (49), we have the upper bound for the LR speed

[TABLE]

for the $d$-dimensional TFIM. This is a major improvement compared to the previous best bound $v_{\mathrm{Ising}}\leq 8edJ$ (calculated from the method in Ref. Hastings (2010)), as we not only removed unphysical divergence of $v_{\mathrm{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$. [Notice however, that in 1D, a special method is available in a recent paper Ref. Chen and Lucas (2019b), which gives a bound for LR speed $2eJ\approx 5.44J$ for the 1D TFIM, which is about 10% tighter than our bound $4X_{0}J\approx 6.04J$ in the large-$h$ limit. Away from the large-$h$ limit, when $h/J<3.24$, our bound $\min\{2X_{0}\sqrt{Jh},4X_{0}h\}$ remains the tightest. Furthermore, it is possible to combine their methods and our method of commutativity graph, which will tighten all the three velocity bounds for 1D TFIM in Eq. (96) by roughly 10%. The result will be $v_{\text{LR}}\leq\min\{e\sqrt{Jh},2eJ,2eh\}$. Notice also that if we use Theorem 4 of Ref. Chen and Lucas (2019b) in our commutativity graph, we may be able to tighten our velocity bounds in higher dimensions as well. However, in $d>1$ it may be a hard combinatorical problem to compute the sum over all irreducible paths in Theorem 4 of Ref. Chen and Lucas (2019b), so we leave this as a future direction.]

VI.3 Example: SU($N$) FH model at large $U$

As a second example, we consider the SU($N$) FH model. Here, the tightest large-$U$ bound is calculated from the commutativity graph $G_{F}^{\prime}$ which is constructed from $G$ in Fig. 7 by eliminating the interaction term (which is $\propto U$) according to the prescription of this section. The result for the LR speed is $v_{\mathrm{FH}}\leq 4X_{0}NJ$. Combined with the results of Sec. V.4, we have the upper bound

[TABLE]

for the 1D SU($N$) FH model, which is also a significant improvement of the previous best bound $v_{\mathrm{FH}}\leq 8eNJ$, see Fig. 5 for the comparison.

VI.4 Example: Wen’s quantum rotor model at large $J$ or large $g$

As a third example of eliminating large parameters, we consider the spin-$S$ truncated version of Wen’s quantum rotor model, whose Hamiltonian has the same form as Eq. (86), with $\hat{S},e^{i\hat{\theta}}$ (on each site) replaced by the following $(2S+1)\times(2S+1)$ matrices:

[TABLE]

where $I_{2S}$ is the $2S\times 2S$ identity matrix. Ref. Hamma et al. (2009) studied the 2D $S=1$ case, and proved an LR bound with velocity $e\sqrt{2gJ}$, in the topologically-ordered phase $J\ll g$. Here we will prove that in this limit $J\ll g$ the LR speed actually has a much tighter upper bound $v\leq\mathrm{const}.\times JS^{2}$. For illustrative purpose, we will focus on the 2D case for simplicity, but the generalization to arbitrary spatial dimension is straightforward. To prove this, we first notice that the commutativity graph of the Hamiltonian terms $\hat{S}_{\vec{r}\alpha}^{2}-S^{2}/2$ and $\hat{W}_{\vec{r}\alpha\beta}+\mathrm{H.c.}$ is exactly the same with that of 2D TFIM shown in Fig. 3, where the blue circles represent the $S$ terms, with parameter $JS^{2}/2$, and the red triangles represent the $W$ terms, with parameter $g$. Therefore, we can borrow the result from the TFIM with substitutions $J\to JS^{2}/2,~{}h\to g$, and we get $v_{\mathrm{LR}}\leq\min\{4X_{\frac{1}{2}}JS^{2},8X_{0}g\}$, which is much tighter than the result of Ref. Hamma et al. (2009) $e\sqrt{2gJ}S$ in the small $J$ limit for any finite $S$. Combined with the result of Sec. V.5, we obtain

[TABLE]

Our bound indicates that the semiclassical treatment of the rotor model would fail for any finite $S$ for sufficiently small $J$ (with fixed $g$), since the semiclassical approximation predicts emergent light traveling at speed $\sqrt{2gJ}$ which violates our bound for sufficiently small $J$. In the $S\to\infty$ limit, however, our large-$J$ bound diverges, and the method in Sec. V.5 has to be used, which gives $v\leq\mathrm{const.}\sqrt{gJ}$, agreeing with the semiclassical result.

VI.5 A general result on perturbed solvable models with mutually commuting terms

We end this section by making a general statement about a family of perturbed exactly solvable models, whose Hamiltonian has the form

[TABLE]

where $\hat{H}_{0}$ is the solvable part and $J\sum_{j}\hat{V}_{j}$ is a sum of local perturbations. Here $j$ is a label for terms in the Hamiltonian (it does not necessarily mean a single site: $\hat{V}_{j}$ can include multi-body interactions). We make the following assumptions:

(1) $\hat{H}_{0}=\sum_{i}\hat{h}_{i}$ is the sum of mutually commuting, locally-interacting terms, i.e. $[h_{i},h_{j}]=0,~{}\forall i,j$. Primary examples include the $J$ term or the $h$ term of TFIM, the interaction term of BH model or SU($N$) FH model, the $J$ term or the $g$ term of the spin-$S$ truncated Wen’s quantum rotor model, the string-net models Levin and Wen (2005) with Kitaev’s toric code model Kitaev (2003) being a prominent special case, and the exactly solvable fracton models Haah (2011); Vijay et al. (2016); Nandkishore and Hermele (2019).

(2) There is a uniform upper bound for the norm of local perturbations $\|\hat{V}_{j}\|\leq c,\forall j$. This is true for any translation invariant system with a finite local Hilbert space.

Then we have

Theorem 1**.**

There exist constant $c^{\prime}>0$ such that $v_{\mathrm{LR}}\leq c^{\prime}|J|$, i.e. the LR speed of the model Eq. (100) grows at most linearly in the strength of the perturbation.

Notice that the second condition is necessary, since in Wen’s quantum rotor model with infinite $S$, in the limit $J\ll g$, the semiclassical treatment gives artificial light propagating at speed $v\propto\sqrt{Jg}$, which does not satisfy $v\leq c^{\prime}|J|$ for sufficiently small $J$. Theorem 1 can be proved by using the prescription of this section to eliminate all the terms in $\hat{H}_{0}$, the reduced graph $G^{\prime}_{F}$ would only involve vertices representing local perturbations $\hat{V}_{i}$. Locality of the original Hamiltonian $\hat{H}$ implies a uniform upper bound on the degree of every vertex in $G^{\prime}_{F}$, which, combined with condition (2), yields an LR bound with $v_{\mathrm{LR}}\leq c^{\prime}|J|$, based on similar treatment as in Eq. (39). The previous bounds for the Hamiltonian Eq. (100) give an LR speed that stays finite in the limit $J\to 0$ unless $\hat{H}_{0}$ only contains on-site terms (which is not the case for models with topological order), since the previous $v_{\mathrm{LR}}$ is proportional to the operator norm of local Hamiltonians $\hat{h}_{a}$ in the solvable part $\hat{H}_{0}$. Our bound, which vanishes linearly in perturbation strength $J$, is a qualitative improvement.

VII Consequences of tighter 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 Bravyi et al. (2006); Gogolin and Eisert (2016); Eisert et al. (2015), studying equilibrium properties like exponential decay of connected correlation functions Hastings (2004a); Hastings and Koma (2006); Nachtergaele and Sims (2006); Cramer and Eisert (2006); Kliesch et al. (2014), quantifying the entanglement area law Hastings (2007b), and providing error bounds on numerical algorithms Osborne (2006a, 2007a, 2007b); Woods and Plenio (2016); Tran et al. (2019). 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 Hastings (2004a); Hastings and Koma (2006); Nachtergaele and Sims (2006); Cramer and Eisert (2006).

We will achieve this goal by directly replacing the previous LR bound by our bound Eq. (3) in the proof of Ref. Hastings and Koma (2006), and obtain a tighter upper bound on the correlation length $\xi$. Our bound is not only a quantitative improvement, but are qualitatively tighter in several limits, e.g. in the $J/h\to\infty$ limit of TFIM and the $S\to\infty$ limit of Heisenberg XXZ, where our bounds have the same asymptotic scaling as the known exact solutions.

Let us start with the following well-established inequality which relates the ground state correlation to the unequal time commutator Hastings and Koma (2006); Hastings (2010)

[TABLE]

where $c_{0}$ is a constant independent of $\alpha$ and $t$, and we assumed $\langle G|\hat{A}_{X}|G\rangle=\langle G|\hat{B}_{Y}|G\rangle=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 LR bound in Eq. (3), while for $t\geq r/u$ we use the trivial bound $\|[\hat{A}_{X}(t),\hat{B}_{Y}]\|\leq 2\|\hat{A}\|\|\hat{B}\|\equiv c_{1}$. We get

[TABLE]

where $r=d_{XY}$, $C$ is a constant defined in Eq. (42), and in the second step we applied the inequality $e^{-\alpha t^{2}}\leq e^{-\alpha(r/u)^{2}+2\alpha(r/u)t}$ for $t\geq r/u$. If we choose $\alpha^{-1}=\lambda r$, then all three terms above decay exponentially with distance $r$, and the first term is upper bounded by $\frac{c}{\pi u}\max_{0\leq\tau\leq 1/u}e^{-r\tau^{2}/\lambda}u^{r}\tau^{r-1}$. Choosing $\lambda$ to maximize the smallest decay coefficient of the three terms, the lower bound for $\xi^{-1}$ is

[TABLE]

where $W(x)$ is the Lambert $W$ function (the product logarithm) defined as the solution to the equation $We^{W}=x$.

As a first specific example, we applied this bound to the 1D TFIM, and in Fig. 8 we give a comparison between our bound with the previous best bound $\xi^{-1}\geq\max_{\mu>0}\frac{\mu}{1+16Je^{\mu}/\Delta}$ given in Ref. Hastings and Koma (2006), and also with the exact solution $\xi^{-1}\geq|\ln(J/h)|$ given in Ref. Pfeuty and Elliott (1971) 555Theorem 2 of Ref. Hastings (2010) upper bounds the ground state correlation by $c_{1}\exp(\frac{r\Delta}{2v_{\text{LR}}})+c_{2}g(r)$, where $c_{1},c_{2}$ are some positive constants, and $g(r)$ first appeared in Eq. (30) of the same paper, which, as Ref. Hastings (2010) claims, “decays faster than exponential in finite range interacting systems”. If this claim is correct, one would have the bound $\xi^{-1}\geq\frac{\Delta}{2v_{\text{LR}}}$, which is even tighter than our bound in Eq. (103) at large $\Delta/{v_{\text{LR}}}$. But this bound is actually incorrect. It is true that in locally interacting systems, the LR bound $C_{B}(X,t)$ decays faster than exponential at fixed time, as demonstrated in our paper. However, when the bound is written in the form Eq. (30) of Ref. Hastings (2010), the function $g(r)$ only decays exponentially. This can be easily verified in the free fermion solution Eq. (85) in our paper, by noticing that the Bessel function $J_{x}(x/w)$ decays exponentially with $x$ for arbitrary $w>1$, see Ref. Watson (1922). Also, in the case of 1D TFIM, the bound $\xi^{-1}\geq\frac{\Delta}{2v_{\text{LR}}}$ would imply that $\xi^{-1}$ grows at least linearly with $h/J$ as $h/J$ becomes large, while the exact solution only shows a logarithmic growth, see Fig. 8 .. As we see, in this case our bound is always tighter than the previous bound. This is an especially strong improvement in the limit $J\gg h$ where the previous bound gives a constant, while the bound in Eq. (103) gives $\xi\propto\log(J/h)$, in agreement with the exact solution.

As a second example, we apply the bound Eq. (103) to the large-$S$ Heisenberg XXZ model in 1D. Our bound for the LR speed in Eq. (63) gives $u=v\leq 2X_{1}(J+\sqrt{J^{2}+8J_{z}J})$, where $J_{x}=J_{y}=J$. Therefore, our bound for $\xi$ remains finite in the limit $S\to\infty$ provided that the spectral gap $\Delta$ is finite, while the bound given by previous method Hastings and Koma (2006) $\xi\leq(\text{const.})\times S/\Delta$ diverges linearly in $S$.

Finally, let us discuss the implication of Theorem 1 from Sec. VI.5 together with the bound in Eq. (103), in the case of perturbed exactly solvable topologically ordered models with a finite local Hilbert space dimension, such as Kitaev’s toric code model Kitaev (2003), string-net models Levin and Wen (2005), the X-cube model Vijay et al. (2016), and Haah’s cubic code model Haah (2011). When these exactly solvable commuting Hamiltonians are perturbed with $J\sum_{i}\hat{V}_{i}$, in the limit $J\to 0$, the spectral gap $\Delta$ stays finite, while Theorem 1 implies that the LR velocity on the commutativity graph has an upper bound $u\leq c^{\prime}J\ll\Delta$. Therefore, Eq. (103) gives the bound

[TABLE]

i.e. when the perturbation is turned off $J\to 0$, the correlation length $\xi$ vanishes at least as fast as $|\ln J|^{-1}$. The previous LR bounds for this family of perturbed exactly solvable models always give a finite LR velocity when $J\to 0$, so the previous upper bounds on $\xi$ would also be finite when $J\to 0$, qualitatively looser than our bound.

VIII Conclusions

We introduce a method to obtain tighter LR bounds by taking into account the details of the Hamiltonian, and show that this method improves LR bounds both qualitatively and quantitatively. The bounds established in this paper have much tighter LR speeds, especially in three broad scenarios: (1) When the number of local degrees of freedom becomes large, for example in the large-spin TFIM and Heisenberg XYZ model, the $N$-state truncated BH model, the large-$N$ SU($N$) FH model, and Wen’s quantum rotor model. (2) When a set of commuting parameters becomes large, for example in the TFIM at large $J$ and the perturbed toric code model. (3) In large spatial dimension, for example in the TFIM and Wen’s quantum rotor model. In these limits our LR speeds have qualitatively better asymptotic scaling than previous best bounds, and in some cases [e.g. Wen’s quantum rotor model and the classical limit of spin-$S$ Ising, Heisenberg and SU($N$) FH] our bounds are qualitatively tight in the sense that they scale in the same way as the known semiclassical solution. Our bounds also have superexponential decaying tails $e^{-d_{XY}\ln d_{XY}}$ as well as a much tighter short time scaling $O(t^{d_{XY}})$, instead of previous bounds’ exponential decay $e^{-\mu d_{XY}}$ and $O(t)$ scaling. As demonstrated in several examples, the large distance and small time asymptotic behaviors of our bounds are often tightest possible.

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.1 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 $c\left(\frac{u|t|}{d_{XY}}\right)^{d_{XY}}$, 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 a wide range of examples 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 bound leads to a tighter bound on the ground state correlation length in gapped systems, which qualitatively improves previous best bounds in the large-$J$ limit of the 1D TFIM and large-$S$ limit of 1D spin-$S$ Heisenberg XXZ model, 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 Bose (2007); Epstein and Whaley (2017), or schemes that create certain entangled states Clark et al. (2005); Monz et al. (2011); Gao et al. (2013). Since the LR bounds can be used to bound the time needed for these dynamical processes Bravyi et al. (2006), 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 bounds 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. Finally, as we demonstrated at the end of Sec. VI.4, our bounds are tight enough that they may be used as a rigorous validity check for the various approximation methods used in quantum many-body physics, such as the semiclassical (truncated Wigner) approximation Sinatra et al. (2002); Polkovnikov (2010); Davidson et al. (2017), mean-field approximations Altland and Simons (2010); Kamenev (2011), and random phase approximation Negele (2018). An approximation is immediately proved invalid if it predicts a propagation speed $v$ that exceeds our bound, or if it gives ground state correlation length $\xi$ and spectral gap $\Delta$ that does not satisfy the inequality in Eq. (103).

The methods introduced in this paper may also shed light on some of the important unsolved problems in the area of LR bounds. We mention two of them below.

One important open problem is how to extend LR bounds to systems with an infinite local Hilbert space dimension, such as in systems with interacting bosons. Despite extensive effort, up to now there are only some preliminary success in obtaining a finite LR bound in harmonic systems Cramer et al. (2008); Nachtergaele et al. (2009), or bounds for a restricted set of operators and special initial states Schuch et al. (2011) in BH-type models. Our method is very promising in this direction, as already demonstrated in this paper, it leads to a finite LR bound for large-$S$ Heisenberg XYZ model and Wen’s quantum rotor model, and has shown qualitatively better scaling laws of $v_{\mathrm{LR}}$ with number of local states in truncated BH model and SU($N$) FH model, all of which have not been achieved before. In future, we hope that by extending the methods of this paper, one can find a bigger class of nontrivial interacting models with infinite local Hilbert space whose LR speeds remain finite.

Another open problem is to tighten the LR bounds in systems with power-law decaying interactions Foss-Feig et al. (2015). The generalization of our method to this case is straightforward. It may not lead to a qualitative improvement in the most generic case, but when applied to specific models, our method can take advantage of the specific properties of the Hamiltonian to give improved bounds. A promising case is when the long-range interactions are mutually commuting. This includes Coulomb interactions, van der Waals interactions between Rydberg atoms, the $\sigma^{z}_{i}\sigma^{z}_{j}$ interaction in power-law TFIM, all of which are of great theoretical and experimental interests. In Sec. V we have seen in many examples that our commutativity graph method has a qualitative advantage over previous methods especially when there are lots of commuting terms that have a spatial overlap. It is therefore natural to expect that if one combines our commutativity graph method with the current best methods Tran et al. (2019); Chen and Lucas (2019a); Kuwahara and Saito (2020) in treating the long-range part, one can get qualitatively better, model-specific LR bounds in a family of models that are of great physical significance.

Acknowledgements.

We thank Stefano Chessa, Yichen Huang, Andrew Lucas, Michael Foss-Feig, Putian Li, Bhuvanesh Sundar, and Ian White for discussions. This work was supported by the Welch Foundation (C-1872) and the National Science Foundation (PHY-1848304).

Appendix A A more general method

In the examples of Sec. V, there are two places where our method slightly deviates from the general methods introduced in Sec. II: in the SU($N$) FH model Eqs. (77-78) we used a method which is essentially a mixture of the commutativity graph method and the fermion method, and in Wen’s quantum rotor model Eq. (V.5) we started with the Heisenberg equation for $\hat{S}^{z}_{ij}$ and $\hat{W}_{i}$ rather than the Hamiltonian terms $(\hat{S}^{z}_{ij})^{2}$ and $\hat{W}_{i}$ as we did in Sec. II.1. One may therefore be curious whether there is a more general prescription behind these seemingly different approaches. Indeed, as we will see in this section, the commutativity graph method, the fermion method, the modified fermion method in the SU($N$) FH model, and the method for Wen’s quantum rotor model, are just special cases of a more general method to be introduced below. This general method starts with the Heisenberg equation for an arbitrary set of local operators, thereby allowing a high level of flexibility. The downside of this flexibility is that in some cases one may need considerable physical and mathematical intuition to make the best choice. The general prescription is:

Step 1. Choose a local operator basis $\{\hat{\Omega}_{i}\}_{i\in M}$ for the Hamiltonian $\hat{H}$, here $M$ is a collection of indices. By basis, we mean that the Hamiltonian can be written as a sum of products of elements of this basis in at least one way (there can be more than one way, in this case we call the basis overcomplete)

[TABLE]

where $I$ denotes an ordered subset of $M$. For example, in the commutativity graph method, the basis is simply the collection of all individual terms in the Hamiltonian Eq. (4); in the fermion method, the basis is the set of all Majorana operators; in the modified fermion method of FH model, the basis is the set of all Majorana operators and all the $\hat{u}_{j\sigma\sigma^{\prime}}$ terms in Eq. (77), however, since $\hat{u}_{j\sigma\sigma^{\prime}}$ is a quartic product of Majorana operators, this basis is overcomplete; in Wen’s quantum rotor model, the basis is the collection of all $\hat{S}^{z}_{ij}$ and $\hat{W}_{i}$. One can easily verify that in all these cases, the Hamiltonian can be written as a polynomial of the basis elements.

Step 2. Write down the Heisenberg equation of motion for the basis operators, as is done in Eqs. (5, 21, 77, V.5). For the generic Hamiltonian in Eq. (105), we have

[TABLE]

where in the second line we apply the basic commutation relations of the basis operators to simplify the equation, and in the third line we isolate all the terms that are left-proportional or right-proportional to $\hat{\Omega}_{i}$, here $P_{i},Q_{i},R_{i}$ are some polynomial functions of the basis operators, and $P_{i}[\{\hat{\Omega}_{j}\}],Q_{i}[\{\hat{\Omega}_{j}\}]$ are assumed to be Hermitian (we can always absorb the non-Hermitian part into $R_{i}[\{\hat{\Omega}_{j}\}]$).

Step 3. On both sides of the Heisenberg equation, take the commutator with another local operator $\hat{B}$ (note: if both the LHS and $\hat{B}$ are fermionic operators, anti-commutator should be taken). Use Leibniz’s rule to expand the product terms, e.g. $[\hat{a}\hat{b}\hat{c},\hat{B}]=[\hat{a},\hat{B}]\hat{b}\hat{c}+\hat{a}[\hat{b},\hat{B}]\hat{c}+\hat{a}\hat{b}[\hat{c},\hat{B}]$. This is done in Eqs. (6, 22, V.5) [and below Eq. (77), implicitly] of the main text. We write the equation in the form

[TABLE]

where $S_{i}[\{\hat{\Omega}_{j}\},\{\hat{\Omega}_{j}^{B}\}]$ is some polynomial function that is linear in $\hat{\Omega}_{j}^{B}$.

Step 4. Use a change of variable $\hat{\Omega}^{B}_{i}(t)\to\hat{U}_{i}(t)\hat{X}^{B}_{i}(t)\hat{V}_{i}(t)$ to cancel the first two terms on the RHS of Eq. (107), where $i\partial_{t}\hat{U}_{i}(t)=P_{i}[\{\hat{\Omega}_{j}\}]\hat{U}_{i}(t),i\partial_{t}\hat{V}_{i}(t)=\hat{V}_{i}(t)Q_{i}[\{\hat{\Omega}_{j}\}]$ 666One can actually get rid of more general terms of the form $\sum_{a}\hat{P}_{a}\hat{\Omega}^{B}_{i}\hat{Q}_{a}$ for arbitrary Hermitian operators $\hat{P}_{a},\hat{Q}_{a}$ (which are implicitly functions of the basis operators $\{\hat{\Omega}_{j}\}$) on the RHS of Eq. (107), i.e. on the RHS of Heisenberg equation for $i\partial_{t}\hat{\Omega}^{B}_{i}$ one can get rid of any product term that is proportional to $\hat{\Omega}^{B}_{i}$ with coefficients being Hermitian operators. This can be done by going to the Liouville picture, where Eq. (107) can be rewritten as $i\partial_{t}|\Omega^{B}_{i})=\sum_{a}\hat{\Pi}_{a}|\Omega^{B}_{i})+\sum_{j\neq i}\hat{S}_{j}|\Omega^{B}_{j})$, where the operator $\hat{\Omega}_{j}^{B}$ is mapped to a state $|\Omega^{B}_{i})$, $\hat{\Pi}_{a}$ denotes the superoperator representing the map $\hat{\Omega}^{B}_{i}\to\hat{P}_{a}\hat{\Omega}^{B}_{i}\hat{Q}_{a}$, and $\hat{S}_{j}$ denotes the superoperator representing the map $\hat{\Omega}^{B}_{j}\to S_{i}[\{\hat{\Omega}_{j}\},\{\hat{\Omega}^{B}_{j}\}]$ (remember that $S_{i}[\{\hat{\Omega}_{j}\},\{\hat{\Omega}^{B}_{j}\}]$ is linear in $\hat{\Omega}^{B}_{j}$). If $\hat{P}_{a},\hat{Q}_{a}$ are all Hermitian, then $\hat{\Pi}_{a}$ is Hermitian, and the first term can also be rotated away by a unitary rotation in the Liouville picture, which also preserves the norm of the state $|\Omega^{B}_{i})$. The steps 5 and 6 will follow without problem.. Eq. (107) becomes

[TABLE]

Notice also $\|\hat{\Omega}^{B}_{i}(t)\|=\|\hat{X}^{B}_{i}(t)\|$ due to $\hat{U}_{i},\hat{V}_{i}$ being unitary. This is done in Eq. (7) [and below Eq. (77), implicitly] and between Eqs. (V.5, V.5).

Step 5. Take operator norm on both sides of the resulting operator evolution equation (108), and apply the triangle inequalities on the RHS to get $\|S_{i}[\{\hat{\Omega}_{j}\},\{\hat{\Omega}_{j}^{B}\}]\|\leq\sum_{j}S_{ij}\|\hat{\Omega}_{j}^{B}\|$ for suitable non-negative coefficients $S_{ij}$. Then using the fundamental theorem of calculus and Grönwall’s inequality, one can prove that $\|\hat{\Omega}^{B}_{i}(t)\|$ is upper bounded by the solution to the system of linear differential equations

[TABLE]

with initial condition $\bar{\Omega}^{B}_{i}(0)=\|\hat{\Omega}^{B}_{i}(0)\|$. This is done in Eqs. (II.2-10, II.3-24, 78, V.5).

Step 6. Solve the resulting linear differential equations (109), and calculate $v_{\mathrm{LR}}$ from its maximal eigenfrequency using Eq. (31).

The main flexibility allowed by this general method is the choice of basis operators $\{\hat{\Omega}_{i}\}_{i\in M}$. A different choice of basis would typically lead to a different LR speed. We would hope to find a basis that makes the resulting LR speed as small as possible to make the LR bound tighter. We have already demonstrated in the main text that the basis $\{\hat{\gamma}_{i}\}$ used in the minimal Clifford decomposition Eq. (4) satisfying Eq. (II.2) is a good choice of basis for a large class of spin and interacting fermion models, and the Majorana operator basis $\{\hat{c}_{i}\}$ satisfying Eq. (II.3) is good for weakly interacting fermions, in the sense that the resulting LR speeds vastly improves the previous results, and the LR bounds typically have tight large distance and small time exponents. We did encounter a few exceptions where other choices are favorable. For example in Wen’s quantum rotor model we used the basis $\{\hat{S}^{z}_{ij},\hat{W}_{i},\hat{W}^{\dagger}_{i}\}$, and in the large-$N$ SU($N$) FH model we used $\{\hat{c}_{i\sigma},\hat{u}_{j\sigma\sigma^{\prime}}\}$. Also, in the spin-$S$ Heisenberg model, an alternative, simpler choice of basis $\{S^{x}_{i},S^{y}_{i},S^{z}_{i}\}$ leads to $v\leq 4dJ_{m}X_{1}$, which is the same as Eq. (63) in the $S\to\infty$ limit. It may be a good future direction to further explore the role of basis operators and hopefully obtain a general prescription on how to find the best choice of basis in a specific model.