Quantum Scrambling and State Dependence of the Butterfly Velocity

TL;DR

This paper investigates how operator growth and scrambling velocity in quantum many-body systems relate to the state dependence of out-of-time-ordered correlators, supported by theoretical bounds and numerical simulations.

Contribution

It establishes a bound on the state dependence of out-of-time-ordered correlators based on the scrambling velocity in local lattice models.

Findings

The bound is verified in the thermal mixed-field Ising spin chain.

The butterfly velocity exhibits a crossover from high-temperature to low-temperature values.

The low-temperature butterfly velocity differs significantly from the high-temperature value.

Abstract

Operator growth in spatially local quantum many-body systems defines a scrambling velocity. We prove that this scrambling velocity bounds the state dependence of the out-of-time-ordered correlator in local lattice models. We verify this bound in simulations of the thermal mixed-field Ising spin chain. For scrambling operators, the butterfly velocity shows a crossover from a microscopic high temperature value to a distinct value at temperatures below the energy gap.

Figures

Equations

$$\mathcal{C}(\mathbf{x},t;\rho )\equiv \mathrm{tr}\left(\rho [O_1(0,t),O_2(\mathbf{x},0){]}^{†}[O_1(0,t),O_2(\mathbf{x},0)]\right),$$

$$\mathcal{C}(\mathbf{x},t;\rho )\sim e^{-\lambda (|\mathbf{x}-{\mathbf{x}}_0|/v_B-t{)}^{1+p}/t^p}.$$

$$H=\sum _{\mathbf{x}\in \Lambda }{h}_{\mathbf{x}},$$

$$v_{\mathrm{LR}}\equiv \sup \left\{v:\underset{t\to \infty }{\lim }\frac{1}{t}\ln \Vert [O_1(0,t),O_2(vt\mathbf{n},0)]\Vert \ge 0\right\}.$$

$$\Vert [O_1(0,t),O_2(x\mathbf{n},0)]\Vert \le C_{\mathrm{LR}}\Vert O_1\Vert \Vert O_2\Vert e^{(vt-x)/{\xi }_{\mathrm{LR}}}.$$

$$\lambda (\mathbf{v};\rho )\equiv \underset{t\to \infty }{\lim }\frac{1}{t}\ln \mathcal{C}(\mathbf{v}t,t;\rho ).$$

$$v_B(\rho )\equiv \sup \left\{v:\lambda (v\mathbf{n};\rho )\ge 0\right\},$$

$$v_S(\mathbf{v})\equiv \underset{t\to \infty }{\lim }\frac{R(\mathbf{v}t,t)}{t}.$$

$$|\partial_{\beta}\lambda({\boldsymbol{v}};\rho)|\leq\frac{2h}{a}\Big{(}v_{S}({\boldsymbol{v}})-(\xi+\xi_{\text{LR}})\lambda({\boldsymbol{v}};\rho)\Big{)}\,,$$

$$\left|\frac{\partial\lambda({\boldsymbol{v}};\rho)}{\partial\mu_{i}}\right|\leq\frac{2c^{i}}{a}\Big{(}v_{S}({\boldsymbol{v}})-(\xi+\xi_{\mathrm{LR}})\lambda({\boldsymbol{v}};\rho)\Big{)}.$$

$${\partial }_{\beta }\mathcal{C}(\mathbf{x},t;\rho )=-\mathrm{tr}(\rho \tilde{H}{O}^{†}O)=-\mathrm{tr}(\tilde{H}\sqrt{\rho }{O}^{†}O\sqrt{\rho }),$$

$$\tilde{H}=\sum _{|\mathbf{y}-{\mathbf{y}}_0|\le R+\delta }{\tilde{h}}_{\mathbf{y}}+\sum _{|\mathbf{y}-{\mathbf{y}}_0|>R+\delta }{\tilde{h}}_{\mathbf{y}},$$

$$\sum _{|\mathbf{y}-{\mathbf{y}}_0|>R+\delta }\left|\mathrm{tr}({\tilde{h}}_{\mathbf{y}}\sqrt{\rho }{O}^{†}O\sqrt{\rho })\right|\le C^{\prime }\underset{\mathbf{y}\in \Lambda }{\sup }\Vert {\tilde{h}}_{\mathbf{y}}\Vert \Vert O{\Vert }^2{e}^{-\delta /\xi }.$$

$$h\equiv 2\underset{\mathbf{y}\in \Lambda }{\sup }\Vert h_{\mathbf{y}}\Vert .$$

$$|\mathrm{tr}({\tilde{h}}_{\mathbf{y}}\sqrt{\rho }{O}^{†}O\sqrt{\rho })|\le \Vert {\tilde{h}}_{\mathbf{y}}\Vert \mathrm{tr}(\rho O^{†}O)\le h\mathcal{C}(\mathbf{x},t;\rho ).$$

$$h=\underset{\mathbf{y}\in \Lambda }{\sup }\frac{|\mathrm{tr}({\tilde{h}}_{\mathbf{y}}\sqrt{\rho }{O}^{†}O\sqrt{\rho })|}{\mathrm{tr}(\rho O^{†}O)}.$$

$$|{\partial }_{\beta }\mathcal{C}(\mathbf{x},t;\rho )|\le V_{R+\delta }h\mathcal{C}(\mathbf{x},t;\rho )+C^{\prime }h\Vert O{\Vert }^2{e}^{-\delta /\xi }.$$

$$|{\partial }_{\beta }\lambda (\mathbf{v};\rho )|\le h\underset{t\to \infty }{\lim }\frac{V_{R-\xi \lambda (\mathbf{v};\rho )t}}{t}.$$

$$|\partial_{\beta}\lambda({\boldsymbol{v}};\rho)|\leq\frac{2h}{a}\Big{(}v_{S}({\boldsymbol{v}})-(\xi+\xi_{\text{LR}})\lambda({\boldsymbol{v}};\rho)\Big{)}\,.$$

$$O_1(t)=\sum _{n=0}^{\infty }\frac{(it[H,\cdot ]{)}^n}{n!}{O}_1=O_1+it[H,O_1]-\frac{t^2}{2}[H,[H,O_1]]+\dots ,$$

$$H=-J\sum _{i=1}^{N-1}{Z}_i{Z}_{i+1}+h_X\sum _{i=1}^N{X}_i+h_Z\sum _{i=1}^N{Z}_i,$$

$$\lambda (\mathbf{v};\rho )=-\lambda (v/v_B-1{)}^{1+p}\text{for}v\ge v_B.$$

$$\displaystyle a\lambda(\Delta v)^{1+p}\Big{|}\partial_{\mu_{i}}\ln\lambda+\ln(\Delta v)\partial_{\mu_{i}}p$$

$$\le 2c^i\left[v_S(v)+(\xi +{\xi }_{\mathrm{LR}})\lambda (\Delta v{)}^{1+p}\right],$$

$$v_S=0\Rightarrow {\partial }_{{\mu }_i}{v}_B=0.$$

$$h=\underset{t>0,\mathbf{y}\in \Lambda }{\sup }\frac{|\mathrm{tr}({\tilde{h}}_{\mathbf{y}}\sqrt{\rho }{O}^{†}O\sqrt{\rho })|}{\mathrm{tr}(\rho O^{†}O)},$$

$$\left|\ln \frac{v_B(\beta =\infty )}{v_B(\beta =0)}\right|\lesssim {\int }_0^{1/\Delta }d\beta \frac{2h{v}_B^{p/(1+p)}[1+(\xi +{\xi }_{\mathrm{LR}})\lambda /v_B]}{a\lambda (1+p)}{\left(v_S^B\right)}^{1/(1+p)},$$

$$v_B(T=0)\gtrsim v_{\text{LR}}{e}^{-\alpha v_S^{\gamma }/\Delta }.$$

$${\mathcal{P}}_S[O]\equiv {\int }_{\mathrm{supp}U\cap S=\varnothing }dUUO{U}^{†},$$

$$\Vert {\mathcal{P}}_S[O]\Vert \le \Vert O\Vert ,\Vert {\mathcal{Q}}_S[O]\Vert =\Vert O-{\mathcal{P}}_S[O]\Vert \le \Vert O\Vert +\Vert {\mathcal{P}}_S[O]\Vert \le 2\Vert O\Vert ,$$

$$\Vert {\mathcal{Q}}_{\mathbf{x}}^r[O]\Vert \le C\Vert O\Vert e^{-(r-R)/\xi }.$$

Full text

Strongly quantum many-body systems have been important in condensed matter [1, 2] and nuclear physics [3, 4] for some time and are likely to become increasingly important with the ongoing development of quantum information processing technology [5, 6, 7]. It is essential to understand the spatio-temporal dynamics of these systems in highly quantum regimes where semiclassical methods such as the Boltzmann equation are inapplicable.

Significant progress has been made recently by considering quantum scrambling in many-body systems [8, 9, 10, 11, 12, 13]. Quantum scrambling arises when operator growth under Heisenberg time evolution redistributes local information to non-local degrees of freedom. It has been found that scrambling in spatially local systems is characterized by both a rate and a velocity, e.g. [14, 15, 16, 17]. These universal properties are manifested in the so-called out-of-time-ordered correlator (OTOC):

[TABLE]

defined for local operators $O_{1},O_{2}$ in state $\rho$. The OTOC has been found to reveal a ‘light cone’ spread of quantum information, with two state-dependent characteristics: the quantum Lyapunov exponent $\lambda$ and the butterfly velocity $v_{B}$. Just outside the light cone (or ‘butterfly cone’) $|{\boldsymbol{x}}|\gtrsim v_{B}t$ for $t>0$, the OTOC grows as the front is approached according to [18, 19]:

[TABLE]

In systems with many local degrees of freedom (e.g. large $N$ systems) the exponent $p=0$ and the growth is exponential. This case is reminiscent of the classical butterfly effect. In spin lattice systems, generally $p>0$, so that the front broadens as it spreads.

The butterfly velocity is a state-dependent speed of information propagation that is universally present in local systems, plausibly controlling important physical processes such as transport in strongly quantum regimes [20, 21, 22, 23, 24, 25, 26]. The state dependence means that the butterfly velocity is a more powerful probe of dynamics than the widely employed microscopic Lieb-Robinson velocity [27]. In this work we will show that this state dependence (e.g. temperature dependence) is tied to the underlying quantum scrambling of operators.

In quantum field theories that describe a nontrivial (quantum critical) continuum limit of lattice systems, the scaling of the butterfly velocity with temperature is $v_{B}\sim T^{1-1/z}$ in the simplest cases [16, 20]. The dynamical critical exponent $z$ describes the relative scaling of space and time. In this work we will characterize the butterfly velocity in general lattice models, away from critical points and without a large $N$ limit. We will obtain the temperature dependence of the butterfly velocity in quantum spin systems, extending previous infinite temperature results [28, 19]. The temperature dependence of scrambling in classical spin systems has been recently discussed in [29].

In a spatially local system the growth of operators determines a ‘scrambling velocity’ $v_{S}$, defined in (8) below. Our first result (9) states that the change of the velocity-dependent Lyapunov exponent — defined shortly in (6) — with temperature is bounded by the scrambling velocity. This result is rigorous for one-dimensional systems and plausibly true more generally. We verify the bound in numerical simulations of the mixed-field Ising model, focusing on the temperature dependence of the butterfly velocity. In Fig. 2 below we see that the non-interacting transverse field model has a temperature-independent butterfly velocity whereas the velocity is temperature-dependent for the interacting mixed field models. In these curves, the butterfly velocity crosses over from a microscopic infinite-temperature value to a low-temperature value. The temperature scale of the crossover is set by the energy gap.

Three velocities from locality

It will be crucial to understand three different velocities that characterize spatially local quantum systems. Our results will tie these velocities together. The velocities emerge in any lattice $\Lambda$ of spins (or fermions) with a local Hamiltonian

[TABLE]

where $h_{{\boldsymbol{x}}}$ are operators localized near lattice site ${\boldsymbol{x}}$. Translation symmetry is not required.

Lieb-Robinson velocity

The Lieb-Robinson velocity defines an emergent ‘light-cone’ causality from local dynamics on a lattice [27]. It is a state-independent, microscopic velocity set by the magnitude of couplings in the Hamiltonian, and is insensitive to operator growth or lack thereof.

A convenient and powerful definition of $v_{\mathrm{LR}}$ is in terms of space-time rays. That is, consider an operator $O_{2}$ located along the ray ${\boldsymbol{x}}=vt\boldsymbol{n}$ (here $\boldsymbol{n}$ is a unit vector). At large times we can introduce a velocity-dependent exponent $\lambda(v)$ that determines the growth or decay of the norm of the commutator along the ray, $\|[O_{1}(0,t),O_{2}(vt\boldsymbol{n},0)]\|\sim e^{\lambda(v)t}$. Here $O({\boldsymbol{x}},t)$ denotes $O$ translated by a lattice vector ${\boldsymbol{x}}$ in space and a time $t$ with Heisenberg evolution, and $\|\cdot\|$ is the operator norm. The causal light cone defined by $v_{\text{LR}}$ is such that for all $v>v_{\mathrm{LR}}$ the norm decays exponentially at late times, so that $\lambda(v)<0$. Therefore we can define $v_{\text{LR}}$ as the largest velocity such that the norm does not decay along a ray:

[TABLE]

We shall not keep the dependence on direction $\boldsymbol{n}$ and operators $O_{1},O_{2}$ explicit.

For any $v>v_{\mathrm{LR}}$ there are ($v$-dependent) constants $\xi_{\mathrm{LR}},C_{\mathrm{LR}}>0$ such that for all $t,x>0$,

[TABLE]

Intuitively, inequality (5) states that for $v>v_{\text{LR}}$, the norm $\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|$ is exponentially small outside the ray $x=vt$, with a tail of length $\xi_{\mathrm{LR}}(v)$.

Butterfly velocity

The butterfly velocity is defined analogously to the Lieb-Robinson velocity, but using the OTOC instead of the operator norm of the commutator [10, 14]. It therefore depends on the quantum state $\rho$.

The ‘velocity-dependent Lyapunov exponent’ is defined by the late time growth or decay of the OTOC along a ray [18]:

[TABLE]

Analogously to the Lieb-Robinson case, the butterfly velocity can now be defined as

[TABLE]

which is state-dependent. The operator norm bounds the OTOC and hence $0\leq v_{B}(\rho)\leq v_{\mathrm{LR}}$.

Scrambling velocity

The Lieb-Robinson bound (5) implies that the size of an operator can grow at most polynomially in time (as $t^{d}$ in a $d$-dimensional system). In contrast, the growth can be exponential without spatial locality, such as in SYK models [30, 31, 32]. Operator growth under Heisenberg evolution in quantum systems with a local Hamiltonian will therefore define another velocity. We will call this the ‘scrambling velocity’ $v_{S}$. For example, in strongly scrambling models, such as random unitary circuits [33, 34, 35, 36], generic operators quickly grow into a superposition of product operators with radius $\sim v_{\mathrm{LR}}t$. In this case $v_{S}=v_{\mathrm{LR}}$.

More precisely, we define the scrambling velocity as follows. Given local operators $O_{1}$ and $O_{2}$, the commutator $[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]$ will grow along the ray ${\boldsymbol{x}}={\boldsymbol{v}}t$. We are interested in the growth of the operator itself rather than its norm or OTOC. Let $R({\boldsymbol{x}},t)$ be the radius of support of the commutator111The radius of an operator $O$ is the minimal distance $R$ such that $O$ is supported in a ball (centered at an arbitrary site) of radius $R$. Throughout the main text ‘support’ should be understood as up to an exponentially decaying tail. Exponential tails are discussed in detail in the appendices. and define

[TABLE]

This is a velocity-dependent velocity because the growth of the operator can depend on the ray that we follow, just like the exponents in (4) and (6) above. This operator growth is illustrated in Fig. 1.

In the random circuit, let $O_{1}$ and $O_{2}$ be two single-site operators. Inside the Lieb-Robinson cone, i.e. for $|{\boldsymbol{x}}|\leq v_{\mathrm{LR}}t$, the commutator $[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]$ has the same support as $O_{1}(0,t)$ so $R({\boldsymbol{x}},t)=v_{\mathrm{LR}}t$ and $v_{S}({\boldsymbol{v}})=v_{\mathrm{LR}}$ for $|{\boldsymbol{v}}|\leq v_{\mathrm{LR}}$. For general systems and for $|{\boldsymbol{v}}|\leq v_{\mathrm{LR}}$ we expect that $0\leq v_{S}({\boldsymbol{v}})\leq v_{\mathrm{LR}}$. A proof of this statement, along with more precise definitions and technical details, is collected in the appendices.

The definition (8) also captures the absence of scrambling in non-interacting theories. A non-interacting field obeys $\phi({\boldsymbol{x}},t)=\int d{\boldsymbol{y}}\,f({\boldsymbol{y}},{\boldsymbol{x}};t)\phi({\boldsymbol{y}},0)$, for some function $f({\boldsymbol{y}},{\boldsymbol{x}};t)$. Although the support of the operator $\phi({\boldsymbol{x}},t)$ spreads out as $t$ increases, it remains a superposition of local operators. Consider the conjugate pair $(\phi,\pi)$. It follows that $[\phi(0,t),\pi({\boldsymbol{x}},0)]=if({\boldsymbol{x}},0;t)$. This is a $c$-number and its support has radius $R({\boldsymbol{x}},t)=0$. Hence $v_{S}({\boldsymbol{v}})=0$ for any ${\boldsymbol{v}}$.

Even in non-interacting theories, however, more general operators — such as a pair of entangled quasiparticles moving in opposite directions — can have a nonzero scrambling velocity according to the definition (8). Relatedly, simple operators in weakly interacting theories need not have a small scrambling velocity. In this work we will mostly be interested in strongly scrambling systems. The bound we obtain will not, in general, usefully constrain weakly scrambling dynamics.

Scrambling bounds the state dependence of the OTOC

In the following subsections we prove a bound on the temperature dependence of the velocity-dependent Lyapunov exponent (6), in one spatial dimension. We also make an argument that an analogous result holds in higher dimensions. Namely:

[TABLE]

where $\beta$ is the inverse temperature, $a$ the lattice spacing, $\xi$ the correlation length, $\xi_{\mathrm{LR}}$ the microscopic lengthscale in (5), essentially the interaction range, and $h\equiv 2\sup_{{\boldsymbol{x}}\in\Lambda}\|h_{{\boldsymbol{x}}}\|$ for the Hamiltonian in (3). The content of (9) is that the change with temperature of the Lyapunov exponent along a ray is bounded by the rate of growth of the commutator along the ray. Zooming in on the butterfly light cone $v\sim v_{B}$, this bound implies that the growth of the commutator at the butterfly light cone bounds the change of characteristics such as the butterfly velocity. As (for example) the temperature is increased, these growing operators are ‘activated’ and contribute to scrambling.

A generalization, with full proof in the appendices, is as follows: For any Gibbs state $\rho=e^{-\sum_{i}\mu_{i}C^{i}}/{\rm tr}\,e^{-\sum_{i}\mu_{i}C^{i}}$ with mutually commuting conserved charges $C^{i}$, where $\mu_{i}\in{{\mathbb{R}}}$ and $C^{i}=\sum_{{\boldsymbol{x}}\in\Lambda}c^{i}_{\boldsymbol{x}}$ is a sum of local operators, then

[TABLE]

The definition of $c^{i}>0$ is similar to $h$ above: $c^{i}\equiv 2\sup_{{\boldsymbol{x}}\in\Lambda}\|c^{i}_{\boldsymbol{x}}\|$.

Outline of proof in one dimension

The following gives an outline of the proof of (9). The logic is straightforward, but technical complications arise, for example, due to the fact that time evolution generates exponentially decaying tails in space for local operators, so one cannot assume that local operators have strictly finite support. These technical points are addressed in the appendices.

Let $\rho=e^{-\beta H}/{\rm tr}\,e^{-\beta H}$ be a thermal state with inverse temperature $\beta$ and correlation length $\xi$. The steps will be as follows: (i) Differentiate the OTOC with respect to the inverse temperature, (ii) show that the main contribution to this derivative is from operators inside the support of the commutator, and (iii) balance the growth of this contribution, due to the growing size of the commutator along a ray, with the growth or decay of the OTOC. We now outline these steps.

1. (i)

Temperature derivative of the OTOC. Taking the derivative of the OTOC (1) with respect to the inverse temperature gives

[TABLE]

where $O\equiv i[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]$ and $\widetilde{H}\equiv H-{\rm tr}(\rho H)$ is the Hamiltonian with thermal expectation value subtracted out.

The Hamiltonian $H$ in (3) is written as a sum of local terms. We can split this sum up into terms that are inside and outside the support of the commutator $O$ (for some location ${\boldsymbol{x}}$ and time $t$). As in the definition of $v_{S}$, let $O$ be roughly supported in a ball of center ${\boldsymbol{y}}_{0}$ and radius $R$. Then

[TABLE]

where $\delta>0$ can take any value. As for $\widetilde{H}$, $\widetilde{h}_{{\boldsymbol{y}}}\equiv h_{\boldsymbol{y}}-{\rm tr}(\rho h_{\boldsymbol{y}})$. This decomposition can now be inserted into the derivative (11). 2. (ii)

Dominance by operators inside the commutator. We first bound the contribution from outside of the support of the commutator, with $|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|>R+\delta$ in (12). Due to the thermal correlation length $\xi$, the connected correlation function of $\widetilde{h}_{\boldsymbol{y}}$ with $O^{\dagger}O$ will decay exponentially in the distance $|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|$. Thus, for some constant $C>0$ and all ${\boldsymbol{y}}\in\Lambda$ such that $|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|>R$: $|{\rm tr}(\widetilde{h}_{\boldsymbol{y}}\sqrt{\rho}\,O^{\dagger}O\sqrt{\rho})|\leq C\|\widetilde{h}_{\boldsymbol{y}}\|\|O\|^{2}e^{(R-|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|)/\xi}$. Summing over $|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|>R+\delta$, the contribution to (11) from operators outside of the commutator is bounded by

[TABLE]

In $d$ spatial dimensions and for $R+\delta\gg\xi$, $C^{\prime}\sim C\xi(R+\delta)^{d-1}/a^{d}$ from doing the sum over $|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|>R+\delta$ ($a$ is the lattice spacing). There is a technical subtlety in obtaining (13) due to the need to commute factors of $\sqrt{\rho}$ through $\widetilde{h}_{\boldsymbol{y}}$; we deal with this in the appendices.

We can similarly bound the contribution to (11) from operators inside the support of the commutator, with $|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|\leq R+\delta$. As in the main text, define the maximal local coupling in the Hamiltonian as

[TABLE]

Note that $\|\widetilde{h}_{{\boldsymbol{y}}}\|\leq 2\|h_{{\boldsymbol{y}}}\|$, so that

[TABLE]

Notice that the inequality still goes through if we take

[TABLE]

Now, the number of terms in the first sum of (12) is $V_{R+\delta}$, the number of lattice points in a ball of radius $R+\delta$. Therefore, putting together (13) and (15), we can bound the derivative (11) by:

[TABLE]

We will see that in a certain kinematic limit, the final term in (17), from outside of the support of the commutator, is small compared to the other terms. 3. (iii)

Bounding the derivative by the growth of the commutator. The inequality (17) simplifies at late times along a ray ${\boldsymbol{x}}={\boldsymbol{v}}t$. From the definition (6) of the velocity-dependent Lyapunov exponent, $\mathcal{C}({\boldsymbol{v}}t,t;\rho)\sim e^{\lambda({\boldsymbol{v}};\rho)t}$ as $t\to\infty$. We furthermore set $\delta=(-\xi\lambda({\boldsymbol{v}};\rho)+\epsilon)t>0$, with $\epsilon>0$ a small number. This choice is such that the final term in (17) decays exponentially faster than the others as $t\to\infty$. This final term is therefore negligible in this limit. In this way, as $t\to\infty$ the following inequality is obtained:

[TABLE]

This expression bounds the temperature dependence of the Lyapunov exponent in terms of the late time growth of the commutator along a ray. The late time limit in (18) is manifestly finite in one spatial dimension, $d=1$. In one dimension at large radii $V_{r}\approx 2r/a$, where $a$ is the lattice spacing. In this case, the operator growth in (18) is precisely given by the scrambling velocity defined in (8). Thus, in terms of the scrambling velocity we obtain (A more rigorous treatment in the appendices, allowing for exponential tails in the support, shows that $\xi\to\xi+\xi_{\text{LR}}$. We include this shift in the following statement of the bound.)

[TABLE]

Generalization to higher dimensions

In higher dimensions, $V_{r}$ will scale as $r^{d}$ for $d>1$ and hence the late time bound (18) is always trivially true. However, we conjecture that the bound stated in (9) holds for arbitrary dimensions, based on a Lieb-Robinson type argument. One way of understanding the Lieb-Robinson bound is to expand

[TABLE]

and observe that in the expansion, for $[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]$ to be nonzero, a commutator sequence of local terms in $H$ connecting $O_{1}$ and $O_{2}$ is necessary, which starts at order $n\approx|{\boldsymbol{x}}|/R_{H}$ where $R_{H}$ is the range of local terms in $H$. For such a high order term to be significant, $t$ has to be later than $|{\boldsymbol{x}}|/(R_{H}h)$ and this gives an estimate of $v_{\mathrm{LR}}\approx R_{H}h$.

In a proof along these lines it is intuitively clear that outside the Lieb-Robinson cone $|{\boldsymbol{x}}|=v_{\mathrm{LR}}t$, the leading contributions to the commutator $[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]$ come from $O_{1}$ taking commutators with local terms in $H$ (as shown in (20)), via the shortest path from the origin to ${\boldsymbol{x}}$. Hence it is plausible that the operator $[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]$, for $|{\boldsymbol{x}}|\gg v_{\mathrm{LR}}t$, is approximately one-dimensional, along the line connecting [math] and ${\boldsymbol{x}}$. Then the bound (9) is still expected to be true, although possibly with a larger ‘renormalized’ $h$.

Temperature dependence of the butterfly velocity

Numerical results on the mixed field Ising chain

To motivate the general discussion of butterfly velocities, it will be useful to have some explicit numerical results for the temperature dependence of the butterfly velocity at hand. To this end we have studied the mixed field Ising chain with Hamiltonian

[TABLE]

where $X_{i}$, $Y_{i}$ and $Z_{i}$ are Pauli matrices at site $i$. Numerics is done with a straightforward generalization of the Matrix Product Operator (MPO) method discussed in [28, 19] to finite temperatures. Some analytic results on OTOCs in the transverse field model ($h_{Z}=0$) can be found in [37]. In numerics we will have $N=25$. More details can be found in the appendices. Results for the temperature dependence of the butterfly velocity for Pauli $Z$ operators are shown in Fig. 2.

The numerical results in Fig. 2 exhibit the behavior advertised in the introduction, and which we will understand in detail below. The transverse field Ising model ($h_{Z}=0$) is dual to free fermions via a Jordan-Wigner transformation. The longitudinal field $h_{Z}$ introduces interactions. We expect interactions to induce scrambling dynamics and hence a nontrivial temperature dependence of the butterfly velocity, and this is what the figure shows.

The temperature-independent butterfly velocity of the transverse field model deserves some elaboration. There are two points to make. Firstly, the transverse field model is special in its duality to a non-interacting integrable system, where $v_{S}=0$ for the commutator of fermion creation and annihilation operators, for example. For interacting integrable systems, typically $v_{S}>0$ and the butterfly velocity is state-dependent [38]. Indeed, we have verified numerically that the butterfly velocity is temperature-dependent in such models. Interacting integrable systems are scrambling, even while they are not chaotic.

Secondly, in the transverse field model, Pauli $Z$’s in the spin frame are dual to nonlocal fermion chains by the Jordan-Wigner transformation. Due to this nonlocality, our inequality doesn’t apply in the fermion frame. In fact, even local operators describing small numbers of quasiparticles in a non-interacting theory can have $v_{S}>0$ by our definition because entangled pairs of quasiparticles moving in opposite directions technically lead to a linearly growing radius of support for the operator. We believe that it may be possible to overcome this technical complication in the future with an improved definition of the scrambling velocity, such that $v_{S}=0$ for spatially separated but entangled non-scrambling operators. Indeed, we shall now argue that the butterfly velocity is temperature independent for all local operators in a non-interacting system.

In a non-interacting theory the propagation of quasiparticles is independent of the state they are propagating in, due to the absence of interactions between them. While the quasiparticles may have a nontrivial dispersion and hence temperature-dependent average velocity, any local operator includes modes of all wavevectors and, in particular, maximal velocity modes. Thus we expect $v_{B}$ is independent of the state. Therefore, the temperature-independence of the butterfly velocity observed in our numerics is indeed symptomatic of the non-interacting integrability of the system.

Bounding the butterfly velocity

The temperature dependence shown in Fig. 2 can be understood from the connections between the OTOC and scrambling velocity that we have described. The ‘light front’ form (2) for the OTOC implies that the velocity-dependent Lyapunov exponent is

[TABLE]

This precise form for $\lambda({\boldsymbol{v}};\rho)$ is conveniently explicit, but the only qualitatively essential aspect for our results is the presence of a ‘butterfly cone’. As we explained above, in general $\lambda$, $v_{B}$ and $p\geq 0$ are state-dependent. Therefore, the $\partial_{\mu_{i}}$ derivative in (10) will act on each of these quantities. Substituting the specific form (22) for $\lambda({\boldsymbol{v}};\rho)$ into (10), for $v>v_{B}$, leads to the following slightly complicated expression:

[TABLE]

where $\Delta v\equiv v/v_{B}-1>0$ is a dimensionless measure of how far the velocity is outside the butterfly cone. A simple consequence of (Bounding the butterfly velocity) follows, when there is no scrambling. Suppose that $v_{S}(v)=0$. In that case, taking $\Delta v\to 0^{+}$, the leading term on the left side of (Bounding the butterfly velocity) is the last one. It follows that

[TABLE]

Hence $v_{B}$ is constant for operators that do not scramble. We noted above, however, that this result is not directly applicable to the transverse field Ising chain.

Increasing variation of $v_{B}$ with temperature is observed in Fig. 2 as integrability is increasingly broken by turning on $h_{Z}$ in the mixed field Ising model. The crossover temperature in Fig. 2 is set by the energy gap $\Delta$ (of order $J$ for $h_{Z}=0.1\sim 0.5J$), as we now explain. Intuitively, one might expect $v_{B}$ to cease varying at temperatures $T\ll\Delta$. This is what is seen in the numerical data. We can argue for this by improving an aspect of the proof outlined previously. As we note there, the proof still goes through if we take $h$ in (9) to be instead given by

[TABLE]

where $O\equiv i[O_{1}(0,t),O_{2}({\boldsymbol{v}}t,0)]$ and $\widetilde{h}_{{\boldsymbol{y}}}\equiv h_{{\boldsymbol{y}}}-{\rm tr}(\rho h_{{\boldsymbol{y}}})$. This is not an especially tractable expression in general, but it can be evaluated for a gapped system at zero temperature, where $\rho\equiv|0\rangle\langle 0|$. In that case $h=\sup_{{\boldsymbol{y}}\in\Lambda}\langle 0|\widetilde{h}_{{\boldsymbol{y}}}|0\rangle=0$, where now $\widetilde{h}_{{\boldsymbol{y}}}\equiv h_{{\boldsymbol{y}}}-\langle 0|h_{{\boldsymbol{y}}}|0\rangle$. Hence in gapped systems at low temperatures, we may set $h\approx 0$ in the bound (9). It follows that $\partial_{\beta}v_{B}\to 0$ when $T\to 0$ in a gapped system, consistent with the finite low temperature butterfly velocities seen in Fig. 2.

The numerical results in Fig. 3 substantiate the above argument, suggesting that $\partial_{\beta}v_{B}$ decays exponentially as $\beta\Delta\to\infty$. In Fig. 3 the bound has furthermore been written as a bound on the derivative of the butterfly velocity, and is found to be most constraining at intermediate temperatures and with strong scrambling, where it is within an order of magnitude of the true value.

Our bound combined together with numerics leads to a consistent picture of the temperature dependence of the butterfly velocity in chaotic spin systems with a gap $\Delta$. Stronger scrambling allows for stronger temperature dependence of $v_{B}$, which furthermore approaches a constant at $T\ll\Delta$. These facts explain the crossover features of the curves in Fig. 2. More quantitatively, the overall variation $v_{B}(\beta=0)/v_{B}(\beta=\infty)$ can be bounded by integrating our bound from $\beta=0$ to $\beta\Delta\sim 1$ (assuming that there are no intervening thermal phase transitions). For small $v_{S}(v)$, this integration can be done explicitly, leading to a bound on the change in the butterfly velocity from infinite to zero temperature. For notational convenience let $v_{S}^{B}\equiv v_{S}(v_{B})$. At small $v_{S}^{B}$ one may take $\Delta v\sim(v_{S}^{B}/v_{B})^{1/(1+p)}$ in (Bounding the butterfly velocity) and the leading term on the left hand side is again the final one, which integrates to

[TABLE]

to leading order in $v^{B}_{S}\to 0$. Typically $v_{B}(\beta=0)\sim v_{\text{LR}}$. Schematically we can therefore write

[TABLE]

Here $\alpha$ is a dimensionful constant, $\gamma$ a dimensionless constant and we have singled out the $v_{S}$ and $\Delta$ dependences. It follows that (i) as $v_{S}\to 0$, $\ln v_{B}$ can vary as a power $v_{S}^{\gamma}$ of the scrambling velocity, and (ii) if the gap $\Delta\to 0$, $v_{B}$ may approach zero at $T=0$. Indeed, power law butterfly velocities $v_{B}\sim T^{1-1/z}$, with $z$ the dynamical critical exponent, are found in strongly chaotic gapless holographic models [16, 20].

Final comments

In summary, we have shown how locality of quantum dynamics ties operator growth to the butterfly velocity. This connection arises because the growth of the spatial support of the commutator right outside the butterfly cone bounds the change of the butterfly velocity with e.g. temperature. The butterfly velocity is state-dependent and therefore gives a richer characterization of the finite temperature dynamics than is possible from the microscopic Lieb-Robinson velocity alone. We have demonstrated these ideas explicitly in numerical studies of quantum chaotic lattice models at finite temperature. Looking forward, we hope that the methods we have developed can be used to bound other important quantities that underpin quantum many-body systems, in particular the thermalization length and time, as well as transport observables such as the thermal diffusivity.

Acknowledgements

It is a pleasure to acknowledge Jacob Marks for helping with numerics and Daren Chen for reading the proofs. We are grateful to Vedika Khemani and Xiao-Liang Qi for insightful comments on an earlier version. SAH is partially funded by DOE award de-sc0018134. XH is supported by a Stanford Graduate Fellowship. Computational work was performed on the Sherlock cluster at Stanford University, with the ITensor library for implementing tensor network calculations.

Appendices

This appendix contains six sections: section A sets up notations and backgrounds for discussions that follow. In section B we review the Lieb-Robinson, Araki and correlation length bounds used in our proof. Precise definitions for Lieb-Robinson, butterfly and scrambling velocities are given in section C and we prove several inequalities regarding them. Section D collects technical lemmas for exponentially local operators and section E gives a rigorous proof of the general results. Details of numerical implementations and data analysis are presented in section F.

Appendix A Notation

In this section we introduce notations and concepts necessary for a rigorous proof of our result. The bound will be formulated for a lattice222Technically the infinite lattice should be thought as the limit of a sequence of increasing finite subsystems. We will not delve into subtleties related to this point. $\Lambda$ of spins in $d$ spatial dimensions, and rigorously proved for $d=1$. There are isomorphic finite-dimensional Hilbert spaces ${{\cal H}}_{\boldsymbol{x}}$ associated to each lattice site ${\boldsymbol{x}}\in\Lambda$ and denote ${{\cal B}}_{\boldsymbol{x}}$ as the space of linear operators acting on ${{\cal H}}_{\boldsymbol{x}}$. An operator $O$ is said to be supported on a subset $S\subset\Lambda$ if $O\in\bigotimes_{{\boldsymbol{x}}\notin S}{{\mathbb{C}}}I\otimes\bigotimes_{{\boldsymbol{x}}\in S}{{\cal B}}_{{\boldsymbol{x}}}$, i.e. $O$ is a sum of product operators that are identity outside $S$. The minimal set that $O$ is supported on is called the support of $O$, denoted as $\operatorname{supp}O$.333Note $\operatorname{supp}O=\emptyset$ if and only if $O=cI$ for some $c\in{{\mathbb{C}}}$.

To better characterize the spatial distribution of operators, define superoperators ${{\cal P}}_{S}$ and ${{\cal Q}}_{S}\equiv\mathrm{Id}-{{\cal P}}_{S}$ such that ${{\cal P}}_{S}$ is the projection onto the subspace $\bigotimes_{{\boldsymbol{x}}\notin S}{{\mathbb{C}}}I\otimes\bigotimes_{{\boldsymbol{x}}\in S}{{\cal B}}_{{\boldsymbol{x}}}$. That is, ${{\cal P}}_{S}$ projects onto operators supported on $S$ (so ${{\cal P}}_{S}[O]=O$ if $O$ is supported on $S$). More explicitly

[TABLE]

where the integral is Haar averaging over unitaries outside $S$. However, note ${{\cal Q}}_{S}$ is not the projection onto operators supported on $\Lambda-S$. Consider an example of two sites $\Lambda=\{1,2\}$ and an operator $O=O_{1}\otimes O_{2}$, where neither $O_{1}$ nor $O_{2}$ is the identity. By definition, $0={{\cal P}}_{1}[O]={{\cal P}}_{2}[O]\neq{{\cal Q}}_{1}[O]={{\cal Q}}_{2}[O]=O$.

Henceforth if the subscript $S=\{{\boldsymbol{x}}\}$ is a single-element set, ${{\cal P}}_{\{{\boldsymbol{x}}\}}$ and ${{\cal Q}}_{\{{\boldsymbol{x}}\}}$ are written as ${{\cal P}}_{{\boldsymbol{x}}}$ and ${{\cal Q}}_{{\boldsymbol{x}}}$ for short. Also define the superoperator ${{\cal P}}_{T}^{r}$ with a superscript $r>0$ as ${{\cal P}}_{S}$ for $S=\{{\boldsymbol{y}}\in\Lambda:\exists\,{\boldsymbol{x}}\in T,|{\boldsymbol{x}}-{\boldsymbol{y}}|<r\}$, i.e. projection onto operators supported within a distance $r$ from the set $T$, and ${{\cal Q}}_{T}^{r}\equiv\mathrm{Id}-{{\cal P}}_{T}^{r}$.

From (28) we have the following inequalities:

[TABLE]

as $\|U\|=\|U^{\dagger}\|=1$. Also ${{\cal P}}_{S}[I]=I$, ${{\cal Q}}_{S}[I]=0$ for any $S\subset\Lambda$. Unless otherwise specified, $\|O\|$ will always denote the operator norm, i.e. the maximal singular value of $O$.

We will be interested primarily in operators that are “exponentially local”, denoted as ${{\cal B}}({\boldsymbol{x}},R;\xi,C)$. We say $O\in{{\cal B}}({\boldsymbol{x}},R;\xi,C)$ with ${\boldsymbol{x}}\in\Lambda$, $R,C\geq 0$ and $\xi>0$, if for any $r\geq R$,

[TABLE]

Intuitively, this means $O$ is supported on the ball of radius $R$ and centered at ${\boldsymbol{x}}$, up to an exponential tail of lengthscale $\xi$. Operators supported on a finite number of sites (called “finitely supported”) are of course exponentially local as well. We shall assume the Hamiltonian is a sum of finitely supported hermitian terms:

[TABLE]

which also defines $R_{H}>0$ and $\alpha$ labels different couplings in the Hamiltonian. Translational invariance is not necessary but $\|h^{\alpha}\|\equiv\sup_{{\boldsymbol{x}}\in\Lambda}\|h^{\alpha}_{\boldsymbol{x}}\|$ should be bounded.

A Gibbs state is a density matrix of the form

[TABLE]

for some $\mu_{i}\in{{\mathbb{R}}}$ and

[TABLE]

In the proof it is not required that $[C^{i},C^{j}]=0$. With only one $i$, with $\mu$ the inverse temperature and with $C=H$, $\rho$ is the thermal density matrix.

Appendix B Review of locality bounds

In this section we review some established locality bounds. First is the Lieb-Robinson bound in local lattice systems [27, 39, 40, 41]. This both bounds the spread of support of a local operator by the distance $v|t|$, where $t$ is the real time of Heisenberg evolution, and also implies an emergent causality with $v$ acting as the “speed of light”. For a discussion of the relation between (i) and (ii) in the following theorem, see section 3 of [42].

Theorem 1** (Lieb-Robinson).**

There exist $v,\xi_{\mathrm{LR}},C_{\mathrm{LR}}>0$, dependent on lattice geometry and Hamiltonian, such that

(i) for any $t\in{{\mathbb{R}}}$, $r>0$ and operator $O$,

[TABLE]

where $|\partial S|$ is the number of lattice links (say, between $x$ and $y$) such that $x\in S$ but $y\notin S$;

(ii) for any $t\in{{\mathbb{R}}}$, operators $O_{1}$ and $O_{2}$,

[TABLE]

where $d=\min\{|{\boldsymbol{x}}-{\boldsymbol{y}}|:{\boldsymbol{x}}\in\operatorname{supp}O_{1},{\boldsymbol{y}}\in\operatorname{supp}O_{2}\}$ is the distance between the support of $O_{1}$ and $O_{2}$.

In this bound $v\sim\sum_{\alpha}|J_{\alpha}|\|h^{\alpha}\|R_{H}$, recall (31), i.e. coupling times range of local terms in the Hamiltonian, and $\xi_{\mathrm{LR}}\sim R_{H}$. So quantities in the Lieb-Robinson bound are set by microscopic scales, to be differentiated from the butterfly velocity, which is an analog of a “renormalized” Lieb-Robinson velocity in thermal states [16].

Next is the Araki bound [42, 43, 44] extending the Lieb-Robinson bound to complex times. Note the theorem is specific to one dimension [44] and $l_{\mathrm{A}}(\mu_{i})$ may be exponential in $|\mu_{i}|$; in this sense the restriction is weaker for complex time evolution:

Theorem 2** (Araki).**

In one dimension, for any Gibbs state $\rho$ as defined in (32) but with $\mu_{i}\in{{\mathbb{C}}}$, there exist $l_{\mathrm{A}}(\mu_{i}),C_{\mathrm{A}}(\mu_{i}),\xi_{\mathrm{A}}>0$, dependent on lattice geometry and charges $C^{i}$, such that for any finitely supported operator $O$ and $r\geq l_{\mathrm{A}}(\mu_{i})$,

[TABLE]

where $|\operatorname{supp}O|$ is the number of sites in $\operatorname{supp}O$.

Note, however, from the proof of the Araki bound (e.g., Theorem 3.1 of [44]) one can see that there are Araki inequalities as stated in Theorem 2 for arbitrarily small $\xi_{\mathrm{A}}$, at the expense of a possibly large $l_{\mathrm{A}}$. Later in the proof of our bound only $\xi_{\mathrm{A}}$ enters the final expression; hence at that time one can take $\xi_{\mathrm{A}}\to 0$ as a large $l_{\mathrm{A}}$ doesn’t affect the result.

Originally the Araki bound is only stated for finitely supported operators but it is straightforward to generalize it to exponentially local ones. Such generalization will be useful in proving our bound, so a proof is given in section D.

Finally we would like to introduce some exponential clustering theorems: for particular kinds of states, equal-time connected correlations decay exponentially in space. More precisely for a state (density matrix) $\rho$, the correlation length of $\rho$ is the $\xi>0$ that is optimal with respect to the following property: there exists $C>0$ and a function $l_{0}(\cdot)>0$ such that for any operators $O_{1}$ and $O_{2}$ (supported on sets $S$ and $T$) sufficiently far apart, i.e., $d\geq l_{0}(\delta)$,

[TABLE]

where $\delta\equiv\min\{|\partial S|,|\partial T|\}$ is the number of lattice links crossing the boundary of $S$ or $T$, and $d\equiv\min\{|{\boldsymbol{x}}-{\boldsymbol{y}}|:{\boldsymbol{x}}\in S,{\boldsymbol{y}}\in T\}$ is the distance between two sets. Note that in this statement, $O_{1}$ and $O_{2}$ could be any, not necessarily local, operators.

Existence of a finite $\xi>0$ with the property stated around (38) has been proved for (i) one-dimensional Gibbs states [43] (restricted to local operators $O_{1}$ and $O_{2}$), (ii) $\rho=|0\rangle\langle 0|$ where $|0\rangle$ is the unique ground state of a gapped Hamiltonian [40, 45], and (iii) thermal states $\rho\propto\exp(-\beta H)$ in general dimensions at sufficiently high temperatures [46] (clearly $\xi\to 0$ when $\beta\to 0$). Of course the Hamiltonians associated with these states must be local, as in (31) above. It is plausible that the correlation length $\xi$ as defined around (38) is finite for Gibbs states $\rho$ in general systems with local dynamics and away from phase transitions.

Appendix C Definitions of velocities

In this section we define precisely the (possibly anisotropic) Lieb-Robinson, butterfly and scrambling velocities introduced in the main text and prove the bound $v_{B},v_{S}\leq v_{\mathrm{LR}}$. For definiteness fix a class of local operators, denoted as $\mathcal{O}$; for example, $\mathcal{O}$ could be all single-site operators with unit norm, localized at origin. The Lieb-Robinson bound Theorem 1 (ii) can be stated for such operators along any particular direction $\boldsymbol{n}$:

Theorem 3** (Operator-dependent anisotropic Lieb-Robinson).**

For any direction $\boldsymbol{n}$ and operator $O_{1},O_{2}\in\mathcal{O}$, there exist $v$, $\xi_{\mathrm{LR}}$, $C_{\mathrm{LR}}>0$, dependent on $\boldsymbol{n}$, $O_{1}$, $O_{2}$, lattice geometry and Hamiltonian, such that for any $t>0$, $x>0$,

[TABLE]

From Theorem 3 one immediate candidate for defining the Lieb-Robinson velocity is

[TABLE]

that is, the smallest velocity with a Lieb-Robinson inequality. However such a definition shows some disadvantages in numerical or experimental applications: it is inaccurate to fit data to exponential tails because the theorem only states an inequality (not an equality), and in fact in many lattice systems of interest the tail is observed to be sub-exponential (e.g., Gaussian) [18, 19]; also it is impractical, if not impossible, to decide whether such $\xi_{\mathrm{LR}}$ and $C_{\mathrm{LR}}$ exist for all times, from only a finite number of data points.

A more convenient definition is found in the original Lieb-Robinson paper [27]

[TABLE]

We will assume that the limit exists and is a continuous function of $v$. By definition $v_{\mathrm{LR}}^{(2)}$ gives a causality “lightcone” outside which (for $x/t>v$) the commutator vanishes exponentially at late times.

It is relatively easy to see that $v^{(1)}_{\mathrm{LR}}\geq v^{(2)}_{\mathrm{LR}}$:

Proposition 1**.**

For any direction $\boldsymbol{n}$ and operators $O_{1},O_{2}\in\mathcal{O}$, we have $v^{(1)}_{\mathrm{LR}}(\boldsymbol{n};O_{1},O_{2})\geq v^{(2)}_{\mathrm{LR}}(\boldsymbol{n};O_{1},O_{2})$.

Proof.

Let $v>0$ belong to the set in (40), i.e., there exist $\xi,C>0$ such that for all $x,t>0$, $\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|\leq C\|O_{1}\|\|O_{2}\|\min\{1,e^{(vt-x)/\xi}\}$. Then, for any $v^{\prime}>v$, $\lim_{t\to\infty}t^{-1}\ln\|[O_{1}(0,t),O_{2}(v^{\prime}t\boldsymbol{n},0)]\|\leq\lim_{t\to\infty}t^{-1}\ln(C\|O_{1}\|\|O_{2}\|e^{(v-v^{\prime})t/\xi})=(v-v^{\prime})/\xi<0$, and hence any $v^{\prime}>v$ is not contained in the set in (41). Therefore the supremum $v^{(2)}_{\mathrm{LR}}$ is at most $v$. This is true for any $v>0$ in the set of (40), hence $v^{(2)}_{\mathrm{LR}}\leq v^{(1)}_{\mathrm{LR}}$. ∎

Conversely to show that $v_{\mathrm{LR}}^{(1)}\leq v_{\mathrm{LR}}^{(2)}$, we need the following lemma:

Lemma 1**.**

For any positive functions $f(x,t)$ and $g(x,t)$, if limits

[TABLE]

exist, are uniform for $v\in[v_{0},\infty)$, and $\lambda_{f}(v)+a<\lambda_{g}(v)$ for some $a>0$ and all $v\geq v_{0}$, then there is $t_{0}>0$ that

[TABLE]

Proof.

Because the limits (42) are uniform, for any $\varepsilon>0$ there is $T(\varepsilon)>0$ such that for any $t\geq T(\varepsilon)$ and $v\geq v_{0}$, $\ln f(vt,t)/t<\lambda_{f}(v)+\varepsilon$, $\ln g(vt,t)/t>\lambda_{g}(v)-\varepsilon$. Now choose $\varepsilon=a/2$ and $t_{0}=T(a/2)$, we have $\ln f(vt,t)/t<\lambda_{f}(v)+a/2<\lambda_{g}(v)-a/2<\ln g(vt,t)/t$ hence $f(vt,t)<g(vt,t)$, for all $t\geq t_{0}$, $v\geq v_{0}$. ∎

Proposition 2**.**

$v^{(1)}_{\mathrm{LR}}(\boldsymbol{n};O_{1},O_{2})\leq v^{(2)}_{\mathrm{LR}}(\boldsymbol{n};O_{1},O_{2})$, given the limit in (41) is uniform for all $v>v^{(2)}_{\mathrm{LR}}(\boldsymbol{n};O_{1},O_{2})$.

Proof.

We would like to prove the proposition in the following two steps:

Step one: For any $v>v^{(2)}_{\mathrm{LR}}$, we show that (i) implies (ii), and (ii) implies (iii), where

1. (i)

$\lim_{t\to\infty}t^{-1}\ln\|[O_{1}(0,t),O_{2}(v^{\prime}t\boldsymbol{n},0)]\|<0$ for any $v^{\prime}\geq v$; 2. (ii)

$\exists\,\varepsilon,\xi>0$ that $\lim_{t\to\infty}t^{-1}\ln\|[O_{1}(0,t),O_{2}(v^{\prime}t\boldsymbol{n},0)]\|\leq(v-v^{\prime})/\xi-\varepsilon$ for any $v^{\prime}\geq v$; 3. (iii)

$\exists\,C,\xi>0$ that $\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|\leq C\|O_{1}\|\|O_{2}\|\min\{1,e^{(vt-x)/\xi}\}$ for $x,t>0$.

Step two: By definition (41) we have for any $v>v^{(2)}_{\mathrm{LR}}$, (i) holds for $v$; so (iii) is true for $v$ as well, and $v$ should be in the set on the right-hand side of (40) hence $v^{(1)}_{\mathrm{LR}}\leq v$. This shows that $v^{(1)}_{\mathrm{LR}}\leq v^{(2)}_{\mathrm{LR}}$.

So now it remains to prove that (i) $\Rightarrow$ (ii) and (ii) $\Rightarrow$ (iii):

(i) $\Rightarrow$ (ii): For clarity let’s denote $\lambda(v)\equiv\lim_{t\to\infty}t^{-1}\ln\|[O_{1}(0,t),O_{2}(vt\boldsymbol{n},0)]\|$, then (i) says that $\lambda(v^{\prime})<0$ for any $v^{\prime}\geq v$ and to arrive at (ii) we hope to find $\varepsilon,\xi>0$ such that $\lambda(v^{\prime})\leq(v-v^{\prime})/\xi-\varepsilon$ for all $v^{\prime}\geq v$.

Before construction of $\varepsilon$ and $\xi$, it is remarkable that there is a restriction on $\lambda(v^{\prime})$ from Theorem 3: the Lieb-Robinson bound states that there are some $C_{0}$, $v_{0}$, $\xi_{0}>0$ such that $\lambda(v^{\prime})\leq\lim_{t\to\infty}t^{-1}\ln(C_{0}\|O_{1}\|\|O_{2}\|e^{(v_{0}-v^{\prime})t/\xi_{0}})=(v_{0}-v^{\prime})/\xi_{0}$ for all $v^{\prime}>0$.

We shall construct $\varepsilon>0$ first. Note $(v-v^{\prime})/\xi\leq 0$ for $v^{\prime}\geq v$, hence it is required that $\lambda(v^{\prime})\leq-\varepsilon$ for all $v^{\prime}\geq v$. So we may choose $\varepsilon=\inf_{v^{\prime}\geq v}(-\lambda(v^{\prime})/2)\geq 0$. To show that $\varepsilon>0$, we have to check that $-\lambda(v^{\prime})>0$ is bounded from zero on $[v,\infty)$. The only concern is $\lambda(v^{\prime})$ may be arbitrarily close to zero when $v^{\prime}\to\infty$; but this is not possible because from the previous paragraph $-\lambda(v^{\prime})\geq(v^{\prime}-v_{0})/\xi_{0}\to\infty$ as $v^{\prime}\to\infty$. Hence $\varepsilon>0$ is well-defined in this way.

Then to satisfy $\lambda(v^{\prime})\leq(v-v^{\prime})/\xi-\varepsilon$ for all $v^{\prime}\geq v$, choose ($\xi_{0}$ is there for future convenience) $\xi\equiv\max\{\xi_{0},\sup_{v^{\prime}\geq v}(v-v^{\prime})/(\lambda(v^{\prime})+\varepsilon)\}$ (as constructed in the last paragraph the denominator is always negative). The task is then to show that $\xi<\infty$; similarly the only place things could go wrong is when $v^{\prime}\to\infty$, but in that limit $|\lambda(v^{\prime})+\varepsilon|\geq|\lambda(v^{\prime})|/2\geq(v^{\prime}-v_{0})/2\xi_{0}$ hence $\lim_{v^{\prime}\to\infty}(v-v^{\prime})/(\lambda(v^{\prime})+\varepsilon)\leq 2\xi_{0}$ is bounded. So $\xi>0$ is well-defined as well and (ii) is proved.

(ii) $\Rightarrow$ (iii): We would like to apply the Lemma 1 for $f(x,t)=\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|$ and $g(x,t)=\|O_{1}\|\|O_{2}\|e^{(vt-x)/\xi}$. Note in this case $\lambda_{f}(v^{\prime})=\lambda(v^{\prime})\leq(v-v^{\prime})/\xi-\varepsilon=\lambda_{g}(v^{\prime})-\varepsilon$ for any $v^{\prime}\geq v$. Then by the lemma there is $t_{0}>0$ such that $\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|\leq\|O_{1}\|\|O_{2}\|e^{(vt-x)/\xi}$ for all $x\geq vt$ and $t\geq t_{0}$. Hence for (iii) to hold it suffices to choose that $C\equiv\max\{2,\sup_{0<x<vt\text{ or }0<t<t_{0}}f(x,t)/g(x,t)\}$. As before we have to check that the supremum is not infinite. We will discuss the three cases (a) $0<x<vt$, (b) $0<t<t_{0}$ with $x\geq v_{0}t$, and (c) $0<t<t_{0}$ with $0<x<v_{0}t$ separately.

For $0<x<vt$, $f(x,t)/g(x,t)=\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|/\|O_{1}\|\|O_{2}\|e^{(vt-x)/\xi}$ is less than $\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|/\|O_{1}\|\|O_{2}\|\leq 2$. So indeed $f(x,t)/g(x,t)$ is bounded in this region.

For $0<t<t_{0}$ and $x\geq v_{0}t$, $f(x,t)/g(x,t)=\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|/\|O_{1}\|\|O_{2}\|e^{(vt-x)/\xi}$ can be bounded using the Lieb-Robinson Theorem 3: there is some $C_{0},v_{0},\xi_{0}>0$ such that $\|[O_{1}(0,t),O_{2}(x\boldsymbol{n},0)]\|\leq C_{0}\|O_{1}\|\|O_{2}\|e^{(v_{0}t-x)/\xi_{0}}\leq C_{0}\|O_{1}\|\|O_{2}\|e^{(v_{0}t-x)/\xi}$ (by construction $\xi\geq\xi_{0}$) so $f(x,t)/g(x,t)\leq C_{0}e^{(v_{0}-v)t/\xi}$ which is a bounded function for $0<t<t_{0}$.

Finally for $0<t<t_{0}$ and $0<x<v_{0}t$, $f(x,t)/g(x,t)$ is bounded because it is continuous and the region is bounded. Hence we’ve shown that $C>0$ is well-defined and with $\xi$ appearing in (ii), (iii) is true. ∎

Henceforth the Lieb-Robinson velocity will be defined as $v_{\mathrm{LR}}\equiv v_{\mathrm{LR}}^{(1)}=v_{\mathrm{LR}}^{(2)}$. The technical uniformity condition is true for known examples. The same proof shows the equivalence of two definitions of the butterfly velocity. For future use only the definition corresponding to $v_{\mathrm{LR}}^{(2)}$ is recorded:

[TABLE]

where the OTOC $\mathcal{C}_{O_{1}O_{2}}({\boldsymbol{x}},t;\rho)$ is defined in (1). As the velocity-dependent quantum Lyapunov exponent is defined as in (6), an equivalent definition of $v_{B}$ reads:

[TABLE]

As expected, the butterfly velocity in any state is bounded by the Lieb-Robinson velocity:

Proposition 3**.**

$v_{B}(\boldsymbol{n};O_{1},O_{2},\rho)\leq v_{\mathrm{LR}}(\boldsymbol{n};O_{1},O_{2})$* for any $O_{1},O_{2}\in\mathcal{O}$, density matrix $\rho$ and direction $\boldsymbol{n}$.*

Proof.

This follows from definition (41) and (44), and $\mathcal{C}_{O_{1}O_{2}}({\boldsymbol{x}},t;\rho)\leq\|[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]\|^{2}$. ∎

Finally the scrambling velocity can be precisely defined in the language of exponentially local operators, defined around (30). Let $O\equiv i[O_{1}(0,t),O_{2}({\boldsymbol{v}}t,0)]$, then444To make sure the limit exists, we have used the limit superior $\operatorname*{\overline{lim}}$ and the limit inferior $\operatorname*{\underline{lim}}$.

[TABLE]

where the smallest ball, with radius $R$ and centered at ${\boldsymbol{x}}$, is understood as roughly the “support” of the commutator $O$. The quantities $\xi$ and $C$ characterize the exponential tail that we neglected in the main text. Clearly $v_{S}\geq 0$ and decreases with increasing $\xi$.

For any triple $(v,\xi\equiv\xi_{\mathrm{LR}},C_{\mathrm{LR}})$ from Theorem 1, we now show that $v_{S}({\boldsymbol{v}};\xi)\leq v$. Thus we have an upper bound of $v_{S}$ by velocities with a Lieb-Robinson inequality. Note the $\xi$-dependence of $v_{S}$ was omitted in the main text. More precisely, if $O_{2}({\boldsymbol{v}}t,0)$ is within the “support” of $O_{1}(0,t)$, for scrambling systems at late times we would expect $\|[O_{1}(0,t),O_{2}({\boldsymbol{v}}t,0)]\|$ to equilibrate to a nonzero constant value; if so, $v_{S}\leq v$:

Proposition 4**.**

Given ${\boldsymbol{v}}$, $\xi>0$, $O_{1},O_{2}\in\mathcal{O}$, if for any $t>0$, $O_{1}(0,t)\in{{\cal B}}(0,vt;\xi,C)$ for some $v>|{\boldsymbol{v}}|$, $C>0$ and $\operatorname*{\underline{lim}}_{t\to\infty}\|[O_{1}(0,t),O_{2}({\boldsymbol{v}}t,0)]\|>0,$ then $v_{S}({\boldsymbol{v}};O_{1},O_{2},\xi)\leq v$.

Proof.

Let $O(t)\equiv[O_{1}(0,t),O_{2}({\boldsymbol{v}}t,0)]$, $c\equiv\operatorname*{\underline{lim}}_{t\to\infty}\|O(t)\|>0$. As $|{\boldsymbol{v}}|<v$, ${{\cal Q}}_{0}^{r}[O_{2}({\boldsymbol{v}}t,0)]=0$ for $r\geq vt$ at late times. Then $O(t)=[{{\cal P}}_{0}^{r}[O_{1}(0,t)],{{\cal P}}_{0}^{r}[O_{2}({\boldsymbol{v}}t,0)]]+[{{\cal Q}}_{0}^{r}[O_{1}(0,t)],{{\cal P}}_{0}^{r}[O_{2}({\boldsymbol{v}}t,0)]]$. But the first term is supported in the ball of radius $r$ centered at origin, so $\|{{\cal Q}}_{0}^{r}[O(t)]\|=\|{{\cal Q}}_{0}^{r}[{{\cal Q}}_{0}^{r}[O_{1}(0,t)],{{\cal P}}_{0}^{r}[O_{2}({\boldsymbol{v}}t,0)]]\|\leq 4\|{{\cal Q}}_{0}^{r}[O_{1}(0,t)]\|\|{{\cal P}}_{0}^{r}[O_{2}({\boldsymbol{v}}t,0)]\|\leq 4C\|O_{1}\|\|O_{2}\|e^{(vt-r)/\xi}$, where we have used the definition (30) that for all $t>0$ and $r\geq vt$, $\|{{\cal Q}}_{0}^{r}[O_{1}(0,t)]\|\leq C\|O_{1}\|e^{(vt-r)/\xi}$ with the inequalities (29).

So there is a time $t_{0}>0$ that for all $t>t_{0}$, $\|O(t)\|\geq c/2$ as well as $\|{{\cal Q}}_{0}^{r}[O(t)]\|\leq 4C\|O_{1}\|\|O_{2}\|e^{(vt-r)/\xi}$ for all $r\geq vt$. Hence $\|{{\cal Q}}_{0}^{r}[O(t)]\|\leq C^{\prime}\|O(t)\|e^{(vt-r)/\xi}$, for all $t>t_{0}$ and $r\geq vt$, if we choose $C^{\prime}=8C\|O_{1}\|\|O_{2}\|/c$. That is, $O(t)\in{{\cal B}}(0,vt;\xi,C^{\prime})$ for $t>t_{0}$ hence by definition (46), $v_{S}({\boldsymbol{v}};O_{1},O_{2},\xi)\leq v$. ∎

All velocities can be maximized over direction $\boldsymbol{n}$ to recover their isotropic definitions, or over $O_{1},O_{2}\in\mathcal{O}$ to remove the operator dependence.

Appendix D Bounds for exponentially local operators

In this section we collect some lemmas and generalize Theorem 2 and the exponential clustering condition (38) to exponentially local operators. Readers are encouraged to review sections A and B. The following inequality will be useful: for any $A,B\geq 0$ and $k,\gamma>0$,

[TABLE]

where $\lceil x\rceil$ denotes the least integer greater than or equal to $x$. To show this, by doing the summation exactly it is easy to check that for any $A,B\geq 0$, $\gamma>0$ and integer $m\geq 1$,

[TABLE]

and the inequality (47) follows because if $m=\lceil k\rceil$, $m\leq k+1$ in the linear factor and $k\leq m$ implies that $e^{-\gamma m}\leq e^{-\gamma k}$ as well.

The following lemma bounds the product of two exponentially local operators:

Lemma 2**.**

Let $O_{1}\in{{\cal B}}({\boldsymbol{x}},R;\xi_{1},C_{1})$ and $O_{2}\in{{\cal B}}({\boldsymbol{x}},R;\xi_{2},C_{2})$, then for any $r\geq R$,

[TABLE]

Proof.

Note that for any $r>0$, $O_{1}O_{2}={{\cal P}}_{\boldsymbol{x}}^{r}[O_{1}]{{\cal P}}_{\boldsymbol{x}}^{r}[O_{2}]+O_{1}{{\cal Q}}_{\boldsymbol{x}}^{r}[O_{2}]+{{\cal Q}}_{\boldsymbol{x}}^{r}[O_{1}]{{\cal P}}_{\boldsymbol{x}}^{r}[O_{2}]$, and ${{\cal Q}}_{\boldsymbol{x}}^{r}[{{\cal P}}_{\boldsymbol{x}}^{r}[O_{1}]{{\cal P}}_{\boldsymbol{x}}^{r}[O_{2}]]=0$. So by (29) and (30), for $r\geq R$,

[TABLE]

∎

Next is the Araki bound (cf. Theorem 2) for exponentially local operators:

Theorem 4**.**

For any one-dimensional Gibbs state $\rho$ as defined in (32) with $\mu_{i}\in{{\mathbb{C}}}$ and operator $O\in{{\cal B}}({\boldsymbol{x}},R;\xi,C)$, there exists $C^{\prime}(\mu_{i},\xi,C)>0$ (dependent on lattice geometry and $C^{i}$ as well) such that for all $r\geq R+l_{\mathrm{A}}(\mu_{i})+a$,

[TABLE]

Here $l_{\mathrm{A}}(\mu_{i})$ and $\xi_{\mathrm{A}}$ are those appearing in the Araki bound, and $a$ is the lattice spacing.

Proof.

For the first inequality, let $m\equiv\lceil(R+a)/a\rceil$. Decompose $O={{\cal P}}_{\boldsymbol{x}}^{(m-1)a}[O]+\sum_{n\geq m}O_{n},$ where $O_{n}\equiv{{\cal P}}_{\boldsymbol{x}}^{na}{{\cal Q}}_{\boldsymbol{x}}^{(n-1)a}[O]={{\cal P}}_{{\boldsymbol{x}}}^{na}[O]-{{\cal P}}_{{\boldsymbol{x}}}^{(n-1)a}[O]$. Then by Theorem 2 with (29) and (30), for $n\geq m$,

[TABLE]

Also by Theorem 2, $\|{{\cal P}}[O]\|\leq\|O\|$ and $m\leq(R+a)/a+1=R/a+2$,

[TABLE]

Sum (D) with (47) (where $A=2$, $B=1$, $k=(R+a)/a$ and $\gamma=a/\xi$) to get the bound

[TABLE]

Denote $C_{1}(\mu_{i},\xi,C)=3C_{\mathrm{A}}(\mu_{i})+5C_{\mathrm{A}}(\mu_{i})C(1-e^{-a/\xi})^{-2}$, so that

[TABLE]

For the second inequality, expand $O={{\cal P}}_{\emptyset}[O]+\sum_{n\geq 0}O_{n},$ where ${{\cal P}}_{\emptyset}[O]$ is proportional to identity and $O_{n}\equiv{{\cal P}}_{\boldsymbol{x}}^{na}{{\cal Q}}_{\boldsymbol{x}}^{(n-1)a}[O]={{\cal P}}_{{\boldsymbol{x}}}^{na}[O]-{{\cal P}}_{{\boldsymbol{x}}}^{(n-1)a}[O]$. Because ${{\cal Q}}_{\boldsymbol{x}}^{r}[I]=0$,

[TABLE]

Let $\delta\equiv\alpha(r-l_{\mathrm{A}}(\mu_{i})-R-a)\geq 0$ for any $0<\alpha<1$ and split the sum (57) into two parts: $0\leq na<R+\delta+a$ and $na\geq R+\delta+a$. Apply Theorem 2 for the first part (also note $\|O_{n}\|\leq 2\|O\|$ by (29)):

[TABLE]

and further with inequalities (29) and definition (30) for the second part:

[TABLE]

Overall, sum (58) as geometric series after applying $n\leq k$ and sum (59) with (47) (where $A=2$, $B=1$, $k=(R+\delta+a)/a$ and $\gamma=a/\xi$):

[TABLE]

where in the second inequality we have replaced $ka=R+\delta+a$ in the exponents and applied the bound $k\leq(r-l_{\mathrm{A}}(\mu_{i}))/a$ (because $\alpha\leq 1$) in the prefactors. Now

[TABLE]

if one chooses $\alpha=\xi/(\xi_{\mathrm{A}}+\xi)$ to equate the exponents and $C_{2}(\mu_{i},\xi,C)=2C_{\mathrm{A}}(\mu_{i})(1-e^{-a/\xi_{\mathrm{A}}})^{-1}+6CC_{\mathrm{A}}(\mu_{i})(1-e^{-a/\xi})^{-2}$.

Finally it suffices to choose $C^{\prime}(\mu_{i},\xi,C)\equiv\max\{C_{1}(\mu_{i},\xi,C),C_{2}(\mu_{i},\xi,C)\}$. ∎

Observe that the operator $\rho O\rho^{-1}$ as stated in (52), is not exponentially local explicitly (due to the prefactor that is linear in $r$). To work around this the following corollary of Theorem 4 is particularly useful:

Corollary 1**.**

For any $\varepsilon>0$, there is a $\widetilde{C}^{\prime}(\mu_{i},\xi,C,\varepsilon)$ such that

[TABLE]

Proof.

First note that for $\zeta(\xi)\equiv\xi_{\mathrm{A}}+\xi$,

[TABLE]

so it suffices to find $\widetilde{C}^{\prime}(\mu_{i},\xi,C,\varepsilon)$ such that for all $x\equiv r-l_{\mathrm{A}}(\mu_{i})-a\geq 0$,

[TABLE]

which clearly exists. ∎

Finally we generalize inequality (38) to exponentially local operators as well; for future use we will work in one dimension only:

Theorem 5**.**

Let $\rho$ be a one-dimensional state with $\xi$, $C$ and $l_{0}(\cdot)>0$ as stated around (38). If $O_{1}\in{{\cal B}}({\boldsymbol{x}},R_{1};\xi_{1},C_{1})$, $O_{2}\in{{\cal B}}({\boldsymbol{y}},R_{2};\xi_{2},C_{2})$ and $|{\boldsymbol{x}}-{\boldsymbol{y}}|\geq l_{0}(2)+R_{1}+R_{2}$,

[TABLE]

Proof.

Let $\Delta\equiv|{\boldsymbol{x}}-{\boldsymbol{y}}|-l_{0}(2)-R_{1}-R_{2}\geq 0$, and define $r\equiv R_{1}+\alpha_{1}\Delta$ and $s\equiv R_{2}+\alpha_{2}\Delta$ for $\alpha_{1},\alpha_{2}>0$ and $\alpha_{1}+\alpha_{2}<1$. Denote $c(O_{1},O_{2})\equiv{\rm tr}(\rho\,O_{1}O_{2})-{\rm tr}(\rho\,O_{1}){\rm tr}(\rho\,O_{2})$ for convenience and observe $|c(O_{1},O_{2})|\leq 2\|O_{1}\|\|O_{2}\|$. Then

[TABLE]

By inequality (38), (note $\delta=2$ if $S$ and $T$ are intervals in (38) and $\|{{\cal P}}[O]\|\leq\|O\|$)

[TABLE]

and by definition (30),

[TABLE]

Now choose $\alpha_{1}=\xi_{1}/(\xi+\xi_{1}+\xi_{2})$ and $\alpha_{2}=\xi_{2}/(\xi+\xi_{1}+\xi_{2})$ so that the exponents with $\Delta$ are all equal. Sum them up to get (65). ∎

Appendix E Proof of the bound

In this section we give a proof of the bounds stated in the main text. To avoid clutter of notations, all quantities in this section may depend on lattice geometry, Hamiltonian $H$ (31) and charges $C^{i}$ (33) implicitly.

Theorem 6**.**

For any one-dimensional Gibbs state $\rho$ as defined in (32) with correlation length $\xi_{\mathrm{cor}}$ (read around (38) for a definition), $\varepsilon,\delta>0$, any operators $O_{1},O_{2}$ and ${\boldsymbol{x}}\in\Lambda$, $t>0$, there exist $A(\mu_{i},\xi,C,\varepsilon)$, $B(\mu_{i})>0$ such that

[TABLE]

and

[TABLE]

where $a$ is the lattice spacing and $\xi_{\mathrm{A}}$ is defined in Theorem 2. The inverse temperature is denoted as $\beta$ and $J_{\alpha}$ labels couplings in the Hamiltonian (31). Denote $O\equiv i[O_{1}(0,t),O_{2}({\boldsymbol{x}},0)]$; $R$, $\xi$ and $C$ are such that $O\in{{\cal B}}({\boldsymbol{y}}_{0},R;\xi,C)$ for some ${\boldsymbol{y}}_{0}\in\Lambda$. Finally

[TABLE]

where $\widetilde{c}^{i}_{\boldsymbol{y}}\equiv c^{i}_{\boldsymbol{y}}-{\rm tr}(\rho\,c^{i}_{\boldsymbol{y}})$, and same for $h^{\alpha}$ with $c^{i}_{\boldsymbol{y}}$ replaced by $h^{\alpha}_{{\boldsymbol{y}}}$. And if $C^{i}$ commute with each other, $c^{i}$ can be chosen as

[TABLE]

Proof.

We start with proving (70). By definition (1) and (32),

[TABLE]

where for any operator $C$, $\widetilde{C}\equiv C-{\rm tr}(\rho\,C)$. Now recall $C^{i}$ is a sum of local terms (33):

[TABLE]

for any $S(r)\equiv\{{\boldsymbol{y}}\in\Lambda:|{\boldsymbol{y}}-{\boldsymbol{y}}_{0}|\leq r\}$. For any ${\boldsymbol{y}}\in\Lambda$, by definition of $c^{i}(s)$,

[TABLE]

The inequality (76) is good enough for the terms in $S(r)$. For the remaining terms with ${\boldsymbol{y}}$ away from ${\boldsymbol{y}}_{0}$ we have a better estimate because connected correlation decays when operators are far apart. There is a technical complication due to the fact that the factors of $\rho$ are separated by – and do not necessarily commute with – the $\widetilde{c}^{i}_{{\boldsymbol{y}}}$. For this reason we need to use the Araki bound to show that operators remain sufficiently local under conjugation by the density matrix. Indeed by Lemma 2, $O^{\dagger}O\in{{\cal B}}({\boldsymbol{y}}_{0},R;\xi,4C)$ and from Theorem 4 and Corollary 1, there is $C_{1}(\mu_{i},\xi,C,\varepsilon)>0$ and $l(\mu_{i})>0$ such that for any $0\leq s\leq 1$,

[TABLE]

[TABLE]

Hence by Theorem 5, because ${\rm tr}(\rho\,\widetilde{c}^{i}_{\boldsymbol{y}})=0$, for any $0\leq s\leq 1$,

[TABLE]

where $C_{2}$ is defined in terms of the prefactor $C_{\mathrm{cor}}(\mu_{i})$ in (38) as, using (77),

[TABLE]

Now bound the sum (75) by choosing $r=R+l(\mu_{i})+a+R_{H}+l_{0}(2)+\delta$ and apply (76) for ${\boldsymbol{y}}\in S(r)$ and (E) for ${\boldsymbol{y}}\notin S(r)$, (denote $\zeta\equiv\xi_{\mathrm{cor}}+\xi_{\mathrm{A}}+\xi+\varepsilon$)

[TABLE]

and use the inequality (80) and $\|\widetilde{c}^{i}_{{\boldsymbol{y}}}\|\leq 2\|c^{i}_{{\boldsymbol{y}}}\|$ to reduce to the form (70).

Proving (71) is essentially the same except in the first step:

[TABLE]

there is an additional term due to coupling dependence of $O_{1}(0,t)$. By definition,

[TABLE]

and (71) follows from the Cauchy-Schwartz inequality for the inner product $\langle O_{1},O_{2}\rangle\equiv{\rm tr}(\rho\,O_{1}^{\dagger}O_{2})$.

Finally if $C^{i}$ commute with each other, the first step (74) can be replaced with

[TABLE]

and the same proof goes through with $c^{i}$ as in (73). It is bounded by $2\sup\|c^{i}_{\boldsymbol{y}}\|$ because $\sqrt{\rho}\,O^{\dagger}O\sqrt{\rho}$ is a positive operator and for any operator $S$ and positive operator $T$, $|{\rm tr}ST|\leq\|S\|{\rm tr}\,T$. ∎

The theorem, as stated, seems complicated; but the physics is much clearer in terms of the velocity-dependent Lyapunov exponent (6):

Corollary 2**.**

For $v_{S}({\boldsymbol{v}};O_{1},O_{2},\xi)$ defined in (46),

[TABLE]

Proof.

Divide both sides of (70) by $t\,\mathcal{C}_{O_{1}O_{2}}({\boldsymbol{x}},t;\rho)$, choose

[TABLE]

${\boldsymbol{x}}={\boldsymbol{v}}t$ and take the limit $t\to\infty$ (assuming the limit and derivative commute):

[TABLE]

Finally let $\varepsilon,\xi_{\mathrm{A}}\to 0$ to conclude555Regarding the limit $\xi_{\mathrm{A}}\to 0$ we refer readers to the discussions following Theorem 2.. ∎

The operator $O$ must decay at large distances at least as quickly as the rate set by $\xi_{\mathrm{LR}}$ (appearing in any triple $(v,\xi_{\mathrm{LR}},C_{\mathrm{LR}})$ with a Lieb-Robinson bound Theorem 1). Therefore we take $\xi=\xi_{\text{LR}}$ in the main text. We have already noted in section C that this then defines a $v_{S}({\boldsymbol{v}};\xi)\leq v$.

The coupling dependence of $\lambda_{O_{1}O_{2}}({\boldsymbol{v}};\rho)$ can be bounded in the same way:

Corollary 3**.**

If $\mathcal{C}_{O_{1}O_{2}}({\boldsymbol{v}}t,t;\rho)\sim\kappa_{1}^{2}e^{\lambda_{O_{1}O_{2}}({\boldsymbol{v}};\rho)t}$ and

[TABLE]

for $\kappa_{1},\kappa_{2}>0$ and $\|h^{\alpha}\|\equiv\sup_{{\boldsymbol{y}}\in\Lambda}\|h^{\alpha}_{\boldsymbol{y}}\|$ at $t\to\infty$,

[TABLE]

If we assume that the growth rate of the OTOC does not depend on choices of operators, i.e., the growth rate in (88) is $\lambda_{O_{1}O_{2}}({\boldsymbol{v}};\rho)$, the same as that of $\mathcal{C}_{O_{1}O_{2}}({\boldsymbol{v}}t,t;\rho)$, this corollary shows that divergence of $\partial_{J}v_{B}$ at zero temperature pinpoints quantum phase transitions at which the system becomes gapless. Indeed, if to the contrary the system is gapped, as observed in Fig. 3 and discussed in the main text, the first term on the right side of (89) is expected to vanish at zero temperature so the right-hand side of (89) should be finite, contradicting the divergence of $\partial_{J}v_{B}$ via an inequality similar to (Bounding the butterfly velocity). Cusps of scrambling characteristics are indeed observed at quantum critical points in e.g. [47, 48].

Appendix F Numerical details

Our method is a generalization of the Matrix Product Operator (MPO) approach to calculating the butterfly velocity, presented in [19], to finite temperature states. The algorithm is implemented with the ITensor library, with operators $O_{1}(0,t)$, $O_{2}({\boldsymbol{x}},0)$ and thermal density matrix $\rho$ represented as MPOs and evolved with a Time-Evolving Block Decimation (TEBD) method (for MPOs). For general quantum systems the thermal entanglement entropy is expected to be extensive. We find in practice that the MPO representation of thermal states works at sufficiently high but finite temperatures (in our case, $0\leq\beta J\leq 3$). Numerical truncation $\varepsilon$ in the MPO is set to $\varepsilon=10^{-14}$ and maximal bond dimension is denoted as $\chi=256$. We will only investigate the mixed field Ising model with hopping $J$ and external fields $h_{X}$ and $h_{Z}$ as defined in (21), and probe the OTOC with Pauli $Z$ operators ($O_{1}=O_{2}=Z$ in (1)). Scrambling characteristics are then determined by least-squares fitting of $\ln\mathcal{C}$ at the wavefront to the expression (2).

The wavefront is determined as follows. First, due to numerical truncation with $\varepsilon=10^{-14}$ only data with $\ln\mathcal{C}>-22$ will be used. This delimits the right end $r$ of the wavefront; the default left end $l_{0}$ is then defined as the position where $\partial_{x}\ln\mathcal{C}$ is half the value at $r$. To eliminate the arbitrariness of $l_{0}$ a hyperparameter $\delta>0$ is introduced and the left end $l\equiv r-(r-l_{0})\delta$. When $\delta=1$, $l=l_{0}$ and when $\delta=0$, $l=r$; hence $\delta$ tunes the range of the wavefront, ending at $r$.

As a sanity check our implementation is verified against Exact Diagonalization (ED), which may be regarded as the MPO approximation with no bond dimension restrictions ($\chi=\infty$). The result is shown in Fig. 4. From the figure the MPO algorithm matches with ED perfectly at times before maximal bond dimension restriction is reached and starts to deviate afterwards. However, as shown in the figure, the wavefront dynamics is well captured by the MPO approximation, even after the bond dimension is saturated inside the butterfly cone. Such effectiveness of MPO (at least at infinite temperature) is observed in [19] and explained by the fact that at the wavefront the operator $O_{1}(0,t)$ is less complex, so only a smaller bond dimension is necessary.

A careful error analysis is necessary to extract reliable information from the nonlinear fit to the five parameters $(C,\lambda,x_{0},v_{B},p)$, appearing in (2). Here $C$ is the prefactor. Three major causes of systematic errors are identified: finite bond dimension $\chi$, a finite time range $[t_{0},t_{1}]$ of data and inaccuracy of the functional form (2). The convergence with respect to bond dimensions is verified: for all data used the difference in $\ln\mathcal{C}$ between $\chi=256$ and $\chi=512$ is less than $0.05$ and our main results do not depend on such a small difference. Also the fitting as presented in Fig. 5 is visually reasonably good, even for the chaotic Hamiltonian $h_{Z}=0.3J$ at low temperature $\beta J=3$.

The effect of a finite range of data and inaccuracy of the functional form is quantitatively manifested as dependence on the hyperparameters $\delta$ and $t_{0}$. Since the butterfly velocity is defined in the late time limit, $t_{0}$ should not be too small; but because only data up to time $t_{1}$ are available, $t_{0}$ cannot be arbitrarily large either. Moreover, larger $t_{0}$ means less data and more significant numerical instability. In Fig. 6, dependence on $\delta$ and $t_{0}$ of the fitted butterfly velocity for $\beta J=3$ and $h_{Z}=0.4J$ is shown. We will work with the values $\delta=1.0$, $Jt_{0}=1.5$ and $Jt_{1}=4.4$.

With this choice of hyperparameters, we produce the figures in the main text. Errors are estimated via slightly tuning hyperparameters. Details are summarized in Fig. 7, with fitted values of $p$ and $\lambda$ given as well. From the plot errors are estimated to be within a scale of $0.05$, $0.05$ and $0.5$ for $v_{B}(\beta)/v_{B}(0)$, $p$ and $\lambda/J$ respectively.

The correlation length $\xi$ is extracted with MPO numerics as well, as the inverse spatial decay rate of connected two-point correlations ${\rm tr}(\rho Z_{15}Z_{15+x})-{\rm tr}(\rho Z_{15}){\rm tr}(\rho Z_{15+x})$ in an $N=50$ chain with operator insertions at sites 15 and $15+x$, where $x=0,1,\ldots,20$. The exponential fit is remarkably good with correlation lengths at different temperatures and longitudinal fields shown in Fig. 8. Given the correlation length $\xi$ along with $p$ and $\lambda$ from Fig. 7, the bound is evaluated (with error estimates) in Fig. 9. In evaluating the inequality (Bounding the butterfly velocity) we have used $v_{S}\leq v$ for $v=3Ja$ and $\xi_{\mathrm{LR}}=a$ (cf. section C), where a Lieb-Robinson inequality with $(v,\xi_{\mathrm{LR}})=(3Ja,a)$ is verified in numerics and $a$ is the lattice spacing.