---
bibliography:
- 'sykgraphbib.bib'
---

$$\begin{tikzpicture}
\draw (\textwidth, 0) node[text width = \textwidth, right] {\color{white} easter egg};
\end{tikzpicture}$$

$$\begin{tikzpicture}
\draw (0.5\textwidth, -3) node[text width = \textwidth] {\huge  \textsf{\textbf{Fast scrambling on sparse graphs}} };
\end{tikzpicture}$$

$$\begin{tikzpicture}
\draw (0.5\textwidth, 0.1) node[text width=\textwidth] {\large \color{black} \textsf{Gregory Bentsen,}$^{\color{{IndianRed3}} \mathsf{a}}$ \textsf{Yingfei Gu}$^{\color{{IndianRed3}} \mathsf{b}}$  \textsf{and Andrew Lucas}$^{\color{{IndianRed3}} \mathsf{a}}$};
\draw (0.5\textwidth, -0.5) node[text width=\textwidth] {$^{\color{{IndianRed3}} \mathsf{a}}$ \small{\textsf{Department of Physics, Stanford University, Stanford, CA 94305, USA}}};
\draw (0.5\textwidth, -0.6) node[text width=\textwidth, below] {$^{\color{{IndianRed3}} \mathsf{b}}$ {\small \textsf{Department of Physics, Harvard University, Cambridge, MA 02138 USA}}};
\end{tikzpicture}$$

$$\begin{tikzpicture}
\draw (0, -13.15) node[right, text width=0.5\paperwidth] {\texttt{ajlucas@stanford.edu}};
\draw (\textwidth, -13.1) node[left] {\textsf{April 2, 2019}};
\end{tikzpicture}$$

$$\begin{tikzpicture}
\draw[very thick, color={IndianRed3}] (0.0\textwidth, -5.75) -- (0.99\textwidth, -5.75);
\draw (0.12\textwidth, -6.25) node[left] {\color{{IndianRed3}}  \textsf{\textbf{Abstract:}}};
\draw (0.53\textwidth, -6) node[below, text width=0.8\textwidth, text justified] {\small Given a quantum many-body system with few-body interactions, how rapidly can quantum information be hidden during time evolution?   The fast scrambling conjecture is that the time to thoroughly mix information among $N$ degrees of freedom grows at least logarithmically in $N$.   We derive this inequality for generic quantum systems at infinite temperature,  bounding the scrambling time by a finite decay time of local quantum correlations at late times.   Using Lieb-Robinson bounds, generalized Sachdev-Ye-Kitaev models, and random unitary circuits, we propose that a logarithmic scrambling time  can be achieved in most quantum systems with sparse connectivity.   These models also elucidate how quantum chaos is not universally related to scrambling:  we construct random few-body circuits with  infinite Lyapunov exponent but logarithmic scrambling time.    We discuss analogies between quantum models on graphs and quantum black holes, and suggest methods to experimentally study scrambling with as many as 100 sparsely-connected quantum degrees of freedom.     };
\end{tikzpicture}$$

$$\begin{tikzpicture}
\draw[very thick, color={IndianRed3}] (0.0\textwidth, -5.75) -- (0.99\textwidth, -5.75);
\end{tikzpicture}$$

Introduction
============

Despite the unitarity of time evolution in quantum mechanics, it is still possible for an isolated many-body quantum system to thermalize and dissipate.  The origin of thermalization in such systems is the spreading of quantum information and entanglement.  Information which was once stored in small regions of space, or in simple operators, can only be accessed by measuring a finite fraction of all $N\gg 1$ degrees of freedom at late times.  Our inability to perform such complicated measurements is why these closed quantum systems still appear thermal. 

Understanding the process of thermalization in quantum systems is of great interest for a number of reasons.  Firstly, a quantitative understanding of quantum thermalization, together with the emergence of hydrodynamics, remains an incredibly hard problem [@polkovnikov] with experimental applications in nuclear, condensed matter and atomic physics.  Secondly, as quantum information has been hidden – but not destroyed – under unitary time evolution, constraints on thermalization also imply bounds on quantum information processing and protection [@hayden07]. 

A natural question is whether thermalization can occur arbitrarily fast.  Remarkably, evidence from black hole physics has led to the conjecture that fundamental bounds on the “scrambling” of quantum information do exist: the time $t_{\mathrm{s}}$ before quantum information can be thoroughly lost to a ‘local’ observer obeys [@susskind08] $$t_{\mathrm{s}} \gtrsim  \frac{\log N}{\lambda_{\mathrm{s}}},  \;\;\;\; (N\rightarrow \infty) \label{eq:introbound}$$ with $\lambda_{\mathrm{s}}\lesssim N^0$. This conjecture comes from a thought experiment where black holes appear to clone quantum states (in violation of a theorem [@zurek; @dieks]) unless $t_{\mathrm{s}}$ obeys (\[eq:introbound\]).  Three immediate questions arise.  Can this conjecture be understood on general quantum mechanical grounds, without appealing to black holes?  If so, what is $\lambda_{\mathrm{s}}$?  In which systems can the bound (\[eq:introbound\]) be saturated?  The purpose of this paper is to explore and resolve such questions in a broad class of quantum systems. 

We first clarify the conjecture.  For simplicity, consider a Hilbert space $\mathcal{H}$ which is a tensor product of $N$ two-dimensional Hilbert spaces $\mathcal{H}_i$. We refer to each $\mathcal{H}_i$ as a single *degree of freedom* (DOF).  We shall be concerned with the class of Hamiltonians $H$ acting on $\mathcal{H}$ that are $k$-local in the computer science sense: namely, they involve products of at most $k$ Pauli matrices $\sigma^\alpha_i$: $$\begin{aligned}
H &= \sum_{i=1}^N \sum_{\alpha=1}^3  h_i^\alpha \sigma_i^\alpha + \sum_{i,j=1}^N \sum_{\alpha,\beta=1}^3  h_{ij}^{\alpha \beta} \sigma_i^\alpha \sigma_j^\beta  +  \cdots \notag \\
&+ \sum_{i_1\cdots i_k, \alpha_1\cdots \alpha_k} h_{i_1\cdots i_k}^{\alpha_1 \cdots \alpha_k} \sigma_{i_1}^{\alpha_1}\cdots \sigma_{i_k}^{\alpha_k}.   \label{eq:introH}\end{aligned}$$ We will always take $k$ to be finite in the thermodynamic limit $N\rightarrow  \infty$. The coefficients above are normalized so that the Hamiltonian is extensive.  We define the scrambling time as the time after which two initially unentangled halves of a system become nearly maximally entangled.  Let $A \subset \lbrace 1,\ldots, N \rbrace$ be a subset of $N^\prime=pN$ of the degrees of freedom, with $0<p{\leqslant}\frac{1}{2}$, such that $\mathcal{H} = \mathcal{H}_A \otimes \mathcal{H}_{A^c}$. Given a time-dependent pure state $|\Psi(t)\rangle$, we define the reduced density matrix of $A$: $$\rho_A(t) = \underset{A^{\mathrm{c}}}{\mathrm{tr}} |\Psi(t)\rangle\langle \Psi(t)|,$$ and the von Neumann entanglement entropy $$S_A(t) =  -\underset{A}{\mathrm{tr}}\left[\rho_A(t) \log \rho_A(t)\right].$$ If the initial state $|\Psi(0)\rangle  = |\Psi_A\rangle \otimes |\Psi_{A^{\mathrm{c}}}\rangle$ has no entanglement between $A$ and its complement, then we define the scrambling time as the smallest time when [@susskind08] $$S_A(t_{\mathrm{s}}) {\geqslant}S_A^{\mathrm{max}} - a, \label{eq:SAmax}$$ regardless of the choice of initial unentangled state and for any bipartition $A,A^{\mathrm{c}}$. Here $S_A^{\mathrm{max}} = N^\prime \log 2$, and $a \propto N^0$ is a fixed constant offset as $N\rightarrow\infty$. Since $\log N\approx \log N^\prime$, we will see that the choice of $a$ or $p$ should not affect (\[eq:introbound\]) at leading order in $N$. Extending these definitions to quantum systems where $\mathcal{H}_i$ are higher-dimensional is straightforward.  Another generalization of the fast scrambling conjecture in more quantum information theoretic terms was given in [@winter]. 

Fast Scrambling from Local Correlation Functions
================================================

We are now ready to answer the first of the questions posed previously: can the fast scrambling conjecture be proved for arbitrary quantum systems using purely quantum mechanical arguments?  We sketch a proof here by demonstrating that scrambling is ultimately limited by the slow growth of entanglement at late times.  Specifically, we shall find that $t_{\mathrm{s}} \gtrsim \log N$ because a fixed, generic $k$-local Hamiltonian cannot readily entangle two subsystems that are already nearly maximally entangled.  This slowdown of entanglement growth at late times can be related to the decay of two-point correlators, which generically decay *exponentially*.  While entanglement may grow rapidly at early times, this exponential time-dependence implies that the approach to maximal entanglement at long times is not arbitrarily fast. 

We make this intuition precise by considering an initially unentangled set of $N$ qubits evolving under an arbitrary $k$-local Hamiltonian, and observe the approach to scrambling from the point of view of individual qubits $i$. (These arguments are easily generalized to DOF with Hilbert spaces of arbitrary dimension.)  Close to the scrambling time $t_s$ a typical qubit $i$ will be highly entangled with the rest of the system, so that its reduced density matrix $\rho_i$ is nearly maximally mixed: $$\rho_i(t) = \frac{1}{2} I_i + \frac{1}{2} \sum_{\alpha=1}^3 \epsilon_i^{\alpha}(t) \sigma^{\alpha}_i
    \label{eq:SingleSiteRho}$$ where $\epsilon_i^{\alpha}(t)$ are parameters encoding the approach of the density matrix toward the maximally mixed state $\rho_i(t) \rightarrow \frac{1}{2} I_i$. Since the qubits are initially unentangled, $\rho_i(0)$ is pure: as a consequence, $\max_{\alpha} (\epsilon_i^\alpha(0)) {\geqslant}\frac{1}{\sqrt{3}}$. The fast scrambling conjecture follows from two assertions: firstly, that the functions $\epsilon_i^{\alpha}(t)$ must decay to a small value of order $\mathrm{O}(1/\sqrt{N})$ at the scrambling time $t=t_{\mathrm{s}}$; and secondly, that the decay of $\epsilon_i^\alpha(t)$ can occur no faster than exponentially quickly: $\epsilon_i^\alpha(t) \gtrsim \mathrm{e}^{-\lambda_2 t} $ as $t\rightarrow \infty$, for a finite rate $\lambda_2 \propto N^0$. Together, these statements imply the fast scrambling conjecture: it is impossible to scramble at time $t \ll \lambda_2^{-1}\log N $. 

The first statement is proved by considering the formal definition of scrambling (\[eq:SAmax\]), which requires that the entanglement entropy $S_A$ of any subsystem $A$ be nearly maximal.  By the property of sub-additivity, the entanglement entropy of $A$ can be no larger than the sum of the entanglement entropies of its constituent qubits: $$S_A(t) {\leqslant}\sum_{i \in A} S_i(t)
    \label{eq:Subadditivity}$$ where $S_i$ denotes the entanglement entropy of qubit $i$ with the other $N-1$ qubits.  Intuitively, (\[eq:Subadditivity\]) holds because the right hand side measures entanglement between different qubits in $A$ that will not be counted in $S_A$. Using (\[eq:SingleSiteRho\]), we can then quantify the growth of single-qubit entanglement $S_i(t)$ in terms of the decaying functions $\epsilon_i^{\alpha}(t)$. Taylor expanding $S_i(t) = -\mathrm{tr}[\rho_i(t)\log \rho_i(t)]$, $$S_i(t) = \log 2 - \frac{1}{4} \sum_{\alpha=1}^3 \left( \epsilon_i^{\alpha}(t) \right)^2 + \mathrm{O}\left( \epsilon^3 \right). 
    \label{eq:EntropyGrowth}$$ This expansion is sensible for $\epsilon_i^\alpha(t) \ll 1$, at a time when the initially unentangled qubits become significantly entangled.  Combining (\[eq:Subadditivity\]) and (\[eq:EntropyGrowth\]), we obtain: $$S_A(t) < N' \log 2 - \frac{1}{4} \sum_{i \in A} \sum_{\alpha=1}^3 \left( \epsilon_i^{\alpha}(t) \right)^2 + \mathrm{O}\left( \epsilon^3 \right).   \label{eq:SAt}$$ In order for $S_A(t)$ to be nearly maximal as defined in (\[eq:SAmax\]), we require that $\sum_{i,\alpha} \left( \epsilon_i^{\alpha}(t) \right)^2 < 4 a$ or that $\epsilon_i^{\alpha}(t) {\leqslant}\sqrt{\frac{8a}{N}}$ for at least half of the qubits $i$. When $N\rightarrow \infty$, $\mathrm{O}(\epsilon^3)$ corrections to (\[eq:SAt\]) are negligible. 

It remains to demonstrate that the functions $\epsilon_i^{\alpha}(t)$ cannot decay arbitrarily quickly, or equivalently, that entanglement cannot build up in the system arbitrarily quickly.  This is the main physical content of the fast scrambling conjecture.  To show this, we relate the decay of the functions $\epsilon_i^{\alpha}(t)$ to the decay of two-point correlation functions, and argue on general grounds that such correlators cannot decay arbitrarily quickly.  We note that other recent work [@magan; @stanford1802; @shenker1803; @lensky] has also discussed contributions to scrambling arising from the decay of two-point functions. 

Using (\[eq:SingleSiteRho\]), we may write: $$\epsilon_i^{\alpha}(t)^2 = \left( \mathrm{tr} \left[ \rho_i(t) \sigma_i^{\alpha} \right] \right)^2 = \left( \mathrm{tr} \left[ \rho \ \sigma_i^{\alpha}(t) \right] \right)^2
    \label{eq:EpsToCorrs}$$ where $\rho$ is the initial (separable) state of the entire many-qubit system and $\mathrm{tr}$ denotes a trace over the entire system.  We expect that single-site correlation functions for the initial state $\rho$ can decay no faster than they decay for the initial state $\rho_i(0) \otimes \rho^{\infty}_{-i}$, where $\rho^{\infty}_{-i} = I_{-i} / 2^{N-1}$ is the maximally-mixed state of all qubits other than qubit $i$. That is, for a constant $C^\prime\propto N^0$, $$\left( \mathrm{tr} \left[ \rho \ \sigma_i^{\alpha}(t) \right] \right)^2 {\geqslant}C^\prime \left( \mathrm{tr} \left[ \left( \rho_i(0) \otimes \rho^{\infty}_{-i} \right) \sigma_i^{\alpha}(t) \right] \right)^2. 
    \label{eq:DecayCompareToInfTemp}$$ The intuition behind this inequality is that the state $\rho^{\infty}_{-i}$ serves as an infinite-temperature bath that relaxes qubit $i$ faster than any other initial state.  A more formal perspective is found in *Supporting Information*.  Combining (\[eq:EpsToCorrs\]) and assuming (\[eq:DecayCompareToInfTemp\]) holds exactly as an inequality: $$\epsilon_i^{\alpha}(t)^2 {\geqslant}\left( \frac{1}{2^N} \sum_{\beta=1}^3 \mathrm{tr} \left[ \sigma_i^{\alpha}(t) \sigma_i^{\beta} \right] \epsilon_i^{\beta}(0) \right)^2 {\geqslant}C \mathrm{e}^{-2\lambda_2 t}  \label{eq:2lambda2}$$ where $\lambda_2$ is the decay rate of infinite-temperature two-point correlators, and $C$ is an O(1) constant.  Crucially, the decay rates $\lambda_2$ will be *finite* for generic qubit $i$, with generic Hamiltonian.  In *Supporting Information*, we find that $\lambda_2 \propto N^0$ whenever the Hamiltonian $H$ is extensive and $k$-local (with $k$ finite), and when the density of eigenvalues of $H$ is a smooth function in the large-$N$ limit. 

Combining (\[eq:SAt\]) and (\[eq:2lambda2\]) with the definition of the scrambling time (\[eq:SAmax\]), we find that $2\lambda_2 {\geqslant}\lambda_{\mathrm{s}}$. Entanglement grows slowly at late times whenever local two-point functions decay exponentially (which generically does occur).  Hence we derive the logarithmic lower bound (\[eq:introbound\]) on scrambling.  The origin of the logarithmic divergence in the scrambling time is the slow saturation to a maximally entangled state at late times. 

Out-of-Time-Ordered Correlators and Operator Growth
===================================================

We now turn to our final question: when might (\[eq:introbound\]) be saturated?  The exponential decay of two-point functions is necessary but not sufficient for fast scrambling: for example, a system made of two decoupled fast scramblers is not itself a fast scrambler (information cannot spread between decoupled systems).  To answer this question, therefore, we now turn to another probe of fast scrambling: the growth of operators [@lensky; @lashkari; @hosur].  In particular, we bound the growth of operators, and $t_{\mathrm{s}}$, by how the Hamiltonian connects the $N$ degrees of freedom. 

Recall that in the Heisenberg picture, an operator $\mathcal{O}$ evolves in time as $\mathcal{O}(t) = \mathrm{e}^{\mathrm{i}Ht} \mathcal{O} \mathrm{e}^{-\mathrm{i}Ht}$. Suppose that at time $t=0$, $\mathcal{O}$ acts only on qubits in $A$: we may store quantum information by preparing the system in the mixed state $\rho(0) = |\Psi_{\mathcal{O}}\rangle\langle \Psi_{\mathcal{O}}| \otimes 2^{-|A^{\mathrm{c}}|} I_{A^{\mathrm{c}}}$, where $|\Psi_{\mathcal{O}}\rangle$ is an eigenvector of $\mathcal{O}$. Whereas initially $S_A[\rho(0)] = 0$, in a scrambling quantum system $S_A[\rho(t)]$ becomes large at late times: equivalently, almost every eigenvector of $\mathcal{O}(t)$ must be highly entangled between $A$ and $A^{\mathrm{c}}$, even though $\mathcal{O}$ acts only on $A$.  This is only possible if $\mathcal{O}(t)$ itself is a complicated operator.  We may parameterize the growth of the operator $\mathcal{O}(t)$ by writing $$\mathcal{O}(t) = \sum_{R\subseteq\lbrace 1,\ldots, N\rbrace}   a_R(t) {\widetilde}{\mathcal{O}}_R(t),  \label{eq:Ovt}$$ where ${\widetilde}{\mathcal{O}}_R(t)$ is an operator that acts non-trivially on qubit $i$ if and only if $i\in R$, with $\lVert {\widetilde}{\mathcal{O}}_R \rVert = 1$, where $\lVert \mathcal{O}\rVert^2 = 2^{-N} \mathrm{tr}(\mathcal{O}^\dagger\mathcal{O})$ is the operator norm.  Note that the real coefficients $a_R(t)$ obey $\sum_R a_R(t)^2 = 1$. Whereas at $t=0$, $a_R(t)\ne 0$ only when $R\subseteq A$, at the scrambling time almost all weight $a_R(t)$ is on large operators where $R$ contains $\mathrm{O}(N)$ qubits in both $A,A^{\mathrm{c}}$.

One way to effectively describe the growth of operators without listing every $a_R(t)$ is to ask what fraction of the operator $\mathcal{O}(t)$ acts non-trivially on qubits $j$ in $A^{\mathrm{c}}$. Expanding the operators $\mathcal{O}(t),{\widetilde}{\mathcal{O}}_R(t)$ over the complete basis of tensor products of Pauli matrices, we find that the overlap with qubit $j$ is [@lucas1809] $$\frac{1}{8}\sum_{\alpha=1}^3\lVert [\mathcal{O}(t),\sigma^\alpha_j] \rVert^2 = \sum_{R: j\in R} a_R(t)^2.   \label{eq:aRbound}$$ Using (\[eq:aRbound\]), we can easily find an operator $\sigma^\alpha_j$ for which $$\lVert [\mathcal{O}(t),\sigma^\alpha_j] \rVert^2 {\geqslant}\frac{8}{3} \sum_{R: j\in R} a_R(t)^2. 
\label{eq:aRbound2}$$ Since operator growth is bounded from above by (\[eq:aRbound2\]), if we can place an upper bound on the time $t_*$ at which $\lVert [\mathcal{O}(t_*),\sigma_j^{\alpha}] \rVert$ becomes O(1), then we bound the scrambling time.  We achieve this by generalizing [@lashkari; @sims; @hastings] the Lieb-Robinson theorem [@liebrobinson] to more general Hamiltonians such as (\[eq:introH\]).  Historically, the Lieb-Robinson theorem was used to show that commutator norms in lattice models can only be large inside of an emergent “lightcone.”  Here, we will strengthen the Lieb-Robinson theorem and show that the growth of operators is intimately tied to the structure of the interaction graph of the Hamiltonian. 

For simplicity, we will focus on 2-local Hamiltonians, leaving a more general discussion to *Supporting Information*.  To any 2-local Hamiltonian $$H = \sum_{(u, v)\in E} H_{uv},    \label{eq:HvHuv}$$ we may associate a discrete, undirected graph $G = (V,E)$, where DOF live on the vertices $v \in V$ and interact pairwise via couplings $H_{uv}$ if and only if the pair $u,v$ is connected by an edge $e = (u,v) \in E$. The connectivity of the graph $G$ is described by the *adjacency matrix* $$A_{uv} = \left\lbrace \begin{array}{ll} 1 &\  (u,v) \in E \\  0 &\ (u,v)\notin E\end{array}\right.$$ and the *degree matrix* $D_{uv} = k_v {\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{uv}$, where the *degree* $k_v = \sum_u A_{uv}$ of a given vertex $v$ counts the number of vertices it is connected to.  Denoting $$\frac{c}{K} = \max_{u,v} \left( \lVert H_{uv} \rVert \right),  \label{eq:maincdef}$$ where $K=\frac{1}{N}\sum_v k_v$ is the mean degree, and assuming $c\propto N^0$ by extensivity, our generalized Lieb-Robinson bound is $$\frac{\lVert [\mathcal{O}_u,  \mathcal{O}_v(t)] \rVert}{2\lVert \mathcal{O}_u \rVert\lVert \mathcal{O}_v \rVert} {\leqslant}\exp\left[\frac{2 c|t|}{K}(D+A)\right]_{uv} \label{eq:LR1st}$$ where $\mathcal{O}_v$ is any operator with support only on vertex $v$. This result is derived by using the triangle inequality: $\lVert [\mathcal{O}_u, \mathcal{O}_v(t)] \rVert {\leqslant}\lVert [\mathcal{O}_u, \mathcal{O}_v] \rVert + t \lVert [\mathcal{O}_u, [H,\mathcal{O}_v]] \rVert + \mathrm{O}(t^2)$. Commutators between operators acting on disjoint sets cancel, so if $u\ne v$, $\lVert [\mathcal{O}_u, [H,\mathcal{O}_v]] \rVert = \lVert [\mathcal{O}_u, [H_{uv},\mathcal{O}_v]] \rVert $.  Furthermore, $\lVert [\mathcal{O}_u, [H_{uv},\mathcal{O}_v]] \rVert {\leqslant}4 \lVert \mathcal{O}_u \rVert \lVert H_{uv} \rVert \lVert \mathcal{O}_v \rVert$. These bounds can be repeated at every order in $t$ and resummed to give (\[eq:LR1st\]): see *Supporting Information*.  Crucially, since the result (\[eq:LR1st\]) bounds the operator norm $\lVert \cdot \rVert$, it bounds operator growth at *infinite temperature* since $\lVert \mathcal{O}\rVert^2 = \mathrm{tr}(\rho^{\infty}\mathcal{O}^\dagger \mathcal{O})$ where $\rho^{\infty} = 2^{-N} I_N$ is the infinite-temperature thermal state of all $N$ qubits. 

It is instructive to average (\[eq:LR1st\]) over all choices of $u$ and $v$: $$\sum_{u,v} \frac{1}{N^2} \frac{\lVert [\mathcal{O}_u,  \mathcal{O}_v(t)] \rVert}{2\lVert \mathcal{O}_u \rVert\lVert \mathcal{O}_v \rVert}  {\leqslant}\frac{1}{N} \exp\left[\frac{4ck_{\mathrm{max}}}{K} |t| \right] \label{eq:expgrowth}$$ where $k_{\mathrm{max}}$ is the maximal degree in the graph.  Whenever $k_{\mathrm{max}}/K$ is finite, we see that operator growth is constrained to be exponentially fast.  The form of the bound (\[eq:expgrowth\]) is particularly useful because it naturally bounds the growth of *out-of-time-ordered correlators* (OTOCs) such as $\mathrm{tr}(\rho [\mathcal{O}_u,  \mathcal{O}_v(t)]^2) \sim \frac{1}{N} \mathrm{e}^{\lambda_{\mathrm{L}}t}$ which measure the growth of many-body chaos with Lyapunov exponent $\lambda_L$ [@shenker13; @douglasweak; @patel; @debanjan; @erez; @scopelliti].  Our Lieb-Robinson bound (\[eq:expgrowth\]) establishes that the growth rate of infinite-temperature OTOCs can be no faster than $$\lambda_{\mathrm{L}} {\leqslant}\frac{8ck_{\mathrm{max}}}{K},$$ thereby placing a bound on quantum many-body chaos at infinite temperature.  Moreover, since (\[eq:aRbound2\]) implies that scrambling may only occur once $\lVert [\mathcal{O}_u,  \mathcal{O}_v(t)]\rVert$ has grown to be O(1), we conclude from (\[eq:expgrowth\]) that $\lambda_{\mathrm{s}} {\leqslant}\lambda_{\mathrm{L}}$. On any graph where $k_{\mathrm{max}}/K$ is finite, we therefore have an alternative proof of the fast scrambling conjecture. 

Our Lieb-Robinson bounds also constrain which graphs are capable of supporting fast scrambling.  Further manipulations to (\[eq:LR1st\]) (see *Supporting Information*) lead to $$\exp\left[\frac{2c|t|}{K}(D+A)\right]_{uv} < \exp\left[\frac{4\mathrm{e} ck_{\mathrm{max}}}{K} |t|  - d_{uv}\right],  \label{eq:LRdistmain}$$ where $\mathrm{e}\approx2.718$ and $d_{uv}$ is the minimal path length between vertices $u$ and $v$. This is the classic Lieb-Robinson bound [@liebrobinson].  Operators spread at most ballistically through the graph, and the timescale required for an arbitrary operator to spread throughout the entire system is limited by the path length between the two most distant sites $\max_{u,v}(d_{uv})$, or the *graph diameter*.  We immediately conclude that only graphs with diameter $\max_{u,v}(d_{uv}) \lesssim \log N$ can support fast scrambling [@lashkari; @hastings].  Fast scrambling is thus impossible for spin models on regular lattices in $D$ spatial dimensions [@prosen], where the diameter $\propto N^{1/D}$. On the other hand, fast scrambling is expected to be rather generic on almost every connected sparse random graph, whose diameter is $\mathrm{O}(\log N)$ for any $K>2$ [@bollobas].  How “random” must a sparse graph of fixed degree $k_v = K$ be to host a fast scrambler?  It is likely that any finite amount of non-locality is sufficient.  For example, consider the small world network [@watts] depicted in Figure 1d, which is constructed by starting with a one-dimensional cyclic graph where each of the $N$ vertices connects to its nearest $K/2$ neighbors on each side, and randomly re-wiring a fraction $p \ll  1/K$ of the edges to connect two randomly chosen vertices.  The typical distance between two vertices in the small world graph is [@newman] $$\frac{1}{N^2}\sum_{u,v=1}^N d_{uv} \ \approx \  \frac{\log N}{4Kp}.$$ Using (\[eq:LRdistmain\]), we see that infinitesimal randomness ($p>0$) is sufficient for logarithmic scrambling.  In order to check our assertion that fast scrambling is achievable on generic sparse graphs, we study a maximally chaotic toy model on a random graph by choosing the $H_{uv}$ in (\[eq:HvHuv\]) to be generalized Sachdev-Ye-Kitaev (SYK) models [@sachdevye; @kitaevunpublished; @sachdev15; @maldacena2016remarks].  The SYK model is a model of $M \gg 1$ Majorana fermions interacting via random all-to-all $q$-body interactions.  To construct our toy model, we select a graph $G$ consisting of $L$ vertices and place one SYK model on each vertex $v$. The local Hilbert space on each vertex has dimension $\dim(\mathcal{H}_v) = 2^{M/2}$. We then couple pairs of vertices $u,v$ by introducing random $q$-body interactions involving $q/2$ fermions from each of the two vertices $u,v$ whenever there is an edge $(u,v)$ in $G$. The result is a model of $N = ML/2$ strongly-interacting degrees of freedom that is both maximally chaotic and analytically tractable.  Our main result for these generalized SYK models is the computation of OTOCs; as discussed earlier, the time $t_*$ when OTOCs become large is a lower bound on $t_{\mathrm{s}}$. We find that whenever the graph Laplacian $\Lambda_{uv} = D_{uv}-A_{uv}$ has a single vanishing eigenvalue as $N\rightarrow \infty$: $$t_* = \frac{\log M + \log L}{\lambda_{\mathrm{L}}^{\mathrm{SYK}}} + \mathrm{O}(1),  \label{eq:tstarSYK}$$ where $\lambda_{\mathrm{L}}^{\mathrm{SYK}}$ is the Lyapunov exponent of the fully connected SYK model, which has previously been computed [@maldacena2016remarks].  As $\log M + \log L \approx \log N$, (\[eq:tstarSYK\]) is the same time that it takes for OTOCs to grow large in a fully connected SYK model of $N$ fermions.  Remarkably, (\[eq:tstarSYK\]) holds whenever the graph has the property that for an arbitrary subset $A\subset V$ with $|A|\gg 1$, the “surface area” of the subregion (i.e. the number of vertices in $A$ connected to vertices not in $A$) is proportional to $|A|$ (see *Supporting Information*, and Figure 1).  Generic random graphs have this property, implying that all-to-all connectivity in $G$ can be redundant and is not necessary for fast scrambling. 

Random Quantum Circuits
=======================

Now, let us return to the question of OTOC growth when $k_{\mathrm{max}}/K$ diverges as $N\rightarrow \infty$. In this regime, (\[eq:LR1st\]) and (\[eq:expgrowth\]) do not rule out rapid OTOC growth in a time $t_*$ independent of $N$ (or worse).  A simple example is the Hamiltonian $$H = \sum_{i=1}^{N-1} h_{iN}^{\alpha\beta} \sigma^\alpha_i \sigma^\beta_N. \label{eq:sigmazstar}$$ For this model, (\[eq:LR1st\]) implies that $$\frac{\lVert [\mathcal{O}_1,  \mathcal{O}_2(t)] \rVert}{2\lVert \mathcal{O}_1 \rVert\lVert \mathcal{O}_2 \rVert} {\leqslant}\frac{N-1 + \mathrm{e}^{2cNt} - N \mathrm{e}^{2ct}}{N(N-1)}.    \label{eq:superNLR}$$ In this case, the bound (\[eq:LR1st\]) allows the possibility that OTOCs are large by the time $t_* \sim N^{-1}\log N$.

Is (\[eq:LR1st\]) accurately capturing incredibly rapid operator growth, or is it simply a lousy upper bound on an exponentially growing commutator norm?  We can answer this question explicitly in a solvable quantum dynamical system: a random unitary circuit (RUC) [@hayden07; @fawzi; @nahum16; @nahum; @tibor].  In the RUC, the time evolution operator $\mathrm{e}^{-\mathrm{i}Ht}$ is replaced by a product of $\lceil Nt \rceil$ unitary operations $U$, each of which act on a finite number $k$ of the DOF (to mimic $k$-locality).  The simplifying assumption is that we uniformly average over all allowed local unitary operations, which reduces the resulting dynamics to a classical stochastic process [@nahum16; @nahum; @tibor]. 

For example, consider a 2-local RUC in the limit where the local Hilbert space dimension $d \rightarrow \infty$. (A full discussion of RUC dynamics at finite $d$ is presented in *Supporting Information*.)  A 2-local unitary $U_{uw}$ acting on vertices $u$ and $w$ almost surely acts on tensor product operators as follows, where we use $\mathcal{O}$ to represent any non-identity operator:

\[eq:RUCdynamics\]$$\begin{aligned}
U_{uw}^\dagger (I_u \otimes I_w \otimes \mathcal{O}_{\mathrm{rest}}) U_{uw} &= I_u \otimes I_w \otimes \mathcal{O}_{\mathrm{rest}}, \\
U_{uw}^\dagger (I_u \otimes \mathcal{O}_w \otimes \mathcal{O}_{\mathrm{rest}}) U_{uw} &= \mathcal{O}_u \otimes \mathcal{O}_w \otimes \mathcal{O}_{\mathrm{rest}}, \\
U_{uw}^\dagger (\mathcal{O}_u \otimes I_w \otimes \mathcal{O}_{\mathrm{rest}}) U_{uw} &= \mathcal{O}_u \otimes \mathcal{O}_w \otimes \mathcal{O}_{\mathrm{rest}}, \\
U_{uw}^\dagger (\mathcal{O}_u \otimes \mathcal{O}_w \otimes \mathcal{O}_{\mathrm{rest}}) U_{uw} &= \mathcal{O}_u \otimes \mathcal{O}_w \otimes \mathcal{O}_{\mathrm{rest}}.\end{aligned}$$

Here $\mathcal{O}_{\mathrm{rest}}$ denotes operators acting on the remaining $N-2$ vertices.  Suppose that the allowed $u$ and $w$ correspond to a uniformly chosen edge from a given graph $G$. From (\[eq:RUCdynamics\]), if an arbitrary operator ${\widetilde}{\mathcal{O}}$ acts either on $u$ or $w$, then $U^\dagger_{uw} {\widetilde}{\mathcal{O}} U_{uw}$ acts on both $u$ and $w$. A series of random unitary operations $U_{uw}$ applied in this manner will spread information from an initial vertex $v$ across the entire system, similar to the way in which an operator $\mathcal{O}_v(t)$ initially localized on vertex $v$ grows under time evolution as in (\[eq:Ovt\]).  This is reminiscent of a spreading ‘infection,’ where we say a vertex $u$ is ‘infected’ if and only if a non-identity operator $\mathcal{O}_u$ is present on $u$. From (\[eq:RUCdynamics\]), infections spread whenever we choose an edge that connects infected and uninfected vertices.  The resulting infection dynamics is equivalent to the well-known SI model [@vespignaniRMP].  Therefore, we can use the theory of infection spreading on general graphs to understand the spread of operators in the RUC. 

Of particular interest is the fact that infections spread faster than exponentially on heterogeneous networks in the SI model [@vespignani].  A heterogeneous network is one where the probability $p(k)$ that vertex $v$ has degree $k_v=k$ has a heavy tail, with the variance of $p(k)$ diverging as $N\rightarrow \infty$. An example is the interaction graph of (\[eq:sigmazstar\]).  An RUC on this graph will be ‘super-chaotic,’ with $\lambda_{\mathrm{L}}=\infty$. To model the commutator in (\[eq:superNLR\]), we start with vertex 1 infected at $t=0$. At time $t_0$ (whose probability distribution is $\mathbb{P}(t_0 > s) = \mathrm{e}^{-s}$) the infection will jump to the central vertex $N$. Starting from the time $t_0$ when $N$ is infected, the probability that vertices $2,\ldots,N-1$ are infected is $1-\mathrm{e}^{t_0-t}$ ($t>t_0$).  At time $t=t_0 + \log 2$, vertex 2 is infected with probability $\frac{1}{2}$; hence $\lVert [\mathcal{O}_1(t), \mathcal{O}_2]\rVert$ is large.  Thus, this RUC exhibits an operator growth time $t_* \propto N^0 \ll \log N$. While our bound (\[eq:superNLR\]) is very weak, the conclusion that operator growth is faster than on a regular graph, where all vertices have the same degree, remains. 

In agreement with our earlier discussions, however, the scrambling time in the RUC is ultimately limited by the slow growth of entanglement entropy at late times.  To demonstrate this, suppose we start in a tensor product state $|\Psi(0)\rangle = \otimes_v |\Psi_v\rangle$. Choose an arbitrary subset $A$ of half of the vertices, with $S_A(t=0)=0$. To generate nearly maximal entanglement $S_A(t) {\geqslant}\frac{N}{2}\log d -1 $, at least one unitary must hit every vertex $v\in A$; otherwise, the reduced density matrix $\rho_A = |\Psi_v \rangle \langle \Psi_v | \otimes \cdots$ and we obtain $S_A {\leqslant}(\frac{N}{2}-1)\log d$. In the RUC, the probability that we have not chosen a vertex $v$ by time $t$ is $\exp[-\frac{k_v}{K}t]$, where $k_v$ is the degree of $v$ and $K$ is the mean degree.  At least half of the vertices $v$ have $k_v {\leqslant}2K$. Therefore, the time $t_A$ it takes to choose every vertex in $A$ obeys $t_A {\geqslant}\frac{1}{2}\log N$ with high probability: at shorter times at least one vertex has not yet been chosen.  Since $t_{\mathrm{s}} {\geqslant}t_A$, the fast scrambling conjecture (\[eq:introbound\]) holds, despite the system exhibiting an infinite Lyapunov exponent. 

Discussion
==========

We have shown that the scrambling time in few-body quantum dynamics is always bounded by the decay rates of local operators.  By contrast, there is no universal relationship between fast scrambling and quantum chaos.  In particular, we emphasize that the rapid growth of OTOCs is not always a reliable indicator of fast scrambling.  For instance, the SYK model is highly chaotic, but generates entanglement extremely slowly at late times [@gu1708].  On the other hand, we found a RUC model which is ‘super-chaotic’ ($\lambda_{\mathrm{L}} = \infty$) at early times, yet still obeys the fast scrambling conjecture.  Furthermore, while black holes are believed to be the fastest scramblers in nature [@susskind08], they are ‘less’ chaotic than this RUC: known Hamiltonians related to black holes have a finite Lyapunov exponent at all temperatures.  One reason for this is that microscopic quantum models of black holes [@bfss] have homogeneous, $k$-local, all-to-all connectivity.  Heterogeneous connectivity is required to obtain $\lambda_{\mathrm{L}}=\infty$, as shown by our generalized Lieb-Robinson bounds. 

That scrambling is instead intimately connected to local decay rates is reminiscent of the historical origins of the fast scrambling conjecture in black hole physics.  The ringdown of quasinormal modes in a time $t_{\mathrm{s}} \propto \log N$ is in large part the motivation for the conjecture [@hayden07; @susskind08].  Furthermore, the scrambling time of black holes appears as fast as possible without measurable violations of quantum no-cloning theorems [@hayden07].  Rather than cloning information, black holes rapidly generate entanglement between degrees of freedom.  Our arguments show that the rate of this entanglement growth is limited by the finiteness of decay rates in general quantum systems.  In the context of the AdS/CFT correspondence [@lucasbook], these decay rates are those of black hole quasinormal modes.  This makes the original justification of fast scrambling precise. 

There is another remarkable analogy between black holes and quantum dynamics on random graphs.  The area of a black hole is proportional to the thermal entropy [@hawking], which is in turn proportional to the total number of degrees of freedom $N$. Generic random graphs have a similar property, where a significant fraction of vertices in any connected region $A \subset V$ are located at the *boundary* of $A$ – they connect to the ‘outside.’  This is in contrast to, say, a regular lattice, in which most of the vertices of any particular connected region $A$ reside deep within the volume of $A$ and do not connect to vertices outside $A$.  See *Supporting Information* for a precise statement, and Figure \[fig:graphsfig\] for a pictorial demonstration of this result; see also [@magan2].  As a specific example, on any graph with this property “$\text{area}\propto N$,” our generalized SYK model is just as chaotic as the fully connected model, where the fully connected model is itself a crude model of a black hole [@kitaevunpublished; @maldacena2016remarks]. 

![ A selection of different graphs of relevance in this paper: (a) a regular lattice; (b) a complete graph; (c) a sparse random graph which is locally treelike; (d) a small world graph. ( a) cannot host fast scramblers due to the large graph diameter, but (b)-(d) can host fast scramblers.  A subset $A$ of vertices in each graph above is denoted in red; vertices in $A^{\mathrm{c}}$ are denoted in blue.  Light red vertices denote vertices away from the ‘area’ of $A$, and dark red vertices are at the boundary of $A$ and $A^{\mathrm{c}}$. In the thermodynamic limit, on graphs (b)-(d), a finite fraction of all edges starting in $A$ end in $A^{\mathrm{c}}$: thus area is proportional to the number of DOF inside the region.[]{data-label="fig:graphsfig"}](graphsnewfig_v2.pdf){width="6in"}

Outlook
=======

The most important extension of our results is to finite temperature quantum dynamics.  We conjecture that at all temperatures and for all initially unentangled states, the scrambling time is limited by the decay of thermal two point functions.  A proof likely requires more sophisticated methods than what we have developed here: see [@lensky] for some preliminary directions.  More interestingly, our explicit construction of random circuits with $\lambda_{\mathrm{L}}=\infty$ is surprising given that, under mild assumptions, Lyapunov exponents generally are bounded from above by $\lambda_{\mathrm{L}} {\leqslant}2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}k_{\mathrm{B}}T/\hbar$ at finite temperature $T$ [@stanfordbound], where $\hbar = h/2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}$ and $k_{\mathrm{B}}$ and $h$ are Boltzmann’s and Planck’s constants, respectively.  Our work does not contradict [@stanfordbound], however: the RUC is not generated by a fixed Hamiltonian, and should be considered as a model of infinite temperature dynamics where the bound of [@stanfordbound] is trivial.  However, a necessary condition to obtain $\lambda_{\mathrm{L}} {\leqslant}2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}k_{\mathrm{B}}T/\hbar$ is that OTOCs stay sufficiently small ($\mathrm{O}(\frac{1}{N})$) at early times.  This condition fails on heterogeneous graphs. 

Our discussion of scrambling on sparse graphs is also relevant for experimental efforts to realize highly chaotic quantum systems.  Simulations of Hamiltonians on sparse graphs can leverage, for instance, the non-local interactions naturally afforded by cavity quantum electrodynamics (QED) setups.  Cavity QED systems have long provided the ability to engineer strong non-local interactions between atoms, where virtual photons in the cavity mode mediate correlated spin-flips between positionally fixed atoms whose ground states serve as pseudo-spin-$\frac{1}{2}$ degrees of freedom [@Miller2005; @Walther2006; @Ritsch2013].  While single-mode cavities naturally generate uniform all-to-all couplings between spins [@Sorensen2002; @Strack2011; @Leroux2010; @Hosten2016], additional experimental tools may be used to engineer non-local interactions on more general graphs. 

Spin models on sparse graphs, for instance, can be obtained in a single-mode cavity by suppressing interactions with a magnetic field gradient and selectively enhancing interactions between pairs of atoms at particular distances by driving the ensemble with multi-frequency light [@Hung2016].  Similarly, by coupling an ensemble of atoms to a multimode cavity and controlling their positions with, e.g., optical tweezers [@Endres2016], one may obtain spin models whose interactions correspond to a complete graph, but whose strength and sign are disordered, similar to random couplings in the SYK model [@Gopalakrishnan2011; @Vaidya2018].  While these tools would enable the implementation of scrambling dynamics, several groups have also demonstrated techniques to measure OTOCs via interferometric protocols [@Swingle2016; @Garttner2017; @Meier2017; @Li2017], thereby allowing one to experimentally observe and characterize the growth of operators in these systems.  Owing to the presence of photon losses in cavity setups, one is always restricted by the effects of dissipation which limit the allowed time evolution of the system before decoherence sets in.  With state-of-the-art cavity systems deep in the strong coupling regime, however, we expect that it is possible to realistically observe scrambling in systems of as many as 100 spins in the near future. 

Acknowledgements {#acknowledgements .unnumbered}
================

We thank Geoffrey Bentsen, Patrick Hayden, Yuri Lensky, Xiao-Liang Qi, Steve Shenker and Brian Swingle for helpful discussions.  We also thank Tori Borish and Ana Maria Rey for providing feedback on the manuscript.  This work supported by the Gordon and Betty Moore Foundation’s EPiQS Initiative through Grants GBMF4306 (YG) and GBMF4302 (AL), and by the Office of High Energy Physics of the Department of Energy (GB). 

Bounding the Growth of Entanglement {#sec:subsystem}
===================================

The supporting information contains the technical details of the results announced in the first section.  Henceforth, we will assume that the Hamiltonian is $k$-local on a hypergraph with $L$ vertices, and that the local Hilbert space dimension is $2^M$. Here we are only interested in models where $M$ is strictly finite: e.g. $M=1$;[^1] thus $\log N \approx \log L$.

In this section, we sketch a proof that $\lambda_2$ is finite under rather general circumstances: (*i*) the $k$-local Hamiltonian (\[eq:introH\]) is extensive, and (*ii*) the spectrum of the many-body Hamiltonian approaches a smooth distribution in the thermodynamic limit.  The proof proceeds in two parts.  Let $\rho_v(t)$ be the reduced density matrix for vertex $v$ at time $t$, and define $$S^{(1)}_v(t) = -\mathrm{tr}\left[\rho_v(t)\log \rho_v(t)\right], \;\;\;\; S^{(2)}_v(t) = -\log \mathrm{tr}\left[\rho_v(t)^2\right].$$ First, we show that if $M\log 2 - S^{(2)}_v = \delta$, and $2^M\delta \ll 1$, $$\left| \left(M\log 2 - S^{(1)}_v\right) - \frac{\delta}{2}\right| {\leqslant}\frac{2^{(M-1)/2}}{3} \delta^{3/2} + \mathrm{O}\left(\delta^2\right).   \label{eq:S1S2}$$ Secondly, we show that $$\frac{ \left| \langle  \mathcal{O}_v(t) \mathcal{O}^\prime_v(0)\rangle \right| }{\lVert \mathcal{O}_v \rVert\lVert \mathcal{O}_v^\prime \rVert }{\geqslant}C \mathrm{e}^{-\lambda^*_v t},  \label{eq:lambdastar}$$ where $\mathcal{O}_v$ and $\mathcal{O}^\prime_v$ are local operators on vertex $v$, and $\lambda^*_v$ is a $v$-dependent decay rate which is independent of $N$ for almost all vertices. 

With these two results in hand, we now consider the von Neumann entanglement entropy of a region $A$ consisting of $L^\prime=pL$ of the vertices in the graph, with $0<p{\leqslant}\frac{1}{2}$. Since $$ML^\prime \log 2 - S_A {\geqslant}\sum_{v\in A} \left[ M\log 2 -  S_v^{(1)}\right] \approx \frac{1}{2} \sum_{v\in A} \left[ M\log 2 -  S_v^{(2)}\right],$$ and $$M\log 2 -  S_v^{(2)} \gtrsim \mathrm{e}^{-2\lambda_{2v}t},$$ if half of the vertices within $A$ have decay rates $\lambda_{2v} {\leqslant}\lambda_2$, then such vertices must (on average) have $\delta {\leqslant}\frac{4a}{L^\prime }$ in order for $ML^\prime \log 2 - \sum S_v  {\leqslant}a$. In the thermodynamic limit ($L\rightarrow \infty$ with $M$ finite), we will show that there exists a positive constant $C\propto N^0$ for which $\delta \gtrsim C \mathrm{e}^{-\lambda_2 t}$, which proves (\[eq:2lambda2\]). 

Relating Renyi to von Neumann Entropy
-------------------------------------

First we show (\[eq:S1S2\]).  If $0{\leqslant}p_i {\leqslant}1$ denote the eigenvalues of the reduced density matrix $\rho_v$, then $$S^{(1)}_v  = -\sum_i p_i \log p_i, \;\;\;\; S^{(2)}_v  = -\log \left(\sum_i p_i^2\right).$$ Also observe that $\sum p_i = 1$. Let $$p_i = \frac{1+\epsilon_i}{n},$$ and let $M\log 2 - S^{(2)}_v = \delta$. If $2\delta < 1$, then by definition $$\frac{1}{2^{M}} \sum_i \epsilon_i^2 = \mathrm{e}^\delta-1 < \delta+\delta^2 < 2\delta.   \label{eq:sumeps2}$$ Thus, for all $i$, $$|\epsilon_i| {\leqslant}\sqrt{2^{M+1}\delta}.  \label{eq:epsibound}$$ Now, when $2^{M+1}\delta \ll 1$, $$\begin{aligned}
S^{(1)}_v &= M\log 2 - \frac{1}{2^{M+1}}\sum_i \epsilon_i^2  + \frac{1}{6\cdot 2^M} \sum_i \epsilon_i^3 -  \frac{1}{12\cdot 2^M} \sum_i \epsilon_i^4 + \cdots,\end{aligned}$$ and so $$\left| \left(M\log 2 - S^{(1)}_v \right) - \frac{1}{2} \frac{1}{2^{M}} \sum_i \epsilon_i^2  \right| {\leqslant}\left(\frac{\max(|\epsilon_i|)}{6}  + \frac{\max(|\epsilon_i|)^2}{12}+ \cdots\right)\frac{1}{2^{M}} \sum_i \epsilon_i^2$$ Taylor expanding the equality in (\[eq:sumeps2\]), and using (\[eq:epsibound\]), we obtain (\[eq:S1S2\]). 

Renyi Entropy and the Memory Function
-------------------------------------

Next, we derive an exact, albeit non-local in time, equation of motion for $S^{(2)}_v(t)$. To do so, let us consider the vector space of all Hermitian operators which act on the Hilbert space of vertex $v$. The density matrix $|\rho_v)$ is in this vector space.  A basis set for this vector space is $\lbrace |\mathbb{I}), |T^a_v)\rbrace$, where $T^a_v$ denote the Hermitian generators of $\mathrm{SU}(2^M)$, acting on vertex $v$ alone.  An (orthogonal) inner product for this basis is $$(\mathcal{O}_1|\mathcal{O}_2) \equiv \mathrm{tr}\left(\mathcal{O}_1\mathcal{O}_2\right).$$

Our goal is now to find a linear equation for $|\rho_v(t))$. This can be done, as the Schrödinger equation is linear.  However, as we are integrating out degrees of freedom, the resulting equation will be nonlocal.  This is the essence of the memory function formalism [@mori; @zwanzig; @forster].  Unlike in the usual of the memory function equations (in real time), here we will also need to keep track of additional terms related to the initial conditions of the entire many-body state.  We will ultimately argue that these additional terms cannot generally modify (\[eq:lambdastar\]). 

The key point is as follows.  First, let us temporarily expand the vector space to include the set of all Hermitian matrices.  Our goal is to project the dynamics back on to the set of reduced density matrices on vertex $v$, which is an exponentially small subset of all possible basis vectors.  Let the initial many-body state be $|0)$. Since $\mathrm{tr} \rho = 1$ for all times, we know that $|0)$ and $|\mathbb{I})$ have some overlap, so we may write the initial state as $$|0) = |0_v) + |{\widetilde}0).$$ The first term is the initial condition for $\rho_v$, tensored with the identity on all other sites; the second term makes up the remainder of terms.  Let $\mathcal{L}$ be the Liouvillian: $$\mathcal{L}|\mathcal{O}) =  -|[H,\mathcal{O}]),$$ with $H$ the many-body Hamiltonian.  The time evolution of the density matrix is given by $\mathrm{e}^{-\mathrm{i}\mathcal{L}t}|0)$. Denoting $\mathfrak{p}$ as the projection operator onto $\lbrace |\mathbb{I}), |T^a_v)\rbrace$, and $\mathfrak{q}=1-\mathfrak{p}$, we thus wish compute $$|\rho_v(t)) = \mathfrak{p}\mathrm{e}^{-\mathrm{i}\mathcal{L}t} |0).$$ More specifically, since $$\mathrm{e}^{-S^{(2)}_v(t)} = ( \rho_v(t)|\rho_v(t)),$$ we wish to compute $$\frac{\mathrm{d}S^{(2)}_v}{\mathrm{d}t} = - \frac{2}{( \rho_v(t)|\rho_v(t))} ( \rho_v(t) |\frac{\mathrm{d}|\rho_v(t))}{\mathrm{d}t}.$$ Using the identity $$\mathfrak{p} \mathrm{e}^{-\mathrm{i}\mathcal{L}t}  =  \mathfrak{p} \mathrm{e}^{-\mathrm{i}\mathcal{L}\mathfrak{q}t} - \mathrm{i}\mathfrak{p} \int\limits_0^t \mathrm{d}s \mathrm{e}^{-\mathrm{i}\mathcal{L}\mathfrak{q}(t-s)}\mathcal{L} \mathfrak{p} \mathrm{e}^{-\mathrm{i}\mathcal{L}s},$$ we obtain $$\frac{\mathrm{d}|\rho_v(t))}{\mathrm{d}t} =   - \mathrm{i}\mathfrak{p}\mathcal{L}\mathfrak{q}\mathrm{e}^{-\mathrm{i}\mathfrak{q}\mathcal{L}\mathfrak{q} t} |{\widetilde}0) - \mathrm{i} \mathfrak{p} \mathcal{L}\mathfrak{p} |\rho_v(t))  - \int\limits_0^t \mathrm{d}s \; \mathfrak{p} \mathcal{L}\mathfrak{q} \mathrm{e}^{-\mathrm{i}\mathfrak{q}\mathcal{L}\mathfrak{q}s} \mathfrak{q}\mathcal{L} \mathfrak{p} |\rho_v(t-s)).  \label{eq:generaleqsec2}$$ This is a linear equation (albeit nonlocal in time) with a source, given by the first term on the right hand side.  We may write its solution as the sum of a particular solution and a homogeneous solution: $$|\rho_v(t)) = |{\overline}\rho^0_v(t)) + |{\widetilde}\rho_v(t))$$ The first term above is a homogeneous solution of (\[eq:generaleqsec2\]), obeying $\mathfrak{p}|{\overline}\rho^0_v)= |{\overline}\rho^0_v)$; the second term is a particular solution of (\[eq:generaleqsec2\]) obeying $|{\widetilde}\rho_v(0)) = 0$. For generic initial conditions, we expect the decay of $|\rho_v(t))$ to $\sim |\mathbb{I})$ can be no faster than the decay of the homogeneous term.  So we expect that it suffices to bound the decay of $|{\overline}\rho^0_v(t))$. This is the first point at which we will sacrifice some rigor and simply sketch out a proof of (\[eq:lambdastar\]).  The fact that a generic linear equation decays no faster than the decay of the homogeneous solution is the formal version of the inequality we provided in (\[eq:2lambda2\]). 

It is most natural to describe the dynamics of $|{\overline}\rho^0_v(t))$ via a Laplace transform.  It is convenient to subtract out the identity component: $$|{\overline}\rho_v(t)) \equiv |{\overline}\rho^0_v(t)) - \frac{1}{2^M} |\mathbb{I}).$$ Since $|\mathbb{I})$ is a null vector of $\mathcal{L}$, it decouples from the dynamics and we will omit this basis vector in what follows.  The assumption that entropy saturates to maximal then implies that ${\overline}\rho_v(t)) \rightarrow 0$ as $t\rightarrow \infty$. Letting $$|{\overline}\rho_v(z)) \equiv \int\limits_0^\infty \mathrm{d}t\; \mathrm{e}^{-zt} |\rho_v(t)),$$ we find that $$z|{\overline}\rho_v(z)) - |\rho_v(0)) = -\left(\mathrm{i}\mathfrak{p}\mathcal{L}\mathfrak{p} + \mathcal{K}(z) \right) |{\overline}\rho_v(z)),$$ where the memory function $$\mathcal{K}(z)  \equiv \mathfrak{p} \mathcal{L}\mathfrak{q} \left(z+\mathrm{i}\mathfrak{q}\mathcal{L}\mathfrak{q}\right)^{-1} \mathfrak{q}\mathcal{L} \mathfrak{p}.$$ The inverse Laplace transform then gives us the time-dependent reduced density matrix: $$|{\overline}\rho_v(t))  = \int\limits_{z_0+\mathrm{i}\mathbb{R}}  \frac{\mathrm{d}z}{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\mathrm{i}}  \;   \mathrm{e}^{zt}  \left(z + \mathcal{K}(z) + \mathrm{i}\mathfrak{p}\mathcal{L}\mathfrak{p}\right)^{-1} |\rho_v(0))$$ The real number $z_0$ is chosen somewhere where the contour can be closed: i.e. for $z_0{\leqslant}\mathrm{Re}(z_*)$, where $z_*$ is the complex number with smallest real part such that $z+\mathcal{K}(z)+\mathrm{i}\mathfrak{p}\mathcal{L}\mathfrak{p}$ is a singular matrix.  For the moment, let us assume that $z_*$ is finite.  We will justify this in the next subsection, together with $\mathrm{Re}(z_*)<0$. Standard theorems then give that (for generic initial conditions) as $t\rightarrow \infty$ $$\begin{aligned}
 |{\overline}\rho_v(t)) &\approx\mathfrak{p}_* \int\limits_{z_0+\mathrm{i}\mathbb{R}}  \frac{\mathrm{d}z}{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\mathrm{i}}  \;   \mathrm{e}^{zt}  \left(z + \mathcal{K}(z) + \mathrm{i}\mathfrak{p}\mathcal{L}\mathfrak{p}\right)^{-1}\mathfrak{p}_* |\rho_v(0)) = \int\limits_{z_0+\mathrm{i}\mathbb{R}}  \frac{\mathrm{d}z}{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\mathrm{i}}  \; \frac{\mathrm{e}^{zt}}{z-z_*} C_* |\rho_v(0)) = \mathrm{e}^{z_*t} C_* |\rho_v(0))
 \end{aligned}$$ where $\mathfrak{p}_*$ projects onto all vectors with this most singular ‘eigen’vector, and $C_* $ is a suitable matrix.  Subleading corrections to this equation decay with a faster exponential rate.  We then conclude that (\[eq:lambdastar\]) holds, with $\lambda^*_v = -2\mathrm{Re}(z_*)$ and $C = (\rho_v(0)| C_*^2 |\rho_v(0))$.

Finiteness of the Memory Function
---------------------------------

It remains to justify why $-\infty < \mathrm{Re}(z_*)<0$. First, let us suppose that $|E_\alpha)$ denote the eigenfunctions of $\mathfrak{q}\mathcal{L}\mathfrak{q}$ in the image of $\mathfrak{q}$. We write $$\mathfrak{p}\mathcal{L}\mathfrak{q} = \sum_\alpha   c_{\alpha E_\alpha}  |\psi_\alpha)  (E_\alpha|$$ Let us further assume that the density of states of the many-body Hamiltonian is continuous, and is given by $\rho(E)$. We then find that the density of states of both $\mathcal{L}$ and $\mathfrak{q}\mathcal{L}\mathfrak{q}$ is given by $$\rho_{\mathcal{L}}(E) = \int \mathrm{d}E^\prime \; \rho(E^\prime) \rho(E^\prime+E).$$ Because $\mathfrak{q}\mathcal{L}\mathfrak{q}$ and $\mathcal{L}$ are identical matrices up to an exponentially small number of entries, corrections to their relative continuous spectra are $\propto \mathrm{e}^{-L}$, and such effects are comparable, in the thermodynamic limit, to other “finite size" effects that we neglect.  We further assume that $$\rho_{\mathcal{L}}(E) \mathcal{F}(E)  \equiv \sum_{\alpha} {\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}(E_\alpha - E)  |c_{\alpha E_\alpha}|^2  |\psi_\alpha)  (\psi_\alpha|   \label{eq:FEeq}$$ is a smooth-valued positive definite matrix as a function of $E$, in the thermodynamic limit.[^2] If these two assumptions hold, we obtain $$\mathcal{K}(z) = \int \mathrm{d}E\;   \frac{\rho_{\mathcal{L}}(E) \mathcal{F}(E)}{z-\mathrm{i}E}$$

We first show that $-\infty < \mathrm{Re}(z_*)$. We do this by bounding $\mathcal{F}$, since $$(\alpha| \mathcal{K}(z) |\alpha) < {\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\max_E \left((\alpha|\rho_{\mathcal{L}}(E) \mathcal{F}(E)|\alpha).\right).$$ First, define $$\mathcal{L}_v = \mathcal{L}\mathfrak{p}.$$ For a $k$-local Hamiltonian, using the $\mathrm{L}_\infty$ operator norm (for simplicity): $$\frac{1}{L}\sum_{v\in V} \lVert \mathcal{L}_v  \rVert {\leqslant}\frac{2k\lVert H \rVert}{L} \propto L^0.    \label{eq:boundL}$$ The last step uses extensivity.  Using Markov’s inequality, this means that a finite fraction of vertices have $\lVert \mathcal{L}_v \rVert < M_0$, for an $L$-independent constant $M_0$ (as $M_0 \rightarrow \infty$, the fraction goes to 1).  From (\[eq:FEeq\]), we find $$(\alpha| \mathcal{F}(E)  |\alpha) <  (\alpha| \mathfrak{p} \mathcal{L}_v^2 \mathfrak{p} | \alpha) <  M_0^2  (\alpha|\alpha).$$ Hence$$\mathrm{Re}(z_*) > -M_0^2 \max_E \rho_{\mathcal{L}}(E) > -\infty.$$

To argue that $\mathrm{Re}(z_*)<0$, observe that $$\mathrm{Re}(\mathcal{K}(0)) = {\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\mathcal{F}(0) \int \mathrm{d}E \; \rho(E)^2.$$ $\mathrm{Re}(z_*)<0$ whenever $\mathcal{F}(0)$ is positive definite.  It is positive semidefinite by construction.  We do expect that all eigenvalues are strictly positive in a typical chaotic model.  One subtlety arises from the possiblity of Hamiltonians which have continuous many-body spectrum but are integrable.  These include either the quantum version of “classical" Ising Hamiltonians with integer coefficients, or many-body localized systems [@rahul].  For these systems, we do expect that $\mathcal{F}$ has null vectors or vectors that are so close to null that the above proof fails.  For example, in the many-body localized theories, one can generalize the above procedure to finite subsets of vertices.  Localization of eigenstates will apply that some eigenvalues of $\mathcal{F}$ are extremely small, when the subset of vertices is finite in the thermodynamic limit, but large compared to the localization length. 

Subject to the caveats described above, this completes the argument for the finiteness of $\lambda_v^*$ for most vertices.  Our assumptions seem quite plausible in chaotic many-body systems.  As stated in the main text, it would be interesting if the “proof" sketched here can be made rigorous. 

Possible Generalizations
------------------------

We now briefly comment on a few possible generalizations of our results.  Firstly, we expect that the bounding of entanglement generation at late times by the decay of correlation functions also holds at finite temperature, though a generalization of our “proof" appears non-trivial.  Nevertheless, we present some circumstantial evidence that hints towards a finite temperature generalization of (\[eq:2lambda2\]).  Combining results from [@solodukhin] and [@shenker13], we find that scrambling at late times is set by the decay rate of two point functions at finite temperature in two dimensional holographic conformal field theories.  Another example where the (partial) saturation of entanglement is related to the decay of two-point functions is the Sachdev-Ye-Kitaev model [@gu1708]. 

The presence of conservation laws, including energy conservation, may further slow down the growth of entanglement by preventing the rapid decay of mutual information between vertices within the subregion $A$ [@tibor2; @khemani].  Such effects further increase $t_{\mathrm{s}}$ beyond $\log L / (2\lambda_2)$, consistent with the fast scrambling conjecture. 

The growth of entanglement is extremely rapid at early times, and slows down at late times.[^3] This is consistent with bounds on the rate of entanglement generation, which show that the generation of entanglement is proportional to the perimeter of $A$ [@childs; @bravyi; @marien].  So we can formally bound the growth of entanglement at early times.  For simplicity, we will focus on a 2-local Hamiltonian with operator norms (\[eq:maincdef\]), though we will now let $\alpha,\beta$ run from $1$ to $2^{2M}-1$: [@bravyi; @marien] $$\frac{\mathrm{d}S_A}{\mathrm{d}t} < 18M\log 2 \times \frac{c}{K} \times |E_{A,A^{\mathrm{c}}}| < 18 cN\log 2.$$ In the above equation, $|E_{A,A^{\mathrm{c}}}|$ counts the number of edges between $A$ and $A^{\mathrm{c}}$. On a generic, locally tree-like graph, nearly every vertex is at the edge of $A$ [@bollobas], so we do expect $\mathrm{d}S_A/\mathrm{d}t\sim N$ at early times.  The fact that $\mathrm{d}S_A/\mathrm{d}t \propto N$ at early times was also noted in [@magan]; however, in [@magan], $a\propto N$ in (\[eq:SAmax\]), so that $t_{\mathrm{s}} \propto N^0$ is possible.  We view this $N$-independent time scale as a time scale for local thermalization, but not for scrambling. 

Lieb-Robinson Bounds {#app:LR}
====================

Derivation
----------

In this section we derive Lieb-Robinson bounds for *arbitrary* Hamiltonians of the form $$H = \sum_{S\subseteq \lbrace 1,\ldots,L\rbrace}  J_S  {\widetilde}{H}_S,   \label{eq:arbitraryH}$$ where the sum over $S$ includes all possible non-empty subsets of vertices, $J_S{\geqslant}0$, and ${\widetilde}{H}_S$ is an arbitrary unit norm Hermitian operator which acts non-trivially only on the Hilbert space $\bigotimes_{i\in S} \mathcal{H}_i$. Observe that the requirement that the Hamiltonian is $k$-local amounts to a requirement that $J_S=0$ if $|S|>k$. To derive our bound, we review the useful notion of a factor graph [@montanari].  A factor graph ${\widetilde}{G} = (V,F,E)$ consists of a vertex set $V=\lbrace 1,\ldots, L\rbrace$, a “factor vertex set" $F\subset \mathbb{Z}_2^V$, which consists of the possible $S$ in (\[eq:arbitraryH\]), and an edge set $E \subset V\times F$ containing all edges obeying the rule $$e = (v,S) \left\lbrace \begin{array}{ll}  \in E &\ v\in S \\ \notin E  &\ v\notin S\end{array}\right..$$ See Figure \[fig:factorgraph\]. 

![ An example of a factor graph.  Red circles denote vertices $\in V$, blue squares denote factor vertices $\in F$, and black lines denote edges.  A Hamiltonian which might be associated with this factor graph is $H=\sigma_1^z + \sigma_3^x+\sigma_4^x+\sigma_5^y + \sigma_1^z \sigma^2_x  + \sigma_2^x\sigma_3^y + \sigma_2^x \sigma^y_4 + \sigma^z_3 \sigma^z_5 + \sigma_4^x \sigma_5^y + \sigma_1^z\sigma_2^z\sigma_3^z + \sigma_2^x\sigma_3^x\sigma_4^x\sigma_5^z$.[]{data-label="fig:factorgraph"}](factorgraph.pdf){width="5in"}

As we will see below, it is convenient to take the factor vertex set $F$ to only include the sets $S$ for which $J_S > 0$, together with all $S$ of size 1 (corresponding to single vertices), but this is not strictly necessary.  In this case, the factor vertices denote terms in $H$, and the edges from $F$ to $V$ simply denote which vertices each term “acts" upon.  We can define a natural adjacency matrix $$\mathcal{A}_{vS} = \left\lbrace \begin{array}{ll}  1 &\   (v,S) \in E \\ 0 &\ (v,S)\notin E \end{array}\right. ,$$ and we also define the matrix $$\mathcal{J}_{SS^\prime} = {\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{SS^\prime} J_S.$$

Let us now return to the question of operator growth.  Consider a pair of operators $A_S$ and $B_Q$ which act non-trivially only on the vertices in $S\in F$ and $Q\in F$, respectively.  Define $$C_{SQ}(t) \equiv  \sup_{A_S, B_Q} \frac{\lVert [A_S(t),B_Q]\rVert }{\lVert  A_S\rVert\lVert B_Q\rVert},$$ with the supremum running over all possible Hermitian operators $A_S$ and $B_Q$ acting on the suitable Hilbert spaces.  Without loss of generality below, we take $t>0$ and $\epsilon>0$. Our goal is to find an inequality governing the growth of $C_{SQ}(t)$. Defining $$H_S \equiv  \sum_{S^\prime : S\cap S^\prime \ne \emptyset} J_{S^\prime} {\widetilde}{H}_{S^\prime}.$$ we observe that $$\begin{aligned}
\lVert [A_S(t+\epsilon),B_Q]\rVert &- \lVert [A_S(t),B_Q]\rVert  = \lVert [A_S(t),B_Q] - \mathrm{i}\epsilon  [[H,A_S(t)],B_Q]\rVert - \lVert [A_S(t),B_Q]\rVert   + \mathrm{O}\left(\epsilon^2\right)  \notag  \\
& {\leqslant}\epsilon \lVert [[H,A_S],B_Q(-t)]\rVert + \mathrm{O}\left(\epsilon^2\right) = \epsilon \lVert [[H_S,A_S],B_Q(-t)]\rVert+ \mathrm{O}\left(\epsilon^2\right) \notag \\
&{\leqslant}2\epsilon  \lVert A_S\rVert \lVert [H_S(t),B_Q]\rVert+ \mathrm{O}\left(\epsilon^2\right)\end{aligned}$$ We thus obtain $$\begin{aligned}
C_{SQ}(t+\epsilon) &= \sup_{A_S,B_Q} \left[ \frac{\lVert [A_S(t),B_Q]\rVert }{\lVert  A_S\rVert\lVert B_Q\rVert} + 2\epsilon \frac{\lVert [H_S(t),B_Q]\rVert }{\lVert B_Q\rVert} \right]+ \mathrm{O}\left(\epsilon^2\right)  \notag  \\
&{\leqslant}C_{SQ}(t) + 2\epsilon \sum_{S^\prime : S\cap S^\prime \ne \emptyset} J_{S^\prime} C_{S^\prime Q}(t)+ \mathrm{O}\left(\epsilon^2\right) \notag  \\
&{\leqslant}C_{SQ}(t) + 2\epsilon \sum_{S^\prime} |S^\prime \cap S|  J_{S^\prime} C_{S^\prime Q}(t)+ \mathrm{O}\left(\epsilon^2\right),\end{aligned}$$ which, upon taking $\epsilon  \rightarrow 0$, leads to the differential inequality $$\frac{\mathrm{d}C_{SQ}(t)}{\mathrm{d}t} {\leqslant}2   \mathcal{A}_{Sv} \mathcal{A}_{vS^{\prime\prime}}   \mathcal{J}_{S^{\prime\prime}S^\prime} C_{S^\prime Q}.   \label{eq:dCSQdt}$$ where we have employed the Einstein summation convention above.  In order to obtain (\[eq:dCSQdt\]), we have used the fact that the matrix $(\mathcal{A}^{\mathsf{T}}\mathcal{A} )_{SS^\prime} = |S\cap S^\prime|$ counts the number of vertices which the two sets $S$ and $S^\prime$ share in common.  Observe that both $C_{SQ}(0)$ and $\mathcal{A}^{\mathsf{T}}\mathcal{AJ}$ are non-negative matrices – i.e., all components are non-negative.  We can integrate (\[eq:dCSQdt\]) to obtain $$C_{SQ}(t) {\leqslant}\exp\left[ 2t \mathcal{A}^{\mathsf{T}}\mathcal{AJ}\right]_{SS^\prime} C_{S^\prime Q}(0) = \sum_{n=0}^{\infty} \frac{(2t)^n}{n!}  \left( \mathcal{A}^{\mathsf{T}}\mathcal{AJ}\right)^n_{SS^\prime} C_{S^\prime Q}(0).  \label{eq:CSQt}$$ By definition, we know that $$C_{SQ}(0) {\leqslant}\left\lbrace \begin{array}{ll} 2 &\  S\cap Q \ne \emptyset \\  0 &\ S\cap Q = \emptyset \end{array}\right. { \leqslant}2  (\mathcal{A}^{\mathsf{T}}\mathcal{A} )_{SQ}.$$

Suppose for simplicity that we are interested in the evolution of operators that act only on a single vertex at time $t=0$. In that case, the initial sets $S=\lbrace u\rbrace$ and $Q = \lbrace v\rbrace$, and we may write $\mathcal{A}^{\mathsf{T}}_{Su^\prime} Z_{u^\prime v^\prime} \mathcal{A}_{v^\prime Q} = Z_{uv}$ for any matrix $Z_{uv}$. Using this identity, together with (\[eq:CSQt\]), and denoting $C_{\lbrace u \rbrace \lbrace v\rbrace}(t) = C_{uv}(t)$, we obtain our most general Lieb-Robinson bound:$$C_{uv}(t) {\leqslant}2 \sum_{n=0}^{\infty} \frac{(2|t|)^n}{n!}  \left( \mathcal{AJ} \mathcal{A}^{\mathsf{T}}\right)^n_{uv} =2 \exp\left[ 2|t|  \mathcal{AJ} \mathcal{A}^{\mathsf{T}} \right]_{uv} \label{eq:LRgen}$$ Let us now focus on 2-local Hamiltonians with only $h_{uv}^{\alpha\beta}\ne 0$, as in the introductory section.  Using (\[eq:maincdef\]) we obtain $$\begin{aligned}
\left(\mathcal{AJ} \mathcal{A}^{\mathsf{T}}\right)_{uv} &=  \frac{1}{2}\sum_{u^\prime, v^\prime}  J_{\lbrace u^\prime v^\prime\rbrace } \mathcal{A}_{\lbrace u\rbrace (u^\prime,v^\prime)}\mathcal{A}^{\mathsf{T}}_{ (u^\prime,v^\prime)\lbrace v\rbrace} \notag \\
&{\leqslant}\frac{c}{2K} \sum_{u^\prime v^\prime \in E}\left({\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{uu^\prime} + {\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{uv^\prime}\right)\left({\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{vu^\prime} +{\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{vv^\prime} \right) = \frac{c}{K}(D+A)_{uv} .   \label{eq:LRgen2}\end{aligned}$$ In the sum on the first line, we put a factor of $\frac{1}{2}$ to account for the double counting of $u^\prime v^\prime$ and $v^\prime u^\prime$. Combining (\[eq:LRgen\]) and (\[eq:LRgen2\]) we find (\[eq:LR1st\]). 

Fast Scrambling on “Regular Hypergraphs"
----------------------------------------

In the main text, we described the constraints on fast scrambling in 2-local models.  To generalize this to $k$-local models, we can now consider the following simple argument.  Let us suppose that the number of $J_S>0$ is given by $\mathcal{N}$. For simplicity, let us assume that $$J_S {\leqslant}\frac{LJ_*}{\mathcal{N}},$$ with $J_*>0$, though it is straightforward to generalize the argument.  Let ${\widetilde}{1}_v = (1,\ldots,1)$. Since $$\left(\mathcal{JA}^{\mathsf{T}}\right)_{Sv} \frac{{\widetilde}{1}_v}{L}  {\leqslant}\frac{J_*}{\mathcal{N}} |S| {\leqslant}\frac{kJ_*}{\mathcal{N}}, \label{eq:JAT}$$ and thus $$\left(\mathcal{AJA}^{\mathsf{T}}\right)_{uv} \frac{{\widetilde}{1}_v}{L} {\leqslant}\frac{kJ_*}{\mathcal{N}}  \sum_{S:  u\in S} 1 = \frac{kJ_*}{\mathcal{N}} \mathfrak{n}_u,$$ where $\mathfrak{n}_u$ denotes the number of terms in $H$ that act non-trivially on $u$. This generalizes the degree $k_u$ of the vertex in the 2-local case.  On a hypergraph, where $\mathfrak{n}_u {\leqslant}m \mathcal{N}/L$ for every vertex (note $m{\geqslant}1$), we thus arrive at $$\left(\mathcal{AJA}^{\mathsf{T}}\right)_{uv} \frac{{\widetilde}{1}_v}{L} {\leqslant}kJ_* m \frac{{\widetilde}{1}_u}{L},$$ and hence $$\sum_{u,v} \frac{C_{uv}(t)}{L^2} {\leqslant}\frac{\mathrm{e}^{2kJ_* mt}}{L}.$$ If both $k$ and $m$ are finite in the limit $L\rightarrow \infty$, the exponent is finite and thus $C_{uv}$ only becomes O(1) after a time $t\sim \log L$, consistent with the fast scrambling conjecture. 

Towards the Traditional Bound
-----------------------------

We now derive (\[eq:LRdistmain\]), along with a generalization for arbitrary Hamiltonians.  The observation is simply that we can tighten (\[eq:JAT\]) to $$\left(\mathcal{AJA}^{\mathsf{T}}\right)_{uv} \frac{1_v}{L} {\leqslant}\frac{1_u}{L} \times   \left\lbrace \begin{array}{ll}  (kJ_* m)^n &\   n{\geqslant}d_{uv} \\  0 &\ n<d_{uv}  \end{array}\right. ,$$ where $d_{uv}$, the distance between $u$ and $v$, is the straightforward generalization of distance to the factor graph: half of the number of edges in the factor graph which must be traversed.  From (\[eq:LRgen\]): $$C_{uv}(t) {\leqslant}2\sum_{n=d_{uv}} \frac{(2t)^n}{n!}  \left( \mathcal{AJ} \mathcal{A}^{\mathsf{T}}\right)^n_{uv} < 2\mathrm{e}^{-d_{uv}}\sum_{n=d_{uv}}\frac{(2\mathrm{e}t)^n}{n!} ( kJ_* m)^n < 2\exp\left[2\mathrm{e}kJ_*mt    - d_{uv} \right].$$ Using (\[eq:LRgen2\]), we obtain (\[eq:LRdistmain\]).  More generally, we see that if $m$ (as defined in the previous subsection) and $k$ are finite, fast scrambling is only possible on hypergraphs with diameter $\lesssim \log L$.

Sachdev-Ye-Kitaev Models on a Graph {#sec:SYK}
===================================

In this section, we describe the generalized Sachdev-Ye-Kitaev model [@gu1609; @gu1702] on a random graph.  Let $\chi^v_i$ denote Majorana fermion $i$ on vertex $v\in V$: $\lbrace \chi^v_i, \chi^u_j\rbrace = {\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{ij}{\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}^{uv}$, and $i\in \lbrace 1,2,\ldots,M\rbrace$, and consider the Hamiltonian $$H =\mathrm{i}^{\frac{q}{2}}  \sum_u  \sum_{i_1< \ldots < i_q} J_{i_1,\ldots,i_q}^u \chi^u_{i_1} \ldots \chi^u_{i_q} +\mathrm{i}^{\frac{q}{2}} \sum_{uv  \in E} \sum_{\substack{i_1<\ldots<i_{q/2}\\j_1<\ldots<j_{q/2}}} {J}^{uv}_{i_1,\ldots,i_{q/2} j_1,\ldots,j_{q/2}} \chi_{i_1}^u \ldots \chi_{i_{q/2}}^u \chi_{j_1}^v\ldots \chi_{j_{q/2}}^v .$$ where we assume Gaussian random variables $J^u_{i_1,\ldots,i_q}$ and ${J}^{uv}_{i_1,\ldots,i_{q/2} j_1,\ldots,j_{q/2}}$ are independent and mean zero with the following variants: $$\mathbb{E}\left[\left(J^u_{i_1,\ldots,i_q} \right)^2 \right] = \frac{(q-1)!}{M^{q-1}} \frac{2^{q-1}}{q}  \left(1 - \frac{b k_u}{2}\right) J^2, \;\;\;\; 
\mathbb{E}\left[\left({J}^{uv}_{i_1,\ldots,i_{q/2} j_1,\ldots,j_{q/2}} \right)^2 \right] = \frac{(q/2)!^2}{ q M^{q-1}}  \frac{2^{q-1}}{q} b J^2.$$ The parameter $0{\leqslant}b{\leqslant}\frac{2}{k_{\text{max}}}$ sets the relative coupling strength on the links; ${J}$ sets the effective coupling constant. 

Large $q$ Limit and Saddle Point
--------------------------------

At large $M$, one solves for the correlation functions of the fermions by computing an effective action for the two-point Euclidean time Green’s function $$G_u(\tau_1,\tau_2) = \frac{1}{M}\sum_{i=1}^M \langle \chi^u_i(\tau_1)\chi_i^u(\tau_2)\rangle,$$ and an analogous self-energy $\Sigma_u(\tau_1,\tau_2)$. The result is $$\begin{aligned}
S_{{\operatorname{eff.}}}=\sum_{u,v} \left[ -\log {\operatorname{Pf}}\left( \partial_\tau - \Sigma_{u}  \right)  \delta_{uv} + \frac{1}{2} \int \mathrm{d}\tau_1 \mathrm{d}\tau_2 \left(
\Sigma_u G_u \delta_{uv} - \frac{J^2}{q} ( G_u^q \delta_{uv} -\frac{1}{2} b \Lambda_{uv} G_u^{q/2} G_v^{q/2}  )
\right) \right],  \label{eq:SYKSeff}\end{aligned}$$ where the matrix $\Lambda=D-A$ is the graph Laplacian, which generalizes a discretized $-\nabla^2$ to a general graph.  $\Lambda$ encodes all spatial dynamics of the effective action.  Because $\Lambda$ always has a null vector: $$1_v = \frac{1}{\sqrt{L}}(1,1,\ldots, 1),$$ the above action admits a simple saddle point, where $G_v = G_*$ and $\Sigma_v=\Sigma_*$ do not depend on vertex $v$. In the large $q$ limit, we expand $G_*$ and $\Sigma_*$ to $\mathrm{O}(1/q)$ for this saddle point: [@maldacena2016remarks]

$$\begin{aligned}
G_*(\tau) = \frac{1}{2} {\operatorname{sgn}}(\tau) \left( 1+ \frac{g_*(\tau)}{q} \right), \quad
 \Sigma_*(\tau) = \frac{{\cal J}^2}{q} \mathrm{e}^{g_*(\tau)},\end{aligned}$$

where the saddle point solution $g_*(\tau)$ is given by $$\begin{aligned}
\mathrm{e}^{g_*(\tau)}= \left[ \frac{\cos \frac{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\eta}{2} }{ \cos {\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\eta \left( \frac{1}{2} - \frac{|t|}{\beta} \right) }  \right]^2, \quad \beta {J}= \frac{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\eta}{\cos \frac{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\eta}{2}}\end{aligned}$$ where $\eta$ is a parameter determined by the coupling constant: in the high temperature limit $\beta {J}\ll 1$, $\eta\approx \frac{\beta {J}}{\pi}$; in the low temperature limit $\beta{J}\gg 1$, $\eta\approx 1- \frac{2}{\beta {J}}$.

In the discussion that follows, we will assume that $q$ is large, but not the largest parameter in the problem.  Namely, we will take $M\gg L \gg q$. We assume that the latter inequality may be safely taken, but have not explicitly checked this assumption. 

Out-Of-Time-Ordered Correlators
-------------------------------

By expanding around this saddle point, we can compute the connected piece of the averaged, regularized OTOC $$\sum_{i,j=1}^M \frac{ \mathrm{tr} \left(y \chi_i^u (0) y \chi^v_j(t_1)  y \chi^u_i(0) y \chi_j^v(t_2) \right)_{\text{conn.}} }{ M^2} \equiv \mathcal{F}_{uv}(t_1,t_2),$$ where $y = \rho_\beta^{1/4}$, with $\rho_\beta$ the thermal density matrix at inverse temperature $\beta$. Because of the combined the large $M$ limit and large $q$ limit, we can describe the exponentially growing behavior of $\mathcal{F}_{uv}(t_1,t_2)$ at all temperatures in this model.  Indeed, following [@kitaevunpublished; @maldacena2016remarks], we find that $\mathcal{F}_{uv}$ obeys the following linear equation in the exponential growth regime: $$\begin{aligned}
\mathcal{F}_{uv}(t_1,t_2) = \int \mathrm{d}t_3 \mathrm{d}t_4 \; K^{\mathrm{R}}_{uw}(t_1,t_2;t_3,t_4) \mathcal{F}_{wv} (t_3,t_4), \label{eqn: kinetic equation}\end{aligned}$$ where $K_{uv}^R$ is a retarded kernel with ‘spatial’ dynamics: $$K^{\mathrm{R}}_{uv}(t_1,t_2;t_3,t_4) = \frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}^2 v^2 \mathrm{\Theta}(t_{13})\mathrm{\Theta}(t_{24})}{\beta^2 \cosh^2 \left(\frac{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\eta}{\beta}t_{34}\right)} S_{uv},$$ and[^4] $$S_{uv} = {\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{uv} - \frac{b}{q-1} \Lambda_{uv}.$$ The exponentially growing ansatz has the following form: $$\mathcal{F}_{uv}(t_1,t_2) = \mathrm{e}^{\lambda_{\rm L} \frac{t_1+t_2}{2}} f_{uv}(t_{12}).$$ Because the vertex dependence in $K^{\mathrm{R}}_{uv}$ comes entirely through $S_{uv}$, we can solve for the spatial dynamics by studying (\[eqn: kinetic equation\]) for each individual eigenvector $\phi$ of $S$ with eigenvalue $s$. Suitable combinations of the solution to (\[eqn: kinetic equation\]) for these eigenvectors can be used to construct $\mathcal{F}_{uv}$ for any initial conditions.  Using the explicit form of the retarded kernel, applying a pair of derivatives, $\partial_1 \partial_2 $, to (\[eqn: kinetic equation\]), and denoting $t_{12}= \frac{\beta u}{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}v}$, we find: $$\begin{aligned}
\frac{\lambda_{\rm L}^2 \beta^2 }{4 {\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}^2 \eta^2 } \cdot f_s(u) = \left(\partial_u^2+ \frac{2s}{\cosh^2u } \right) f_s(u)\end{aligned}$$ This equation is the Schrödinger equation in a cosh-potential in one spatial dimension.  The following exact bound state solution is known, together with a corresponding eigenvalue parametrized by $a$: $$\begin{aligned}
f_s(u) \propto \frac{1}{\cosh^{a} u } \quad \Rightarrow \quad \frac{\lambda_{\rm L}^2 \beta^2 }{4 \pi^2 \eta^2 } = a^2, \quad 2s=a (a+1) .\end{aligned}$$ Importantly, we have now fixed the Lyapunov exponent $\lambda_{\mathrm{L}}$. Focusing for simplicity on the OTOC $\mathcal{F}_{uv}(t) = \mathcal{F}_{uv}(t,t)$, we conclude that $$\mathcal{F}_{uv}(t)  = \frac{1}{M}\sum_\lambda  C_\lambda \mathrm{e}^{\lambda_{\mathrm{L}}(s)t} \phi^\lambda_u \phi^\lambda_v,  
\label{eq:FuvC}$$ where $C_\lambda$ are undetermined constants, and $\phi^\lambda_v$ denote eigenvectors of $\Lambda$, and therefore $S$: $$\sum_{v}
S_{uv}\phi^\lambda_v =  s\phi^\lambda_u = \left(1-\frac{b\lambda}{q-1}\right)\phi^\lambda_u.$$ In the large $q$ limit, the eigenvalues of $S_{uv}$, $s$, are close to $1$. Therefore $$\begin{aligned}
\lambda_{\rm L}= \frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} \eta a 
\approx \frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} \eta \left(1-\frac{2b\lambda}{3(q-1)}\right).   \label{eq:lambdaLhighT}\end{aligned}$$

It remains to fix $C_\lambda$ to solve the problem.  Unfortunately, this is highly non-trivial in general, and depends on complicated details of the early time physics in the SYK model, e.g. on the dynamics at early time $0<t<\beta$ where (\[eqn: kinetic equation\]) doesn’t apply.  We do understand, however, the late time dynamics in a number of important limits, which we elucidate below. 

High Temperatures
-----------------

First, we describe the high temperature limit $\beta J \ll 1$. In this regime, $\eta\approx \frac{\beta {J}}{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}$ and $ \lambda_{\rm L}=2J\left(1-\frac{2b\lambda}{3(q-1)}\right)$. It is reasonable to expect that $\mathcal{F}_{uv}(t)$ is approximately local at early times $0<t<\beta$. To zeroth order, we then anticipate that (\[eq:FuvC\]) can be applied at $t=0$, and that $$\sum_\lambda C_\lambda \phi^\lambda_u \phi^\lambda_v \propto {\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{uv}.$$ Since the set of $\phi_u^\lambda$ form an orthonormal basis, we conclude that $C_\lambda$ is a constant, independent of $\lambda$. This implies that $$\mathcal{F}_{uv}(t)\propto \frac{1}{M} \exp\left[2 Jt \left(1-\frac{2b}{3(q-1)} \Lambda\right)\right]_{uv}.   \label{eq:FuvhighT}$$ We can interpret this result in a simple way: $\mathcal{F}_{uv}(t)$ is proportional to a concentration of “infected random walking individuals" located on vertex $u$ at time $t$, given that the only infected individuals were located on vertex $v$ at time $t=0$. The assumptions we make are that infected individuals grow at a constant rate of $2J$, and that infected individuals perform a random walk on $G$: traversing any given edge at a constant rate of $\frac{4bJ}{3(q-1)}$. This connection between the growth of chaos and the spread of infections has been observed for some time, and is also visible in the RUC. 

Graphs with Finite Spectral Gap
-------------------------------

Next, we relax the high temperature assumption, but require that the graph Laplacian $\Lambda$ has a finite spectral gap $\gamma>0$. The gap of $\Lambda$ implies a gap of Lyapunov spectrum, i.e. $$\mathcal{F}_{uv}(t)  = \frac{1}{M} \left(  C_0 \mathrm{e}^{\frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} \eta t} 1_u 1_v
+
C_1 \mathrm{e}^{\frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} \eta (1-\frac{2b}{3(q-1)} \gamma) t } \phi^\gamma_u \phi^\gamma_v + \ldots
\right)$$ In the long time limit, e.g. $t\gg q\beta$, the subleading terms are negligible comparing to the first term, related to the null vector $1_u= \frac{1}{\sqrt{L}}(1,1,\ldots,1)$.[^5] Therefore we have: $$\mathcal{F}_{uv}(t) \approx C_0 1_u 1_v \frac{\mathrm{e}^{\frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} \eta  t}}{M} = C_0  \frac{\mathrm{e}^{\frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} \eta  t}}{N}$$ In the last step we used $N=ML$. The constant $C_0$ is set by the overlap of the initial condition with the spatially uniform component of ${\mathcal{F}}_{uv}$, before the exponential growing regime, and we therefore expect $C_0$ to be an O(1) number.  Thus, for the graphs with finite spectral gap, we obtain $$t_* = \frac{\beta }{2\pi \eta}  \log N 
\label{eqn:tsr}$$ at leading order in $N$, independently of any details of the graph structure. 

As we are more interested in the spatial dynamics on the graph than in the dynamics of a single-site SYK model, we wish to take the limits $M,L\rightarrow \infty$ in such a way that $\frac{\log M}{L}\rightarrow 0$. In this case, it is necessary for the spectral gap $\gamma$ to remain finite even in the thermodynamic limit $L\rightarrow \infty$, in order to obtain (\[eqn:tsr\]).  Rather remarkably, it is a famous result in graph theory that $\gamma$ is strictly finite on any graph where the ‘perimeter’ of any subset of vertices scales proportionally to the number of vertices: see Appendix \[app:graphs\].  Therefore, on a typical sparse and locally treelike graph, the SYK model described above is just as chaotic as it would be on a fully connected graph.  This serves as an explicit example of a chaotic quantum system where some amount of sparsity to the connectivity graph does not affect the time scales of quantum information loss and operator growth. 

Low Temperatures
----------------

Another commonly discussed solvable limit is the low temperature limit $M\gg \beta J \gg 1$. In this limit, we do not need to assume that $q$ is large.  In this limit, the saddle point equation can be approximately solved by a conformal ansatz, where $K^{\mathrm{R}}_{uv}$ takes the following form (for simplicity, we set $q=4$ for this subsection): $$K^{\mathrm{R}}_{uv}(t_1,t_2;t_3,t_4) = 
\frac{3{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}\mathrm{\Theta}(t_{13})\mathrm{\Theta}(t_{24})}
{\beta^2 \cosh \left(\frac{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta}t_{34}\right) \left( \sinh \left(\frac{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta}t_{13} \right) \sinh \left(\frac{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta}t_{24}\right) \right)^{\frac{1}{2}}} 
S_{uv},$$ Similarly to the large $q$ discussion, $K^{\mathrm{R}}_{uv}$ can be diagonalized by the eigenvectors of $S_{uv}$, and eigenfunctions of the temporal part [@maldacena2016remarks; @kitaevunpublished], which determines the Lyapunov exponent through the following equation $$\begin{aligned}
\frac{3}{1+\frac{\beta}{{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}} \lambda_{\mathrm{L}}} \left( 1 - \frac{b\lambda}{3} \right)=1.\end{aligned}$$ We hence obtain $$\begin{aligned}
 \lambda_{\mathrm{L}} = \left( 1 -\frac{b\lambda}{2} \right) \frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta}.  \label{eq:blambdaLlowT}\end{aligned}$$ Therefore the gap $\gamma$ in $\Lambda$ also indicates a gap in the Lyapunov spectrum for the $q=4$ model at low temperature (a similar result also applies to arbitrary $q{\geqslant}4$).  Following the same argument as before, we can ignore the subleading terms and only focus on the null vector $1_u$ and leading exponent $\lambda_{\rm L}= \frac{2\pi}{\beta}$, which saturates the chaos bound [@stanfordbound]: $${\mathcal F}_{uv} (t)\propto \frac{\mathrm{e}^{\frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{ \beta} t}}{ML}.$$ independently of the graph structure. 

We can gain further intuition about the spatial dynamics by using an alternative treatment of the leading exponent in the low temperature limit.  Following [@gu1609; @gu1702], we describe the contribution to the leading exponent by the dynamics of reparametrization modes on every vertex.  The generalization of these works is straightforward but as the computation is rather technical we do not write it explicitly.  We find that the OTOC becomes $$\begin{aligned}
{\mathcal F}_{uv} (t) \propto \frac{1}{M} \left( \frac{\alpha}{\beta {J}}  + \frac{b}{3} \Lambda  \right)_{uv}^{-1} \mathrm{e}^{\frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} t}\equiv \frac{\mathrm{e}^{\frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta} t}}{M} Z_{uv}^{-1}. 
 \label{eq:FuvlowT}\end{aligned}$$ where $\alpha$ is a numerical constant.  The identity matrix is implicitly assumed to be multiplying constants in the expression above.  Note that this calculation assumes that $\lambda_{\mathrm{L}} \approx \frac{2{\text{{{\mbox{        \sbox{\foobox}{$\pi$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}}{\beta}$ for all modes; we see from (\[eq:blambdaLlowT\]) that this is analogous to approximating $b\rightarrow 0$.

At leading order in the small parameter ${1}/{\beta J}$, we reproduce the previous result $Z_{uv}^{-1} \approx \frac{\beta J}{\alpha L}$ independent of indices $u$ and $v$. Regarding the higher order effects, we can interpret $Z_{uv}^{-1}$ in terms of a simple statistical problem: at vertex $u$, we release one random walker per unit time onto the graph.  At any given instant in time, a random walker may walk along an edge of the graph to a neighboring vertex – this occurs with rate $\frac{b}{3}$. Finally, the random walkers die with rate constant $\frac{\alpha}{\beta J}$. Let $n_v$ be the expected number of random walkers on vertex $v$. In steady state, the rate of incoming walkers equals the rate of outgoing/dying walkers: $${\text{{{\mbox{        \sbox{\foobox}{$\delta$}        \hskip\wd\foobox
        \pdfsave
        \pdfsetmatrix{1 0 -.18 1}        \llap{\usebox{\foobox}}        \pdfrestore
}}}}}_{uv} + \frac{b}{3} \sum_{wv\in E} n_w = \left(\frac{\alpha}{\beta J}+\frac{b}{3}k_v\right)n_v.  \label{eq:randomwalk}$$ We conclude that $n_v = Z^{-1}_{uv} {\geqslant}0$. Interestingly, this does not have the interpretation of a diffusing random walker.  The meaning of the different mechanisms for the spread of chaos at low vs. high temperature is not clear to us, although its origins can be straightforwardly understood: the mechanisms responsible for chaos at high [@stanford1802] and low [@maldacena2016remarks] temperatures are very different.  See also the recent discussion in [@swingle1805], which argues that the high temperature behavior is more generic at finite $M$.

Random Unitary Circuit Model {#app:RUC}
============================

Mapping OTOCs to a Classical Stochastic Process
-----------------------------------------------

In this subsection, we review the stochastic process which the computation of OTOCs in the $m$-local RUC maps on to.  These results were all found previously in [@nahum]. 

First, consider a growing operator $\mathcal{O}_v(\mathrm{\Delta} t) = U\mathcal{O}_v U^\dagger$, with $\mathcal{O}_v$ a local operator which acts non-trivially only on vertex $v$. Recall that $\mathrm{\Delta}t=1/L$ in order to recover extensive quantum dynamics.  For simplicity, let us assume that $U$ will act on a subset $A$ of $m$ vertices: namely, $U$ is a random $2^{Mm}\times 2^{Mm}$ unitary matrix.  $U\mathcal{O}_v U^\dagger = \mathcal{O}_v$ if $v\notin A$.  If $v\in A$, after averaging over all $U$ with uniform measure, $U\mathcal{O}_v U^\dagger$ will be an equal superposition of all non-trivial $2^{2Mm}-1$ Hermitian operators acting on the subset $A$.  Observe that the fraction of these operators which act non-trivially on $\ell  {\leqslant}m$ vertices is given by $$p_\ell  =  \frac{(2^{2M}-1)^\ell}{2^{2Mm}-1}  \left(\begin{array}{c} m \\ \ell \end{array}\right).$$ If we fix vertices $u$ and $v$, and average over all $U$s acting on the subset $A$, and also average over local operators $\mathcal{O}_u$ and $\mathcal{O}_v$, we find $$\mathbb{E}\left[ \left\lVert \mathcal{O}_u, U \mathcal{O}_vU^\dagger \right\rVert^2   \right] \propto \left\lbrace \begin{array}{ll} \displaystyle  \sum_{\ell=1}^m \frac{\ell}{m}p_\ell &\ \lbrace u, v\rbrace \subset A \\ 0 &\ \text{otherwise} \end{array}\right..$$ The proportionality constant is related to the normalization of the $\mathcal{O}_v$ and is not needed for our purposes.  The right hand side should be interpreted as the weight in the operator $U\mathcal{O}_v U^\dagger$ which act non-trivially on $u$.

At later time steps, the above procedure generalizes straightforwardly.  For an operator $\mathcal{A} =  \sum \bigotimes \mathcal{O}_i $ which consists of a large number of complicated terms, $U \mathcal{A} U^\dagger$ can be evaluated term-by-term.  In each term $\bigotimes \mathcal{O}_v$, if there exists a vertex $v\in A$ for which $\mathcal{O}_v$ is not the identity, then we replace $\bigotimes_{u\in A} \mathcal{O}_u$ with a random sum of all possible $2^{2Mm}-1$ Hermitian operators with coefficients whose squares sum to one. 

This leads to the following observation: after averaging over all possible circuits, we can compute the proportion of the growing operator $U(t) \mathcal{O}_v U(t)^\dagger$ which acts non-trivially on the subset $S\subset V$ by mapping onto the following stochastic process.  Let $n_u \in  \lbrace 0,1\rbrace$ denote whether a vertex is “infected" or not; at time $t=0$, $n_v=1$ and all other $n_u=0$ for $u\ne  v$. At each time step, pick a random allowed subset $A\subset V$ of $m$ vertices.  If $\sum_{u\in A} n_u = 0$, then do nothing; otherwise – regardless of the microscopic state – with probability $p_\ell$, set $\ell$ vertices in $A$, chosen uniformly at random, to be infected, and the remaining $m-\ell$ vertices to be uninfected. 

In the limit $M\rightarrow \infty$, $p_m=1$. The infection only grows.  This is why, when $m=2$, the RUC maps onto a discrete time analogue of the SI epidemic model.  This is the example which we focused on in the main text. 

2-Local Operator Dynamics: Mean Field Methods and their Breakdown {#sec:RUC2}
-----------------------------------------------------------------

As we saw at infinite $M$ in the main text, the RUC wih 2-local dynamics maps on to the SI epidemic model which has super-exponential infection growth on heterogeneous networks.  In this section, we will show that this effect is an artifact of the $M\rightarrow \infty$ limit for 2-local RUCs.  The following section will give an example of super-exponential operator growth with a 3-local RUC. 

For now, we follow the literature on epidemics on complex networks [@vespignaniRMP] and use a mean field approximation to solve for the dynamics of the stochastic process of the previous section with $m=2$. The key approximation is a closed set of equations for $P_v \equiv \mathbb{P}[n_v=1]$: $$\frac{\mathrm{d}P_v}{\mathrm{d}t} = \frac{2}{K} \sum_{u\sim v} \left[(1-p_1)(1-P_v)P_u - p_1 P_v \right].$$ The first term corresponds to the rate at which the vertex $v$ becomes infected: this occurs with probability $1-p_1$ (as we do not care whether its neighbor gets uninfected in the process!),  and requires the neighbor to be infected while $v$ should be uninfected.  The latter term corresponds to the rate at which a vertex uninfects itself, which occurs at rate $2p_1/K$ per edge, regardless of the state of the neighbor. 

At early times, we may write the growth equation as $$\frac{\mathrm{d}P_k}{\mathrm{d}t}\approx \frac{2}{K} \left[(1-p_1)A_{vu}  - p_1 D_{vu}\right] P_u .$$ The growth rate of the epidemic is thus related to the spectrum of a particular linear combination of adjacency and degree matrices.  In order to make further progress, we resort to a further “degree-based" mean-field description [@vespignaniRMP].  Let $P_k$ denote the probability that $n_v=1$ for a vertex $v$ with degree $k$, and let $\rho_k$ be the probability that a randomly chosen vertex in the graph has degree $k$. We define $$\theta = \sum_k \frac{k\rho_k}{K}  P_k,  \label{eq:thetaPk}$$ which is the probability that a randomly chosen edge points to a node which is infected.  The mean-field approximation is that each node effectively sees each neighbor infected with probability $\theta$. With this approximation, $P_k$ obeys the closed differential equation $$\frac{\mathrm{d}P_k}{\mathrm{d}t} = \frac{2k}{K} (1-p_1) (1-P_k)\theta - \frac{2k}{K} p_1P_k.   \label{eq:dPkdt}$$ The overall factor of $2k/K$ arises because this is the rate at which a random edge which connects to a vertex of degree $k$ is chosen.  The first term in (\[eq:dPkdt\]) counts the rate at which the vertex is uninfected, and an edge is chosen between the given vertex and one of its neighbors: the final factor of $1-p$ counts the probability that the infection spreads to the vertex.  The second term counts the probability that a random edge is chosen, and this causes the infection to decay from the central vertex. 

At early times, we expect that $P_k \sim \mathrm{e}^{\lambda_{\mathrm{L}} t}$: namely, the largest eigenvalue of the linearized equation of motion dominates growth and sets a scrambling time, as defined by OTOCs.  This eigenvalue equation gives that $$\left(K\lambda_{\mathrm{L}} +  2kp_1\right) P_k = 2k(1-p_1)\theta  \label{eq:thetaPk2}$$ when $P_k \rightarrow 0$. Combining (\[eq:thetaPk\]) and (\[eq:thetaPk2\]), we find that $$1 = \sum_k \frac{k\rho_k}{K} \frac{2k(1-p_1)}{K\lambda_{\mathrm{L}} + 2kp_1}.   \label{eq:thetafin}$$

Let us begin with some formal bounds on $\lambda_{\mathrm{L}} $.  A lower bound on $\lambda_{\mathrm{L}} $ can be found from (\[eq:thetaPk2\]) by applying Jensen’s inequality on the convex function $x^2/(x+1)$ (for $x>0$): $$1>   \frac{2K(1-p_1)}{K\lambda_{\mathrm{L}} + 2p_1K}.$$ Hence $$\lambda_{\mathrm{L}}  > 2(1-2p_1).$$ A lower bound can be found by applying Jensen’s inequality to the concave function $x/(x+1)$, assuming the probability distribution on degrees $k\rho_k/K$, instead of $\rho_k$: $$\lambda_{\mathrm{L}}  < 2(1-2p_1) \frac{1}{K^2}\sum_k \rho_k k^2.$$ This latter probability upper bound is proportional to the growth rate of an SI epidemic on a heterogeneous graph [@vespignani]: $\frac{1}{K^2}\sum \rho_k k^2$. Unsurprisingly, infections spread at a reduced rate at finite $M$. This accordingly increases the time for the operator to grow on a regular graph by the factor $(1-2p_1)^{-1}$. As noted in [@nahum], this is analogous to the emergence of a butterfly velocity [@localized], which plays the role of an effective Lieb-Robinson velocity, accessible in correlation functions [@lrbutterfly]. 

In general, we expect that the upper bound above is better for heterogeneous graphs with a finite variance of the degree distribution, as it is known that epidemics do spread faster on heterogeneous networks, and our model reduces to the SI model as $p_1\rightarrow 0$. However, if $\sum_k \rho_k  k^2$ diverges, as it can on scale free graphs with $\rho_k \propto K^{\nu-1}k^{-\nu}$, with $2<\nu {\leqslant}3$, then this upper bound is lousy.  We can estimate the value of $\lambda_{\mathrm{L}}$ on these scale free graphs by observing that $$\begin{aligned}
1 &\propto \int\limits_1^{K\lambda_{\mathrm{L}}/2p_1} \frac{\mathrm{d}k}{k^\nu} K^{\nu-1} \frac{2k^2(1-p_1)}{K^2\lambda_{\mathrm{L}}} + \int\limits^\infty_{K\lambda_{\mathrm{L}}/2p_1} \frac{\mathrm{d}k}{k^\nu}  K^{\nu-1}  \frac{k(1-p_1)}{Kp_1} \propto   \frac{1-p_1}{p_1} \left(\frac{\lambda_{\mathrm{L}}}{2p_1}\right)^{2-\nu},\end{aligned}$$ which gives us that as $p_1\rightarrow  0$ $$\lambda_{\mathrm{L}}  \propto p_1^{-\frac{3-\nu}{\nu-2}} .$$ This suggests that for any finite $p_1$, the “epidemic spreading" dynamics of the RUC is fundamentally different from the SI model: in particular, $\lambda_{\mathrm{L}}$ remains finite. 

Indeed, for 2-local dynamics, we are not able to construct any explicit examples where we can prove that $\lambda_{\mathrm{L}} = \infty$. One way to try is to consider the “star graph" (the interaction graph of (\[eq:sigmazstar\])).  Here (\[eq:thetafin\]) becomes $$1 \approx   \frac{1-p_1}{2\lambda_{\mathrm{L}} + 2p_1} + \frac{ L(1-p_1)}{ 2\lambda_{\mathrm{L}}+2Lp_1 }.$$ This equation is approximately solved by $$\lambda_{\mathrm{L}} \approx L(1-2p_1);$$ namely, the mean-field description implies that $\lambda_{\mathrm{L}}$ diverges in the thermodynamic limit $L\rightarrow \infty$. However, on the star graph, we may also give a more explicit construction of the dynamics.  For simplicity, let us consider dynamics where the infection starts on the central node, 1.  At time steps $\propto 1/L$, we act on the central node (and another) with a random unitary; if the central node is infected, at each time step there is a finite probability $p_1$ that this node is “uninfected".  With high probability, at $t\propto  1/(p_1L)$, the central vertex becomes uninfected.  If $Q\propto p_1^{-1}$ of the other vertices were infected in this time frame, the time it takes for the central vertex to be reinfected is $\approx Q^{-1} \propto p_1$. So we can estimate that the time it takes to infect a finite fraction of vertices as $$t\approx p_1 + \frac{p_1}{2}+ \frac{p_1}{3} +\cdots \sim p_1 \log L.$$ This gives us a finite Lyapunov rate: $$\lambda_{\mathrm{L}} \approx \frac{1}{p_1}.$$ We have not mathematically shown that any 2-local RUC on any graph has a larger Lyapunov exponent than this.  It is possible that the mean field treatment of the dynamics on scale free graphs is also inaccurate. 

A 3-Local Model with $\lambda_{\mathrm{L}}=\infty$ {#sec:RUC3}
--------------------------------------------------

One may ask – is the observation that $\lambda_{\mathrm{L}} < \infty$ when $M$ is finite a general feature of the RUC, or is it an artifact of the 2-local dynamics above.  We now show that it is the latter.  Consider a 3-local RUC on a generalized star hypergraph consisting of $R$ “inner" nodes and $N$ “outer" nodes: we take $R/N \rightarrow 0$ in the thermodynamic limit.  Allowed unitaries in the RUC act on a single outer vertex and two inner vertices: all such pairs are allowed.  See Figure \[fig:3star\]. 

![ A 3-local RUC on a generalized star graph consisting of $R$ inner nodes (solid circles) and $N$ outer nodes (hollow circles).  The black triangle denotes a sample triplet of vertices on which a unitary can act.[]{data-label="fig:3star"}](3localRUC.pdf){width="2in"}

As on the star graph under 2-local dynamics, we can understand 3-local dynamics on the generalized star by focusing on the fraction $s$ of inner vertices which are infected at any given time.  Denoting with $q$ the fraction of outer edges which are infected, the probability that $s$ transitions at each time step is

$$\begin{aligned}
\mathbb{P}\left(s\rightarrow s+\frac{2}{R}\right) &= (p_3+p_2) q(1-s)\left(1-s-\frac{1}{R}\right) , \\
\mathbb{P}\left(s\rightarrow s+\frac{1}{R}\right) &= 2(p_3 + p_2) s\left(1-s\right) + 2N(p_1+p_2) q(1-s)\left(1-s-\frac{1}{R}\right), \\
\mathbb{P}\left(s\rightarrow s-\frac{1}{R}\right)  &= 2p_1 s(1-s) + 2N(p_1+p_2)s\left(s-\frac{1}{R}\right), \\
\mathbb{P}\left(s\rightarrow s-\frac{2}{R}\right)  &= p_1 s\left(s-\frac{1}{R}\right).\end{aligned}$$

At early times $s\propto R^{-1}$ and $q\propto N^{-1}$. Note that the rates of transition are the above probabilities multiplied by $N$ (the total number of vertices, in the thermodynamic limit).  The dominant infection-spreading at these early times are associated with unitaries that act on one infected inner node, and two uninfected nodes: one inner and one outer.  This means that $Rs$ undergoes a biased random walk from 0 to $R$. Let us now ask – given that at some time the random walk is at $Rs=m$, what is the probability $\mathfrak{p}_m$ that the random walk reaches $m=0$ before $m\sim R$?  This probability obeys the recursive equations $\mathfrak{p}_0=1$ and [@redner] $$\mathfrak{p}_m = \frac{p_3+p_2}{p_3+p_2+p_1} \mathfrak{p}_{m+1} + \frac{p_1}{p_3+p_2+p_1} \mathfrak{p}_{m-1}.$$ In the $R\rightarrow \infty$ limit, this set of equations is solved by $$\mathfrak{p}_m \approx \left(\frac{p_1}{p_2+p_3}\right)^m.$$ What this means is that once a single vertex in the inner core is infected, there is a rather high probability that the inner core never becomes uninfected.  Generalizing this logic, we conclude that for any finite $m$, there is a finite probability that $Rs {\geqslant}m$ for all future times.  What this means is that we may approximate $$\frac{\mathrm{d}s}{\mathrm{d}t} \approx \frac{N}{R} (p_3+p_2-p_1)s,$$ and so after an initial time $$t_{\mathrm{init}} \approx \frac{R\log R}{N(p_3+p_2-p_1)},$$ a finite fraction of the nodes inner core are infected.  So long as $(R\log R) /N \rightarrow 0$ in the thermodynamic limit, $t_{\mathrm{init}} \rightarrow 0$.

Once a finite fraction of the inner core is infected, then we may approximate $$\frac{\mathrm{d}q}{\mathrm{d}t} \approx (p_3 + 2p_2 + p_1)s_*(2-s_*)(1-q) - (2p_1+p_2) q,$$ where $$s_* = \frac{p_3+p_2-p_1}{p_3+2p_2+p_1}$$ is the fraction of inner vertices that are infected once $t\gg t_{\mathrm{init}}$. The time it takes for $q$ to be O(1) is finite in thus given by a constant in the thermodynamic limit.  By construction, we have shown that there exist 3-local RUCs with small local Hilbert space dimensions, for which $\lambda_{\mathrm{L}} = \infty$.

2-Local Entanglement Dynamics {#app:dynamics}
-----------------------------

Here we provide a more serious bound on the growth of entanglement for the random unitary circuit than the one found in the main text.  For simplicity, we focus on the 2-local case.  The discussion below follows [@nahum16], which studied the growth of entanglement on regular lattices.  As in the main text, we assume that our initial state is a product state.  To bound the growth of entanglement we will need the following two inequalities for the von Neumann entropy [@chuang]: $$|S_A - S_B| {\leqslant}S_{AB} {\leqslant}S_A +S_B.$$ Here $A$ and $B$ are disjoint subsets of the vertex set $V$. This implies that $$S_A(t) {\leqslant}S_{A+v}(t) + M\log 2, \;\;\;\;\;  S_A(t) {\leqslant}S_{A-v}(t) + M\log 2.   \label{eq:2entropyineqs}$$ Here $A+v=A\cup \lbrace v\rbrace$ and $A-v=A\cap \lbrace v\rbrace^{\mathrm{c}}$. The second useful fact that we will need is that $S_A[|\Psi\rangle] = S_A[ U_{ij}|\Psi\rangle]$ if either $\lbrace i,j\rbrace \subset A$ or $\lbrace i,j\rbrace \cap A = \emptyset$: this is easy to show by writing out the explicit formula for $S_A$ in terms of a partial trace. 

(\[eq:2entropyineqs\]) may be used to get sharper bounds on the growth of entanglement.  In particular, the result can be phrased as optimizing the following “cutting" problem on a graph: choose a (possibly) time-dependent subset $B(t)\subseteq V$ of vertices for which at every discrete time $t$, the unitary $U_{i_t j_t}(t)$ acts either entirely within $B(t)$ or entirely outside $B(t)$. At time steps where the subset $B(t)$ changes, we either remove or add one vertex to $B(t)$. Thus we conclude that $$S_B(t) {\leqslant}\left\lbrace  \begin{array}{ll} S_B(t-\frac{1}{L}) &\   B(t) = B(t-\frac{1}{L}) \\  S_B(t-\frac{1}{L}) + M\log 2 &\ B(t) \ne  B(t-\frac{1}{L}) \end{array}\right..  \label{eq:cutting}$$ For simplicity, we simply write $S_B(t)$ instead of $S_{B(t)}(t)$. We now look for fluctuating subsets $B(s)$ with the boundary condition $B(t)=A$, and arbitrary initial subset $B(0)$. As the initial state is a tensor product, $S_B(0)=0$, and so we conclude from (\[eq:cutting\]) that $$S_A(t) = S_B(t) {\leqslant}M\log 2  \times n_{\mathrm{cut}}[B(t)]$$ where $n_{\mathrm{cut}}$ denotes the number of time steps at which the subset $B$ changes.  See Figure \[fig:cutting\] for an illustration of the algorithm described above. 

![ Bounding the entanglement generation in a 2-local RUC with $L=8$ vertices.  After 5 time steps, we see that the set $B(t)$ has changed during 2 time steps, and so we conclude that $S_A(t) {\leqslant}2M\log 2$.[]{data-label="fig:cutting"}](entanglementcut.pdf){width="3in"}

The entropy $S_A(t)$ is clearly bounded by the choice of fluctuating subset $B(t)$ with the minimal number of cuts.  In the limit of large $M$, $S_A(t)$ was observed to be given *exactly* by the minimal number of cuts (up to a factor of $M\log  2$) [@nahum16]. 

The Cheeger Inequality {#app:graphs}
======================

In this appendix we show that $\Lambda$ has a finite spectral gap on any graph with sufficient “nonlocality".  The result goes by the name of the Cheeger inequality: the derivation below follows [@chung] and we review it here for completeness. 

Let $\phi_v$ be any eigenvector of $\Lambda$, and let $\lambda>0$ be the associated non-zero eigenvalue: $$\lambda \phi_v  = k_v\phi_v - \sum_{u:uv\in E} \phi_u.$$ The spectral gap is given by the smallest possible choice of $\lambda$, and so we will simply look for bounds on $\lambda$. Let $V_+\subset V$ be the vertices where $\phi_v > 0$. Note that $V_+\ne V$ and $V_+\ne  \emptyset$ as $\Lambda$ is a symmetric matrix with orthogonal eigenvectors, and $\Lambda$ has a null vector $(1,\ldots,1)$. We may write $$\lambda = \dfrac{\displaystyle\sum_{v\in V_+}  \phi_v \left(k_v{\widehat}{\phi}_v - \sum_{u:uv\in E} \phi_u\right) }{\displaystyle\sum_{v\in V_+} \phi_v^2}.$$ Defining $${\widehat}{\phi}_v \equiv \left\lbrace \begin{array}{ll}  \phi_v &\  v\in V_+ \\ 0 &\ \text{otherwise} \end{array}\right.$$ we obtain $$\begin{aligned}
\lambda &= \dfrac{\displaystyle\sum_{v}  {\widehat}{\phi}_v \left(k_v{\widehat}{\phi}_v - \sum_{u:uv\in E} \phi_u\right) }{\displaystyle\sum_{v} {\widehat}{\phi}_v^2} > \dfrac{\displaystyle\sum_{v}  {\widehat}{\phi}_v \left(k_v{\widehat}{\phi}_v - \sum_{u:uv\in E} {\widehat}{\phi}_u\right) }{\displaystyle\sum_{v} {\widehat}{\phi}_v^2}  = \frac{\displaystyle\sum_{uv\in E}({\widehat}{\phi}_u-{\widehat}{\phi}_v)^2} {\displaystyle\sum_{v} {\widehat}{\phi}_v^2} \notag \\
&> \frac{\displaystyle\sum_{uv\in E}({\widehat}{\phi}_u-{\widehat}{\phi}_v)^2} {\displaystyle\sum_{v} {\widehat}{\phi}_v^2}\frac{\displaystyle\sum_{uv\in E}({\widehat}{\phi}_u+{\widehat}{\phi}_v)^2} {2\displaystyle\sum_{v}k_v {\widehat}{\phi}_v^2}\end{aligned}$$ Now, using the identity (for real numbers) $a^2+b^2 {\geqslant}2ab$, and thus $$(a-b)^2(c+d)^2 + (a+b)^2(c-d)^2 {\geqslant}2\left|\left(a^2-b^2\right)\left(c^2-d^2\right)\right|,$$ together with the Cauchy-Schwarz inequality on the vectors ${\widehat}{\phi}_v$ and $\sqrt{D}{\widehat}{\phi}_v$, we obtain $$\lambda > \frac{1}{2} \left(\dfrac{\displaystyle \sum_{uv\in E} \left|{\widehat}{\phi}_u^2-{\widehat}{\phi}_v^2\right|}{\sum_v \sqrt{k}_v{\widehat}{\phi}_v^2}\right)^2.$$

Without loss of generality, label the vertices such that $\phi_1{\geqslant}\phi_2{\geqslant}\cdots{\geqslant}\phi_L$. Then $$\sum_{uv\in E} \left|{\widehat}{\phi}_u^2-{\widehat}{\phi}_v^2\right| = \sum_{u=1}^L \sum_{v>u} A_{vu} \sum_{\ell=u}^{v-1} \left({\widehat}{\phi}^2_\ell - {\widehat}{\phi}^2_{\ell+1}\right)$$ Denote with $\mathcal{E}_i$ the number of edges between $\lbrace 1,\ldots,i\rbrace$ and $\lbrace i+1,\ldots,L\rbrace$. Then $$\sum_{uv\in E} \left|{\widehat}{\phi}_u^2-{\widehat}{\phi}_v^2\right|  = \sum_{\ell=1}^{L-1} \mathcal{E}_\ell \left({\widehat}{\phi}^2_\ell - {\widehat}{\phi}^2_{\ell+1}\right)$$ In the mathematics literature, it is more helpful to express $\mathcal{E}_i$ in terms of a constant which we define as $$h_i \equiv \dfrac{\mathcal{E}_i}{\displaystyle \min \left( \sum_{v=1}^i k_v, \sum_{v=i+1}^L k_v\right)}.$$ $h_i$ denotes the fraction of edges that go between $\lbrace 1,\ldots, i \rbrace$ and $\lbrace i+1, \ldots, L\rbrace $.  In fact, the Cheeger constant, defined as $$h_G = \min_{A\subset V} \frac{|E_{A,A^{\mathrm{c}}}|}{\min(|E_A|,|E_{A^{\mathrm{c}}}|)}$$ where $E_{A,A^{\mathrm{c}}}\subset E$ consists of the edges between $A$ and $A^{\mathrm{c}}$, is well known to be *finite* in the thermodynamic limit for many families of random graphs [@chung].  This is a quantitative measure of how “locally treelike" the graph $G$ is.  Clearly, $h_G {\leqslant}h_i$. So $$\sum_{uv\in E} \left|{\widehat}{\phi}_u^2-{\widehat}{\phi}_v^2\right|  {\geqslant}\sum_{\ell} h_G\left({\widehat}{\phi}^2_\ell - {\widehat}{\phi}^2_{\ell+1}\right) \sum_{j{\leqslant}\ell} k_j = h_G \sum_\ell k_\ell {\widehat}{\phi}_\ell^2$$and hence $$\lambda > \frac{h_G^2}{2}  \left(\frac{\displaystyle  \sum_\ell k_\ell {\widehat}{\phi}_\ell^2}{\displaystyle   \sum_v \sqrt{k_v}{\widehat}{\phi}_v^2}\right)^2 > \frac{h_G^2}{2} \frac{k_{\mathrm{min}}^2}{k_{\mathrm{max}}}.$$ On any graph where $h_G>0$ and $k_{\mathrm{min}}/k_{\mathrm{max}}$ is a finite positive constant in the thermodynamic limit, this proves that the spectral gap of $\Lambda$ is strictly positive. 

[^1]: In SYK-like models, the local Hilbert spaces below should be interpreted as few fermion clusters on a given graph vertex, in contrast to our interpretation in (\[eq:tstarSYK\]), and in Appendix \[sec:SYK\]. 

[^2]: It is, by definition, positive semidefinite; we further assume that there are no null vectors. 

[^3]: Although saturation of entanglement entropy after a quench appears to be abrupt in translation invariant theories with a holographic dual [@suh1; @suh2], such models do not saturate (\[eq:introbound\]): the spatial spread of information is “slow". 

[^4]: See [@gu1609; @gu1702] for similar discussions and diagrams for spatially dependent kernels in generalized SYK models on lattices. 

[^5]: Technically, we also require $ML = N \gg \mathrm{e}^q$ for this long time limit to be sensible.