Quantum Scrambling and State Dependence of the Butterfly Velocity

Operator growth in spatially local quantum many-body systems defines a scrambling velocity. We prove that this scrambling velocity bounds the state dependence (for example, the temperature dependence) of the out-of-time-ordered correlator in local lattice models. This amounts to a basic constraint on quantum chaos from locality. Thus, for example, while operators that do not grow with time can have an associated butterfly velocity, that velocity is necessarily temperature-independent. For scrambling operators, in contrast, the butterfly velocity shows a crossover from a microscopic high temperature value to a distinct value at temperatures below the energy gap.

defined for local operators O 1 , O 2 in state ρ. The OTOC has been found to reveal a 'light cone' spread of quantum chaos, with two state-dependent chaos characteristics: the quantum Lyapunov exponent λ and the butterfly velocity v B . Just outside the light cone (or 'butterfly cone') |x| v B t for t > 0, the OTOC grows as the front is approached according to [18,19]: 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 chaotic 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 ∼ 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. This will allow us to understand the temperature dependence of the butterfly velocity in quantum chaotic spin systems, extending previous infinite temperature results [19,28]. The temperature dependence of chaos 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 with temperature is bounded by the scrambling velocity. This result is rigorous for one-dimensional systems and plausibly true more generally. A corollary is that while an operator that does not grow (with v S = 0) formally has a butterfly velocity, this velocity cannot depend on temperature. Such operators are present in non-interacting systems.
In contrast, the butterfly velocity of scrambling operators will typically be temperaturedependent. Turning our result around, we see that measurements of the butterfly velocity as a function of temperature can be used to quantify operator growth in local lattice systems.
We go on to determine numerically the temperature dependence of the butterfly velocity in certain one-dimensional chaotic spin chains. The results are shown in Fig. 2 below. The non-interacting transverse field model is indeed seen to have a temperature-independent butterfly velocity whereas the velocity is temperature-dependent for the chaotic mixed field models. In these curves, the butterfly velocity crosses over from a microscopic infinitetemperature 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 Λ of spins (or fermions) with a local Hamiltonian where h x are operators localized near lattice site 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. O(x, t) denotes O translated by a lattice vector x in space and a time t with Heisenberg evolution, and · is the operator norm. The causal light cone defined by v LR is such that for all v > v LR the norm decays exponentially at late times, so that λ(v) < 0. Therefore we can define v LR as the largest velocity such that the norm does not decay along a ray: v LR ≡ sup v : lim We shall not keep the dependence on direction n and operators O 1 , O 2 explicit.
For any v > v LR there are (v-dependent) constants ξ LR , C LR > 0 such that for all t, x > 0, Intuitively, inequality (5) states that for v > v LR , the norm [O 1 (0, t), O 2 (xn, 0)] is exponentially small outside the ray x = vt, with a tail of length ξ 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 ρ.
The 'velocity-dependent Lyapunov exponent' is defined by the late time growth or decay of the OTOC along a ray [18]: Analogously to the Lieb-Robinson case, the butterfly velocity can now be defined as which is state-dependent. The operator norm bounds the OTOC and hence 0

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'.
The role of the scrambling velocity will be to differentiate non-scrambling (e.g. free) and scrambling (e.g. chaotic) dynamics. For example, consider a free real scalar field φ(x). Because the equation of motion is linear, the solution takes the form φ(x, t) = dy f (y − x, t)φ(y, 0), for some real function f (x, t). Hence although the support of the operator φ(x, t) spreads out as t increases, it remains a superposition of local operators. In contrast, in random unitary circuits [33][34][35][36], generic operators do not stay in such a form for a long time; instead they quickly grow into a superposition of product operators with radius ∼ v LR t. We would like to say that free theories contain operators for which the scrambling velocity v S = 0 but The following definition captures the intuition of the previous paragraph. Given local operators O 1 and O 2 , the commutator [O 1 (0, t), O 2 (x, 0)] will grow along the ray x = vt.
We are interested in the growth of the operator itself rather than its norm or OTOC. Let R(x, t) be the radius of support of the commutator 1 and define 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. Coming back to the free field example, let (O 1 , O 2 ) be conjugate pair (φ, π). From the . This is a c-number and its support has radius R(x, t) = 0. Hence v S (v) = 0 for any v. In the random circuit, let O 1 and O 2 be two single-site operators. Inside the Lieb-Robinson cone, i.e. for |x| ≤ v LR t, For general systems and for |v| ≤ v LR we expect that

Scrambling bounds the state dependence of the OTOC
We have seen that in a free theory the commutator [φ(0, t), π(x, 0)] is a c-number. It follows that (i) this operator does not grow and (ii) its OTOC does not depend on the state ρ.
We will now see that this is an extreme case of a more general relation: the derivative of the OTOC, with respect to parameters such as temperature or chemical potential, is upper bounded by the scrambling velocity (8).
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: where β is the inverse temperature, a the lattice spacing, ξ the correlation length, ξ LR the microscopic lengthscale in (5), essentially the interaction range, and h ≡ 2 sup x∈Λ h 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 ∼ v B , this bound implies that the growth of the commutator at the butterfly light cone bounds the change of chaos characteristics such as the butterfly velocity. As (for example) the temperature is increased, these growing operators are 'activated' and contribute to the spread of chaos.
A generalization, with full proof in the appendices, is as follows: For any Gibbs state and C i = x∈Λ c i x is a sum of local operators, then The definition of c i > 0 is similar to h above:

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 ρ = e −βH /tr e −βH be a thermal state with inverse temperature β and correlation length ξ. 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.
(i) Temperature derivative of the OTOC. Taking the derivative of the OTOC (1) with respect to the inverse temperature gives where where δ > 0 can take any value. As for H, h y ≡ h y − tr(ρh y ). This decomposition can now be inserted into the derivative (11).
(ii) Dominance by operators inside the commutator. We first bound the contribution from outside of the support of the commutator, with |y − y 0 | > R + δ in (12). Due to the thermal correlation length ξ, the connected correlation function of h y with O † O will decay exponentially in the distance |y − y 0 |. Thus, for some constant C > 0 and all y ∈ Λ such that |y − y 0 | > R: Summing over |y − y 0 | > R + δ, the contribution to (11) from operators outside of the commutator is bounded by In d spatial dimensions and for R + δ ξ, C ∼ Cξ(R + δ) d−1 /a d from doing the sum over |y − y 0 | > R + δ (a is the lattice spacing). There is a technical subtlety in obtaining (13) due to the need to commute factors of √ ρ through h 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 |y − y 0 | ≤ R + δ. As in the main text, define the maximal local coupling in the Hamiltonian as Note that h y ≤ 2 h y , so that Notice that the inequality still goes through if we take Now, the number of terms in the first sum of (12) is V R+δ , the number of lattice points in a ball of radius R + δ. Therefore, putting together (13) and (15), we can bound the derivative (11) by: 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. Lyapunov exponent, C(vt, t; ρ) ∼ e λ(v;ρ)t as t → ∞. We furthermore set δ = (−ξλ(v; ρ)+ )t > 0, with > 0 a small number. This choice is such that the final term in (17) decays exponentially faster than the others as t → ∞. This final term is therefore negligible in this limit. In this way, as t → ∞ the following inequality is obtained: 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 ≈ 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 ξ → ξ + ξ LR . We include this shift in the following statement of the bound.)

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  (20)), via the shortest path from the origin to x. Hence it is plausible that the operator approximately one-dimensional, along the line connecting 0 and 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 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 [19,28] 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 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, it is possible that v S > 0 by our definition and the butterfly velocity is state-dependent [38].
Indeed, we have verified numerically that the butterfly velocity is temperature-dependent in such models. 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 free theory can have v S > 0 because entangled pairs of quasiparticles moving in opposite directions technically lead to a linearly growing radius of support for the operator. The reason the corresponding butterfly velocity is temperatureindependent in free theories, despite the existence of such entangled pairs, is as follows. In a free 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 This precise form for λ(v; ρ) 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 λ, v B and p ≥ 0 are state-dependent. Therefore, the ∂ µ i derivative in (10) will act on each of these quantities. Substituting the specific form (22) to the following slightly complicated expression: where ∆v ≡ v/v B − 1 > 0 is a dimensionless measure of how far the velocity is outside the butterfly cone. A simple consequence of (23) follows, however, when there is no scrambling.
Suppose that v S (v) = 0. In that case, taking ∆v → 0 + , the leading term on the left side of (23) is the last one. It follows that Hence v B is constant if the system is not scrambling.
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 ∆ (of order J for h Z = 0.1 ∼ 0.5J), as we now explain.
Intuitively, one might expect v B to cease varying at temperatures T ∆. 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 where . This is not an especially tractable expression in general, but it can be evaluated for a gapped system at zero temperature, in gapped systems at low temperatures, we may set h ≈ 0 in the bound (9). It follows that ∂ β v B → 0 when T → 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 ∂ β v B decays exponentially as β∆ → ∞. 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   (23) and the leading term on the left hand side is again the final one, which integrates to Here α is a dimensionful constant, γ a dimensionless constant and we have singled out the v S and ∆ dependences. It follows that (i) as v S → 0, ln v B can vary as a power v γ S of the scrambling velocity, and (ii) if the gap ∆ → 0, v B may approach zero at T = 0. Indeed, power law butterfly velocities v B ∼ 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.

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.

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 lattice 2 Λ of spins in d spatial dimensions, and rigorously proved for d = 1. There are isomorphic finite-dimensional Hilbert spaces H x associated to each lattice site x ∈ Λ and denote B x as the space of linear operators acting on To better characterize the spatial distribution of operators, define superoperators P S and Q S ≡ Id − P S such that P S is the projection onto the subspace x / ∈S CI ⊗ x∈S B x . That is, P S projects onto operators supported on S (so P S

More explicitly
where the integral is Haar averaging over unitaries outside S. However, note Q S is not the Henceforth if the subscript S = {x} is a single-element set, P {x} and Q {x} are written as P x and Q x for short. Also define the superoperator P r T with a superscript r > 0 as P S for S = {y ∈ Λ : ∃ x ∈ T, |x − y| < r}, i.e. projection onto operators supported within a distance r from the set T , and Q r T ≡ Id − P r T .
2 Technically 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. From (28) we have the following inequalities: We will be interested primarily in operators that are "exponentially local", denoted as B(x, R; ξ, C). We say O ∈ B(x, R; ξ, C) with x ∈ Λ, R, C ≥ 0 and ξ > 0, if for any r ≥ R, Intuitively, this means O is supported on the ball of radius R and centered at x, up to an exponential tail of lengthscale ξ. 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: which also defines R H > 0 and α labels different couplings in the Hamiltonian. Translational invariance is not necessary but h α ≡ sup x∈Λ h α x should be bounded. A Gibbs state is a density matrix of the form for some µ i ∈ R and In the proof it is not required that [C i , C j ] = 0. With only one i, with µ the inverse temperature and with C = H, ρ is the thermal density matrix.

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, ξ LR , C LR > 0, dependent on lattice geometry and Hamiltonian, such that (i) for any t ∈ R, r > 0 and operator O, where |∂S| is the number of lattice links (say, between x and y) such that x ∈ S but y / ∈ S; (ii) for any t ∈ R, operators O 1 and O 2 , In this bound v ∼ α |J α | h α R H , recall (31), i.e. coupling times range of local terms in the Hamiltonian, and ξ LR ∼ 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].
Note the theorem is specific to one dimension [44] and l A (µ i ) may be exponential in |µ i |; in this sense the restriction is weaker for complex time evolution: Theorem 2 (Araki). In one dimension, for any Gibbs state ρ as defined in (32) but with , ξ A > 0, dependent on lattice geometry and charges C i , such that for any finitely supported operator O and r ≥ l A (µ i ), where | supp O| is the number of sites in supp O.
Note, however, from the proof of the Araki bound (e.g., Theorem 3. 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) ρ, the correlation length of ρ is the ξ > 0 that is optimal with respect to the following property: there exists C > 0 and a function l 0 (·) > 0 such that for any operators O 1 and O 2 (supported on sets S and T ) sufficiently far apart, i.e., d ≥ l 0 (δ), where δ ≡ min{|∂S|, |∂T |} is the number of lattice links crossing the boundary of S or T , and d ≡ min{|x − y| : x ∈ S, y ∈ 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 ξ > 0 with the property stated around (38) has been proved for (i)

one-dimensional Gibbs states [43] (restricted to local operators
where |0 is the unique ground state of a gapped Hamiltonian [40,45], and (iii) thermal states ρ ∝ exp(−βH) in general dimensions at sufficiently high temperatures [46] (clearly ξ → 0 when β → 0). Of course the Hamiltonians associated with these states must be local, as in (31) above. It is plausible that the correlation length ξ as defined around (38) is finite for Gibbs states ρ in general systems with local dynamics and away from phase transitions. and Hamiltonian, such that for any t > 0, x > 0,

From Theorem 3 one immediate candidate for defining the Lieb-Robinson velocity is
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 ξ LR and C 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] v We will assume that the limit exists and is a continuous function of v. By definition v LR 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 LR :

Proposition 1. For any direction n and operators
Proof. Let v > 0 belong to the set in (40), i.e., there exist ξ, C > 0 such that for all , and hence any v > v is not contained in the set in (41). Therefore the supremum v (2) LR is at most v. This is true for any v > 0 in the set of (40), hence v LR , we need the following lemma:

Lemma 1. For any positive functions f (x, t) and g(x, t), if limits
Proof. Because the limits (42) are uniform, for any ε > 0 there is T (ε) > 0 such that for LR (n; O 1 , O 2 ), given the limit in (41) is uniform for all v > v (2) LR (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) LR , we show that (i) implies (ii), and (ii) implies (iii), where Step two: By definition (41) we have for any v > v (2) 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 LR . So now it remains to prove that (i) ⇒ (ii) and (ii) ⇒ (iii): says that λ(v ) < 0 for any v ≥ v and to arrive at (ii) we hope to find ε, ξ > 0 such that Before construction of ε and ξ, it is remarkable that there is a restriction on λ(v ) from Theorem 3: the Lieb-Robinson bound states that there are some C 0 , v 0 , ξ 0 > 0 such that we have to check that −λ(v ) > 0 is bounded from zero on [v, ∞). The only concern is λ(v ) may be arbitrarily close to zero when v → ∞; but this is not possible because from the Hence ε > 0 is well-defined in this way.
+ε)} (as constructed in the last paragraph the denominator is always negative). The task is then to show that ξ < ∞; similarly the only place things could go wrong is when v → ∞, but in that limit /ξ for all x ≥ vt and t ≥ t 0 . Hence for (iii) to hold it suffices to choose that C ≡ max{2, sup 0<x<vt 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 ≥ v 0 t, and (c) 0 < t < t 0 with 0 < x < v 0 t separately.
can be bounded using the Lieb-Robinson Theorem 3: there is some C 0 , v 0 , ξ 0 > 0 such that 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 ξ appearing in (ii), (iii) is true.
Henceforth the Lieb-Robinson velocity will be defined as v LR ≡ v where the OTOC C O 1 O 2 (x, t; ρ) is defined in (1). As the velocity-dependent quantum Lyapunov exponent is defined as in (6), an equivalent definition of v B reads: As expected, the butterfly velocity in any state is bounded by the Lieb-Robinson velocity: Proof. This follows from definition (41) and (44), and Finally the scrambling velocity can be precisely defined in the language of exponentially local operators, defined around (30).
where the smallest ball, with radius R and centered at x, is understood as roughly the "support" of the commutator O. The quantities ξ and C characterize the exponential tail that we neglected in the main text. Clearly v S ≥ 0 and decreases with increasing ξ.
. But the first term is supported in the ball of radius r centered at origin, so All velocities can be maximized over direction n to recover their isotropic definitions, or over O 1 , O 2 ∈ O to remove the operator dependence.

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 ≥ 0 and k, γ > 0, where x 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 ≥ 0, γ > 0 and integer m ≥ 1, and the inequality (47) follows because if m = k , m ≤ k + 1 in the linear factor and k ≤ m implies that e −γm ≤ e −γk as well.
The following lemma bounds the product of two exponentially local operators: Let O 1 ∈ B(x, R; ξ 1 , C 1 ) and O 2 ∈ B(x, R; ξ 2 , C 2 ), then for any r ≥ R, Proof. Note that for any r > 0, (29) and (30), for r ≥ R, Next is the Araki bound (cf. Theorem 2) for exponentially local operators: Theorem 4. For any one-dimensional Gibbs state ρ as defined in (32) with µ i ∈ C and operator O ∈ B(x, R; ξ, C), there exists C (µ i , ξ, C) > 0 (dependent on lattice geometry and Here l A (µ i ) and ξ A are those appearing in the Araki bound, and a is the lattice spacing.
Observe that the operator ρOρ −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: Proof. First note that for ζ(ξ) ≡ ξ A + ξ, so it suffices to find C (µ i , ξ, C, ε) such that for all 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 ρ be a one-dimensional state with ξ, C and l 0 (·) > 0 as stated around (38).

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 ρ as defined in (32) with correlation length ξ cor (read around (38) for a definition), ε, δ > 0, any operators O 1 , O 2 and x ∈ Λ, and where a is the lattice spacing and ξ A is defined in Theorem 2. The inverse temperature is denoted as β and J α labels couplings in the Hamiltonian (31).
R, ξ and C are such that O ∈ B(y 0 , R; ξ, C) for some y 0 ∈ Λ. Finally where c i y ≡ c i y − tr(ρ c i y ), and same for h α with c i y replaced by h α y . And if C i commute with each other, c i can be chosen as Proof. We start with proving (70). By definition (1) and (32), where for any operator C, C ≡ C − tr(ρ C). Now recall C i is a sum of local terms (33): for any S(r) ≡ {y ∈ Λ : |y − y 0 | ≤ r}. For any y ∈ Λ, by definition of c i (s), The inequality (76) is good enough for the terms in S(r). For the remaining terms with y away from 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 ρ are separated by -and do not necessarily commute with -the c i 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 † O ∈ B(y 0 , R; ξ, 4C) and from Theorem 4 and Corollary 1, there is C 1 (µ i , ξ, C, ε) > 0 and l(µ i ) > 0 such that for any 0 ≤ s ≤ 1, Hence by Theorem 5, because tr(ρ c i y ) = 0, for any 0 ≤ s ≤ 1, where C 2 is defined in terms of the prefactor C cor (µ i ) in (38) as, using (77), Now bound the sum (75) by choosing r = R + l(µ i ) + a + R H + l 0 (2) + δ and apply (76) for y ∈ S(r) and (79) for y / ∈ S(r), (denote ζ ≡ ξ cor + ξ A + ξ + ε) and use the inequality (80) and c i y ≤ 2 c i y to reduce to the form (70). Proving (71) is essentially the same except in the first step: there is an additional term due to coupling dependence of O 1 (0, t). By definition, and (71) follows from the Cauchy-Schwartz inequality for the inner product . Finally if C i commute with each other, the first step (74) can be replaced with The theorem, as stated, seems complicated; but the physics is much clearer in terms of the velocity-dependent Lyapunov exponent (6): (46), Proof. Divide both sides of (70) by t x = vt and take the limit t → ∞ (assuming the limit and derivative commute): Finally let ε, ξ A → 0 to conclude 5 .
The operator O must decay at large distances at least as quickly as the rate set by ξ LR (appearing in any triple (v, ξ LR , C LR ) with a Lieb-Robinson bound Theorem 1). Therefore we take ξ = ξ LR in the main text. We have already noted in section C that this then defines The coupling dependence of λ O 1 O 2 (v; ρ) can be bounded in the same way: for κ 1 , κ 2 > 0 and h α ≡ sup y∈Λ h α y at t → ∞,  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 ∂ J v B via an inequality similar to (23). Cusps of chaos characteristics are indeed observed at quantum critical points in e.g. [47,48].

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.  (21), and probe the OTOC with Pauli Z operators (1)). Chaos characteristics are then determined by least-squares fitting of ln C at the wavefront to the expression (2).
The wavefront is determined as follows. First, due to numerical truncation with ε = 10 −14 only data with ln 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 ∂ x ln C is half the value at r.
To eliminate the arbitrariness of l 0 a hyperparameter δ > 0 is introduced and the left end l ≡ r − (r − l 0 )δ. When δ = 1, l = l 0 and when δ = 0, l = r; hence δ 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 (χ = ∞). The result is shown in Fig. 4   In the second panel, each curve is a spatial profile of the OTOC at a fixed time. Propagation of a butterfly wavefront is clearly observed. For all χ the agreement with ED is remarkable until the MPO truncation ε = 10 −14 kicks in after ln C drops to approximately −25.
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  (2) to wavefront is marked as solid. Each curve is ln C for a fixed Jt = 0.2, 0.4, . . . , 4.8. The first plot is for β = 0 and h X = 1.05J, h Z = 0 with a fitting v B = 1.95Ja, p = 0.46 to be compared with exact values v B = 2Ja and p = 0.5 (a = 1 is the lattice spacing); the second plot is for βJ = 3, h X = 1.05J, h Z = 0.3J and the best fitting is v B = 1.39Ja with p = 0.65. to the five parameters (C, λ, x 0 , v B , p), appearing in (2). Here C is the prefactor. Three major causes of systematic errors are identified: finite bond dimension χ, 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 C between χ = 256 and χ = 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 βJ = 3.
The effect of a finite range of data and inaccuracy of the functional form is quantitatively manifested as dependence on the hyperparameters δ 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 δ and t 0 of the fitted butterfly velocity for βJ = 3 and h Z = 0.4J is shown. We will work with the values δ = 1.0, Jt 0 = 1.5 and Jt 1 = 4.4.  inequality with (v, ξ LR ) = (3Ja, a) is verified in numerics and a is the lattice spacing.