On entanglement spreading in chaotic systems

We discuss the time dependence of subsystem entropies in interacting quantum systems. As a model for the time dependence, we suggest that the entropy is as large as possible given two constraints: one follows from the existence of an emergent light cone, and the other is a conjecture associated to the"entanglement velocity"$v_E$. We compare this model to new holographic and spin chain computations, and to an operator growth picture. Finally, we introduce a second way of computing the emergent light cone speed in holographic theories that provides a boundary dynamics explanation for a special case of entanglement wedge subregion duality in AdS/CFT.


Introduction
In free field theory, the spreading of entanglement can be understood in terms of freestreaming particles [1,2,3]. In chaotic systems, this picture is no longer accurate [4,5,6]. We would like to propose a replacement, based on the idea that the entropy is as large as possible given two constraints: first, information must remain within an emergent light cone, and second, the rate of change of the entropy is bounded by a certain multiple of the area.
The analysis in this paper will revolve around three speeds that can be assigned to extended quantum systems. We introduce these now: v E the "entanglement velocity" was defined in [7,8,9] by the statement that at early times after a quench, the entropy of a large region grows as d dt S[A(t)] = v E s th area(A), where s th is the equilibrium entropy density.
v LC is the speed of the effective light cone, defined by some fixed and small threshold for commutators [W (t, x), V (0)] 2 β , where V, W are arbitrary operators that do not change the energy density drastically, and the expectation value is taken in the canonical ensemble with inverse temperature β. 1 v B the "butterfly effect velocity" was defined in [11,12] as the speed at which the region where the commutator is order one expands outwards.
One expects a scrambling time delay between the effective light cone (where the commutator is small) and the "butterfly cone" (where it is order one). However, we suspect that generically the slopes are the same, so that v LC = v B . This is true in holography and likely to be true in a weakly coupled field theory, using an analysis along the lines of [13]. We caution the reader that in the spin chain we study below it does not appear to be precisely true. After the first version of this paper appeared, 2 [14] found that in the Bose-Hubbard chain v LC is significantly greater than v B . Although we have been careful to distinguish them in this introduction, we will use v B and v LC somewhat interchangeably in what follows. It would be interesting to refine our analysis in a way that distinguishes between the two velocities.
The main point of this paper will be to associate each of the speeds v E and v LC to a bound on the time dependence of the entropy, and to compare to data. The two bounds interact in a somewhat nontrivial way, and we believe that many coarse-grained features of entropy and information dynamics in systems with v LC = v B can be captured by a simple model in which we simply saturate the combined bound. We will compare this model to new data from holography and from chaotic spin chain numerics. Our work should be understood as a combination of the tsunami picture [7,8,9,4] and the models of [6,15].
The plan of the paper is as follows. In section 2 we motivate the two bounds, one of which is essentially [16]. In section 3 we show with new computations (presented in detail in [17]) that the actual holographic results are close to the combined upper bound. In particular, certain holographic systems saturate entropy as fast as allowed by the bounds, at least for spherical and strip entangling regions. In that section we also discuss bounds and holographic data for a mutual information quantity that tracks the spread of information in a more detailed way than the standard entropy. In section 3.3, we make a similar comparison to numerical data from a simple chaotic spin chain.
In section 4 we discuss an operator growth model that saturates the entropy inequalities. In this model, v B = v LC is the rate at which operators tend to grow, and v E is related to the small probability that they actually don't grow at all.
In section 5, we make a comment that relates the speed v B to a special case of entanglement wedge reconstruction. Finally, we conclude and provide a brief outlook. Further numerical data about the chaotic spin chain, the holographic computation relating v B to entanglement wedge reconstruction in higher derivative gravity theories, and the comparison of the bounds to the quasiparticle model are collected in the appendices.
2 Two constraints on the time dependence of entropy In this section we will argue for two inequality constraints on the time evolution of entropy. We assume that the state of interest has an approximately uniform energy density that allows us to assign an effective temperature T = β −1 . We will consider large entangling regions with characteristic size r β, see [18] for a recent study of small regions. The first constraint comes from the idea that the system will have an emergent light cone at speed v LC . Concretely, this means that operators will approximately commute with other operators if their spatial separation is more than v LC times their time separation. This imposes a speed limit for entanglement [6,16] by an argument exactly parallel to that of [16]. The argument goes as follows: any subsystem A (t ) at time t < t is actually a subsystem of A(t) at a later time t provided that the future "v LC cone" of A (t ) passes entirely through A(t), see figure 1.

A(t)
A'(t') Figure 1: Left: The "v LC cone" for a disk A (t ) passes through A(t). Any information in A (t ) is contained in A(t), hence A (t ) is a subsystem of A(t) and we can use the monotonicity of relative entropy (2.1). Right: The tsunami wavefront for different times corresponding to the case where ∂A is an ellipse (purple). Lighter color means later time.
Monotonicity of the relative entropy then implies is the reduction of the thermal density matrix (at the effective temperature of our state) to region A. It was pointed out by [16] that if the region A is much larger than the thermal scale β, then S(ρ A |ρ (th) We can understand this inequality intuitively. I 1 [A(t)] is a measure of the amount of information contained in the density matrix ρ A (t). The inequality is just the statement that if the subsystem A (t ) is contained in A(t), then it contains no more information than A(t) does. The parameter that enters here is v LC , the slope of the effective light cone. This will depend on the system (and, in general, state) of interest. However, it can be established by an independent calculation, such as the commutator of local operators. If we consider this bound alone, we get where S[A(0)] subtracts the entropy of the initial state to make the LHS finite even in field theory, and vol(tsunami(t)) is the volume covered by a tsunami wave in time t propagating in from ∂A with speed v LC [16], see figure 1. 3 The second inequality we propose is that the time derivative of the entropy of any region A should be bounded by the area: where s th is the equilibrium thermal entropy density at the energy/temperature scale we are working. 4 The parameter v E is defined by requiring equality in (2.4) for large regions and short times after a quench. 5 So the content of (2.4) is that the entropy should never change faster than right after a quench. Unlike (2.2), this is simply a conjecture. There are rigorous bounds of the form (2.4) but with a coefficient that depends on the Hilbert space dimension and operator norms [19,20]. We expect that in practice there might be a better bound that applies even to theories with unbounded operators, provided that the energy is not too high. (2.4) is a guess for what this bound might be. Moreover, (2.4) can be proven in holographic theories [17].
We will now make three comments about the above constraints. First, as noted by [16], in the setup of a global quench, (2.2) implies (2.4) with v E replaced by the light cone speed. (In [16] this was taken to be c.) It follows that v E ≤ v LC . (2.5) Indeed, for the particular example of charge neutral high-temperature states of holographic theories with boundary spacetime dimension d, we have [7,8,9,11,12] which satisfy (2.5). In [17] the null energy condition is used to verify (2.5) for general states in holographic theories. In the spin chain we will discuss in section 3.
Second, if we know the entropies of all subregions at a given time, we can use (2.2) and (2.4) together to give upper bounds on the entropy at later times (or equivalently, lower bounds on the information). We will see that in many cases, the optimal bound involves both conditions. It is interesting to compare this bound with the actual answer for the entropy. In the next section we will make this comparison for holographic theories and for a simple chaotic spin chain.
Third, the bounds inform us about the saturation time t S , when the entropy reaches the thermal value. For general shapes, it is somewhat nontrivial to evaluate the best bound that follows from (2.2) and (2.4). But it is easy to see that we will have the inequalities .
Here r insc is the radius of the largest ball that fits inside A. The first of (2.7) is a better bound for "round" shapes, and is easy to derive from (2.3). The second is better for elongated shapes and follows from integrating (2.4).
3 Comparing the bounds to data

Global quench
We begin by studying the entropy S[A(t)] following a global quench from a sparsely entangled state. Then we set the initial conditions S[A(0)] = 0 for all regions A, and we compute the upper bound on S that follows from (2.2) and (2.4). 6 An interesting test case is to consider A to be a sphere of radius r A . Then we use (2.2) with A chosen to be a sphere of radius 0 ≤ r ≤ r A at time t = t − (r A − r)/v LC > 0. This gives We further bound S[A (t )] using (2.4). All together, To get the best upper bound, we minimize (3.2) over 0 ≤ r ≤ r A , finding  If v E < v LC d−1 , then the v LC constraint plays no role and we have linear growth all the way until saturation, see figure 2.  until saturation for a large sphere. Black/dashed is our upper bound using the speeds (2.6), blue/dotted is the bound from relativistic causality [16] and red/solid is the holographic computation [17]. The dashed lines should be understood as one-parameter (v E ) fits, because v LC is determined by an independent calculation based on commutators (or equivalently, out-of-time order correlation functions).
It is also straightforward to consider shapes other than spheres following the logic of (3.2). For strips only the bound (2.4) is relevant before the saturation time, and we get This bound is saturated for strips in holographic theories for all times, so the saturation time bound (2.7) applicable for elongated regions is also saturated [7,8,9]. It would be interesting to work out the upper bound for other shapes.

The thermofield double state
The increase of entropy after a quench is a result of the delocalization of information under time evolution. We can map this out in more detail by considering another entropy quantity, the two-sided mutual information [21,7] in the thermofield double state The state |T F D(0) is a thermal version of a maximally entangled state, with a high degree of entanglement between large subsystems A R ⊂ R and the corresponding subsystem A L ⊂ L. It will be convenient to think about R as the physical system, evolving in time, and L as a static reference system that keeps track of the movement of information in R. In other words, we view the time evolution as |T F D(t) = e −iH R t |T F D(0) , where H R is the Hamiltonian acting on the R system. The key point is that the mutual information can be understood as a measure of the amount of information originally contained in A R (0) that is now contained in B R (t). Note that this is a statement about only the physical R part of the system; the thermofield double is a tool to make this notion precise.
We will think about I[B R (t), A L ] for the case where A, B are concentric balls, of radius r A and r B . The mutual information is then a function of three parameters: t, r A , r B , which we take to be large compared to the inverse temperature β. It is interesting to fix t, r A and consider I as a function of r B . This tells us how much of the information originally contained in a ball of radius r A is later contained in a concentric ball of radius r B .
In the |T F D(t) state, the first two terms in (3.6) are exactly given by the thermal values. The third term can be computed in holography by applying the Ryu-Takayanagi (RT) formula [22] and the techniques of [17]. We can also compute a bound on this quantity using an analog of (2.2) for mutual information, namely . This inequality, together with (2.4), and the initial conditions I[A L (0), B R ] = 2s th · vol(A ∩ B) lead to the lower bound on I, for any 0 ≤ r ≤ min{r A , r B } such that the expression in parentheses is positive. We get the best bound by maximizing over such r values. In the above expression Ω d−2 is the volume of the unit S d−2 . The maximization can be done explicitly, but rather than giving the answer in the various cases, we will just plot the result against the holographic answer. We show this in figure 4. The holographic answer always lies a little above the lower bound, but the curves are quite close. This means that the holographic answer is almost as small as it can be, given v LC and v E .  As a function of r B , there are two interesting points: the value where the mutual information first becomes nonzero, and the value where it saturates. In both the bound and in holography, these are given in terms of v LC as This is illustrated in figure 5.
) It is clear that the mutual information should saturate for r B ≥ v LC t + r A , since then the emergent light cone implies that B(t) contains A(0) as a subsystem, see figure 5. However, the fact that the mutual information is zero until r B = v LC t − r A is not obvious, and is a striking difference from a free field theory. 8 8 For this to hold, we have assumed that the holographic dual black hole geometry is such that in the quench setup giving the same final state black hole the entropy saturates "continuously". This is not true for the 2 + 1-dimensional case, but the violation of the first equality in (3.9) is invisibly small on figure 5. (The second equation in (3.9) remains valid.) More discussion of these issues is found in footnote 7, section 5, and [17], where examples with large violations of r first nonzero

Spin chain data
In this section we repeat the analysis of the previous sections for a simple chaotic spin chain with the Hamiltonian (3.10) First we consider the evolution of the entropy following a "quench." This was considered for essentially the same model by [23]. By representing the Hamiltonian as a sparse matrix and evolving the state directly, we can avoid exact diagonalization and study a reasonably long chain of n = 26 spins. We take the initial state |Y + , which thermalizes well [24]. The entropy shows clear linear growth until near saturation, with a fitted slope of v E ≈ 1 in lattice units, see figure 6. We note that if we had studied a smaller subsytem, the saturation would have been somewhat sharper. Rényi entropies are plotted in appendix A. We can also consider the mutual information quantity (3.6). In this d = 2 setting, the maximization of (3.8) over r is simple, and we find the lower bound where L A , L B are the respective lengths of the A, B subsystems, which are both assumed to begin at the end of the spin chain. Notice that this bound involves both v E and v LC in a nontrivial way. To compute the mutual information numerically, we use exact diagonalization to study an infinite temperature thermofield double of 14 spins on each side. Note that this is just the standard maximally entangled state. See figure 6 for a comparison to the lower bound from (3.8). Although the agreement is not as good as for holography, the actual mutual information is surprisingly close to the lower bound. It seems that the spin chain result starts to lag behind somewhat as we move to later times. This may be related to the fact that v LC appears numerically to be slightly larger than v B for this system.

An operator growth model
In this section we discuss an operator growth model that saturates inequalities analogous to (2.2), (2.4) for Rényi 2-entropy. The model is not derived from unitary time evolution, but it may help to give an understanding of v E and v B = v LC . The discussion is inspired by the model of [15], but the details are different. 9 The setting is as follows. We imagine that we have a system of N qubits arranged in a regular (d − 1) dimensional lattice. They start in the "quench" pure state |0 |0 ...|0 so that the density matrix of the complete system is We would like to compute tr A ρ A (t) 2 . We can think about computing this by studying the time evolution of an individual term in (4.1) i∈S σ (i) z (t) = subsetsS of sites and Pauli indices Here we are expanding the evolution of a given term in ρ as a sum of products of Pauli 9 The operator counting model of [15] gives v E = v LC , and in d > 2 it gives a saturation time that violates (2.7); our model improves on these aspects.
operators. The Rényi 2-entropy can be obtained from where A c is the complement of A, and we plugged in the definitions (4.1), (4.2). Because the Pauli operators are traceless, ifS 1,2 contains sites outside A, then tracing over A c gives zero. This restricts the sum toS 1,2 ⊂ A, and we can use the orthogonality of Pauli operators to obtain: where in the final step we neglected the off-diagonal contributions. For large regions, this can be justified under the assumption that the coefficients c S→S have random and uncorrelated signs. 10 Fixing S in the final sum corresponds to focusing on a single product operator in the expansion (4.1); the contribution of such a term can be understood as the "probability" that this operator remains inside the region A after time t. Here, in defining the probability we are thinking of the space of operators as a Hilbert space, with basis vectors the different Pauli strings. The Heisenberg evolution of an operator defines a quantum evolution on this state space. A simple model of this the probability is to assume that operators grow at a rate v B = v LC : the probability is 1, if the support of the operator (growing outwards from the initial support with v B for time t ) is inside A, and 0, if the support reached outside. Hence all operators supported inside the "dry region" not reached by the tsunami discussed around (2.3) contribute to the sum (4.4) with weight 1, and the rest of the operators have 0 weight. Then (4.4) becomes where s th is log 2 times the density of sites. This model leads to an entropy saturating the inequality (2.3) for all times. 11 10 In detail, let S andS denote the sizes of the sets from which we take S 1,2 andS respectively. Then there are SS diagonal, and S 2S off-diagonal terms. If we assume independent fluctuating signs for the off-diagonal terms, their contribution is proportional to ±S S , which can be neglected compared to the diagonal terms: Next, we introduce a refinement of the simple operator growth model discussed above. We assume that operators mostly grow outwards, at a rate v B = v LC . However, we also assume that there is some tail in the probability distribution for a given operator to remain the same size. This tail decays exponentially in time with a coefficient proportional to the area of the boundary of the operator: e −γ t·area . As an example, we can compute the probability that an operator initially filling a ball of radius r < r A will be within a concentric ball of radius r A after time t. The optimal strategy is for the operator to pay to stay the same size for a certain time t and then grow outward until time t, which leads to We now use this formula to evaluate (4.4) for the case that A is a ball of radius r A . We imagine that we are evaluating the sum for large regions, so we use a continuum notation. The concentric ball-shaped operators dominate this sum (4.4), and we find where the one is a special treatment of the identity operator, which does not grow. The Rényi 2-entropy is S 2 [A(t)] = − log tr ρ A (t) 2 . For large regions the integral is dominated by the maximum, and we recover the maximization over r of (3.2). In other words, this model saturates the entropy bounds for spherical regions. More specifically, we find that it saturates the bounds with v E = γ s th . Thus, the operator growth model provides a picture of entanglement spread that saturates the bounds. In this model, v B = v LC is the rate at which operators mostly grow, and v E is related to the small probability that operators actually do not grow. For spherical regions we get two complementary pictures of saturation: from the perspective of the bounds saturation occurs when the tsunami covers the whole entangling region, while in the operator picture saturation happens when the support of an initially small operator in the center grows bigger than the region. The agreement between the operator growth model and the bounds reinforces our confidence in the overall picture provided in this paper.

The emergent light cone and the entanglement wedge
In this section, we will give a second discussion of the emergent light cone from the perspective of information spreading. Rather than studying the dynamics of information following a global quench, we will examine the delocalization of a small perturbation to the equilibrium thermal state. The calculation is very simple, but it connects in an interesting way to the somewhat mysterious entanglement wedge [25,26,27] subregion duality [28] in AdS/CFT.
Concretely, the problem is the following: imagine acting with a light local operator W x on the thermal state. Initially, some information about which operator we applied can be recovered from a local measurement at position x. But as time passes, scrambling delocalizes the information over a larger and larger region. We can think about this as the growth of the operator W x (−t), where its size is just the smallest region that contains significant information about the applied operator. We will defineṽ B as the rate of growth as a function of time. One expects to findṽ In a theory with a gravity dual, we can evaluate this size using subregion duality, which asserts that certain subregions of the boundary theory completely describe corresponding subregions of the bulk. The initial thermal state is represented by a static black hole, and acting with W x introduces a small perturbation that falls towards the horizon. The smallest subregion of the boundary that contains significant information about W x after time t is simply the smallest boundary subregion such that the corresponding bulk subregion contains most of the falling particle's wave function. We will follow [25,26,27] and assume that the bulk subregion is the "entanglement wedge." For our purposes, this is just the region of a constant time slice of the bulk that is contained within the RT surface associated to the boundary subregion. Our analysis will focus on the near horizon region, where it is convenient to use a Rindler radial coordinate ρ that measures proper distance from the horizon. The nearhorizon metric of a general static planar black hole is then where β is the inverse temperature and r h is the area-radius of the horizon. In Schwarzschild or Rindler coordinates, the falling wave packet never crosses the horizon, but it is easy to check that as time advances it approaches exponentially, To evaluate the size of W x (t), we need to find the smallest boundary region such that the RT surface extends down to the radius in (5.2). The optimal choice is to take a ballshaped region of radius R, on a constant t slice in the boundary. The RT surface can then be parameterized by ρ(x i ), and by symmetry, it will depend on only one (radial) coordinate in the boundary. Finding the surface exactly would involve a nonlinear ODE, but for large boundary regions, most of the RT surface lies very close to the horizon, where the equation linearizes. 12 Indeed, near the horizon one can check that Solutions to the equation of motion ∂ 2 i ρ = µ 2 ρ arising from (5.3) will vary exponentially as a function of x i . Suppose that the radius of closest approach to the horizon is ρ min , and that it happens at the origin of the x i coordinates. Then the solution in the near-horizon region is [29] ρ( Eventually ρ will exceed β and the surface will exit the near-horizon region. From this point, it will reach the boundary within an order one distance in x. We can therefore determine the size of the operator R in terms of ρ min , to within an order one error, by solving where we ignored prefactor powers and constants in the approximate expression. Now we put these two steps (5.2) and (5.5) together. In order for the falling particle to be within the entanglement wedge, we need ρ min ≤ ρ(t), which implies that R ≥ṽ B t withṽ For the special case of an uncharged black hole in Einstein gravity, we getṽ B = d/2(d − 1). This matches (2.6), so we see that indeedṽ B = v B . We will make four further comments: (1) It can be checked that (5.6) agrees with v B derived from the shock wave computations for a general black hole background in Einstein gravity. We are not certain if v B = v B in general, but in appendix B we show thatṽ B = v B remains true in certain higher derivative theories, where we know how to do both computations. Also, the fact that saturation of I[B R (t), A L ] saturates at the point predicted by v B in the left panel of figure 6 is an indication thatṽ B = v B in the spin chain as well.
(2) Although we have presented this calculation as a second derivation of v B , one can turn the logic around and view these results as a dynamical explanation for a special case of entanglement wedge subregion duality [25,26,27]. Other proposals have been made for the bulk region that is described by a given boundary region, including the "causal wedge" [28,25,30]. That proposal would have led toṽ B = 1.
Let us clarify this point slightly. In order to get an operator deep in the bulk, we can start with a local operator W x and time evolve. To measure this operator in the CFT at a later time t, relativistic causality implies that we need no more than a ball of radius t in the field theory. If this were the only constraint, we would end up with reconstruction only within the causal wedge. However, because of the emergent v B = v LC cone, we can measure W x (−t) on a smaller subregion, of radius v B t. We've seen that the specific value of v B singles out the entanglement wedge. Notice that the boundary domain of dependence of this smaller subregion will nowhere contain W x (−t) as a local operator, but it will approximately contain it as a nonlocal operator. This possibility was previously conjectured to be relevant for subregion duality, at the very end of [31].
(3) We can consider non-spherical regions of the boundary, and ask whether they can be used to reconstruct the operator. It is an easy extension of the above to check that the falling particle will be within the entanglement wedge if and only if the region contains a ball of radius v B t surrounding the initial location of the operator. (More precisely, this is correct up to an imprecision of scale β in the boundary coordinate.) (4) Finally, we can also view the above calculation as giving the saturation time for the entropy of a ball-shaped region following a global quench in certain circumstances. The argument is simple. In the model of a quench discussed in [8,9] an infalling shell of matter creates the black hole, and we get a nontrivial time dependence for the entropy because the Hubeny-Rangamani-Takayanagi (HRT) [32,33] surfaces cross the shell. In the case of "continuous" saturation, saturation of entropy occurs at time t S , when the HRT surface climbs out from the AdS region behind the shell, and just barely touches the infalling matter. At this time, the geometry that the HRT surface is experiencing is that of a static black hole.
Note that at this point, the HRT surface becomes the RT surface that we have been discussing above. The role of the trajectory of the falling particle is played by the infalling null shell. It follows that the saturation time is t S = R/ṽ B . 13 In the holographic cases we have been able to analyze,ṽ B = v B = v LC , so we saturate the first bound in (2.7). It is interesting that the two setups appear superficially different: the falling particle is localized in the boundary theory spatial directions and its back reaction can be neglected, while the null shell is translation invariant and it creates the black hole geometry. For the problem we consider, however, the only thing that matters is that both follow the trajectory (5.2), and that the HRT surface only experiences the static black hole part of the geometry.
However, in certain black hole geometries it can happen that there are multiple extremal surfaces corresponding to a given boundary theory time, and according to the HRT prescription we have to pick the one with minimal area. In such a case, the saturation is "discontinuous": the HRT surface barely touching the null shell doesn't play a role in saturation, t S > R/v B , and the time derivative of the entropy at saturation is discontinuous. The conditions under which this happens is analyzed in detail in [17]. For uncharged black holes, saturation is discontinuous in d = 3 and continuous in d > 3.

Discussion
Our motivation for this work was to clarify the relationship between v E , the "entanglement velocity" and v B , the "butterfly velocity." We have argued that they are associated to two different bounds on the time evolution of the entropy.
These are not universal bounds, in the sense that they depend on parameters v E and v LC that will vary from system to system. But we found evidence from holography and (to a somewhat lesser extent) spin chain numerics that the true evolution of the entropy is well approximated by saturating the combined bound. We suggest that for chaotic systems with v B = v LC , this may be an appropriate replacement for the particle streaming picture of entropy dynamics in free field theories [1,6,3].
Added in v2: While we expect that generically v B = v LC , this does not appear to be precisely true in the spin chain that we studied. In the recent exciting paper [14], it was found that v LC is significantly greater than v B in the Bose-Hubbard chain. We suspect that these violations are due to the special kinematics in one spatial dimension, but the relation between the two velocities (or more generally the shape of the commutator [W (t, x), V (0)] 2 β ) deserves further investigations. When v LC > v B , we have to use v LC in the bound (2.2). Because a significant amount of quantum information spreads with the slower speed v B , the entropy is not expected to lie as close to the combined bounds, as in the case v B = v LC .
It would be nice to have a more precise understanding of v E . In the holographic dual, v E seems to be an interesting quantity, as it is determined by the geometry behind the horizon of a black hole. One hint comes from the operator growth model discussed in section 4, in which v B is related to the typical speed at which operators grow as a function of time, and v E is related to the small probability that they do not grow. Another possibility [7,34], is that a tensor network perspective may be helpful.
It would also be nice to know whether there are subregion shapes for which the holographic answer for the entropy is far from saturating the combined bounds discussed in this paper.

A More spin chain numerics
In this appendix we present some further numerical results in the chaotic spin chain defined by the Hamiltonian (3.10). First we discuss the growth of the Rényi entropies following a quench (see also [35]) from the initial product state |Y + . We evolve the state by directly applying the Hamiltonian to the state as a sparse matrix, avoiding the need for exact diagonalization. This allows us to study a 26 site chain, and we evaluate the entropy of the first 12 sites. The results are plotted in figure 8. If we had studied a smaller subsytem, the saturation (particularly of the higher Rényis) would have been somewhat sharper. We note that v E seems to depend on the Rényi index. This is reasonable given our interpetation (see section 4) of v E as being determined by a small tail in the probability distribution for the evolution of an operator. Although we have not attempted the calculation, it seems likely that v E would depend on the Rényi index in holography as well, since the geometry will be affected by the backreaction of the Rényi brane.
Next, we discuss the computation of v LC and v B in the spin chain. We choose to define v LC by a value of the commutator of 0.1, and we choose to define v B as a value of 1. For these choices we find v B ≈ 1.7 and v LC ≈ 1.85. The results do seem to depend (weakly) on these specific choices. The results are plotted in figure 9.

B Higher derivative corrections
In this appendix we show that the method of computing v B from shock waves andṽ B from the entanglement wedge agrees in four derivative gravity theories to all orders in the higher curvature couplings. These are theories where we not only know the generalization of the RT formula, but also know the equation of motion that it satisfies [36]. The Lagrangian  for the theories we study is In the second line we have reexpressed the couplings in a way that simplifies some of the expressions below. Before we discuss the calculation, we should mention a limitation of our analysis. Higher derivative theories require completion by massive states such as strings [37]. Because the butterfly velocity comes from a high energy scattering problem, it may receive corrections due to Regge behavior of strings, in addition to the effective field theory higher derivative corrections. This did not happen at the order in α studied in [38], but it may happen at the next order. We have not attempted to check this, but we caution the reader that Regge corrections may disturb the agreement we find between v B andṽ B below.

B.1 Higher derivative corrections to v B from shock waves
The general static planar black hole metric can be written as where we put the horizon at z = 1, hence f (1) = 0. h(z) is a positive function, and we allow arbitrary matter to support the geometry. The shock wave profile can be determined in higher curvature gravity in a straightforward manner following the methods developed in [39], and recently concisely summarized in [40]. Let us write the equations of motion corresponding to (B.1) as: where we treat the negative cosmological constant as matter, and H (i) µν can be found e.g. in [41]. The black hole (B.2) solves this equation of motion with matterT µν . We now add a shockwave on top of this background. It will be convenient to work in Kruskal coordinates and look for the solution: where the contributions fromT in δT have to be included because we are using discontinuous coordinates [39]. Plugging this Ansatz into (B.3) we get an equation for the shock wave profile h(x). In terms of the couplings defined in (B.1), the equation of motion is where the coefficients C i (Λ 1 , Λ 2 ) depend on the black hole metric, and we used that B(0) = 1. 14 Note that the Gauss-Bonnet coupling constant Λ GB doesn't appear in the equation. However, v B still receives corrections in Gauss-Bonnet gravity, all of which are encoded in the change of the black hole solution. 15 Using (5.6), the speed v B can be read off from the tradeoff between the exponential growth of scattering of the blueshifted particles e 2π β t and the decay of the shockwave profile e −µ|x| at large x [11,12]. The decay constant µ can be determined by plugging into the fourth order equation (B.5). We get a second order equation for µ 2 : Their explicit form is given by (B.6) 15 E.g. in a Schwarzshild black hole v B is modified to [12] v B (Λ GB ) = N # d 2(d − 1) , This equation has two solutions for µ 2 . We take the branch which gives back the Einstein result µ 2 = (d−1) 2 f 1 , when we take Λ i → 0. (The other branch diverges as we take Λ 2 → 0.) For comparison with the entanglement wedge computation ofṽ B it will be useful to convert from the Kruskal coordinates expression (B.6) to the coordinates used in (B.2): (B.9)

B.2 Higher derivative corrections to the entanglement wedge
We now determine the entanglement wedge in higher derivative gravity. The static RT surface for a sphere is given by the function z(r) = 1− s(r) 2 , where is small corresponding to close approach to the horion. 16 In [36] a generalized area formula for proposed for higher derivative theories (see also [42]), but the equation of motion satisfied by the RT surface has up to now only been shown to follow from this area function in four derivative and Lovelock theories (see [43] for further discussion on the equation of motion). We will concentrate on four derivative gravity as in (B.1), where the generalized area is: We construct the geometric quantities entering this functional below. The Euler-Lagrange equation from (B.10) determines the RT surface. Let us take a spherically symmetric surface z(r). Because it is a codimension-2 surface, it has two orthogonal normal vectors: From these we can construct the quantities needed for (B.10) P µν = g µν − η ab n (a) µ n (b) ν , K aµν = 1 2 P α µ P β ν L n (a) g αβ , K a = g µν K aµν , R ab = n (a)µ n (b)ν n (c)ρ n (d)σ R µνρσ , R abcd = n (a)µ n (b)ν R µν , (B.12) 16 In section 5 Rindler coordinates were used for the same problem, here we found it more convenient to work in z coordinates. Introducing the variable s(r) instead of its square simplifies our equations.
where a, b indices are contracted with η ab and the vectors are written in polar coordinates. 17 Plugging in everything into the area functional (B.10) the equation of motion is obtained by varying. We get a complicated equation of motion, then we zoom in onto the near horizon region by taking z(r) = 1 − s(r) 2 and expanding for small . To leading order in , we get 0 = (1 + C 1 (Λ 1 , Λ 2 )) s (r) − f 1 d − 1 2 + C 2 (Λ 1 , Λ 2 ) s(r) + Λ 2 s(r) s (r) 2 + 2s(r) s (r) s (3) (r) − 2s (r) 2 s (r) s(r) 2 + O 1 r , where we have used the change of basis in the coupling constants (B.1), and there are many terms hidden in O 1 r , including some that depend on Λ GB . These terms however will not play a role in the following, as we only want to determineμ from the large r asymptotic behavior of the solution: s(r) ∼ eμ r r # . (B.14) From (B.13) it follows that we get the same equation forμ 2 as we got for µ 2 , (B.8), and we have to take the same root forμ 2 that we took for µ 2 , as explained below (B.8). Thus, using that v B = f 1 2µ andṽ B = f 1 2μ , we conclude that v B =ṽ B in four derivative gravity to all orders in the higher curvature couplings. We note that the extrinsic curvature terms in (B.10) played an important role in the above computation.

C Bounds and the quasiparticle model
We can also apply the bounds discussed in section 2 to the quasiparticle model [1,2,6,3] of integrable field theories. In these theories there is no chaos, hence no v B cone, and v LC = c. For such systems, we do not expect the time dependence of the entropy to approximately saturate the bounds.
Let us first consider strip geometries of width 2R for d > 2. For early times t < R:

S[A(t)]
= v E s th · area(A) t , (C.1) but for t > R the result (computed in [6]) significantly deviates from the combined bound, in particular we find that t S = ∞. This is in stark contrast with the holographic results, where (C.1) holds for all times until saturation, and t S = R/v E . 18 17 When doing computations one has to be careful about how to compute the extrinsic curvature, as one needs to have a definition of n (a) away from the surface. 18 Another notable difference is that v (free) E < v (holographic) E , one of the main findings of [6].
However, somewhat coincidentally, when we consider spherical geometries the quasiparticle model comes closer to saturating the bounds. The expressions for S[A(t)] in the quasiparticle model with EPR pattern of entanglement are given in [6]. In figure 10 we compare the results of the quasiparticle computation in d = 3, 5 to the combined bound. The curves are actually somewhat close. For example, the saturation time t S = R saturates the first bound in (2.7). 3) using the quasiparticle value of v E = 2 π , 4 3π and with v LC set to c, and red/solid is the quasiparticle result [6].
We conclude that the quasiparticle model obeys bounds discussed in section 2. While the overall shape of the curves on figures 3 and 10 are similar, the finer details of entropy spread in chaotic systems is different from what one gets in the quasiparticle model. For generic shapes, we do not expect the quasiparticle model to approximately saturate the bounds from section 2.