Many-body quantum dynamics slows down at low density

We study quantum many-body systems with a global U(1) conservation law, focusing on a theory of $N$ interacting fermions with charge conservation, or $N$ interacting spins with one conserved component of total spin. We define an effective operator size at finite chemical potential through suitably regularized out-of-time-ordered correlation functions. The growth rate of this density-dependent operator size vanishes algebraically with charge density; hence we obtain new bounds on Lyapunov exponents and butterfly velocities in charged systems at a given density, which are parametrically stronger than any Lieb-Robinson bound. We argue that the density dependence of our bound on the Lyapunov exponent is saturated in the charged Sachdev-Ye-Kitaev model. We also study random automaton quantum circuits and Brownian Sachdev-Ye-Kitaev models, each of which exhibit a different density dependence for the Lyapunov exponent, and explain the discrepancy. We propose that our results are a cartoon for understanding Planckian-limited energy-conserving dynamics at finite temperature.


Introduction
There is a conjectured universal "bound on chaos" [1] in many-body quantum systems: loosely speaking, a suitably defined out-of-time-ordered correlator (OTOC) at finite temperature is constrained to obey (We work in units with ħ h = k B = 1.) Originally, this rather abstract inequality was motivated by observations about quantum gravity [2]; indeed, saturation of (1.2) is believed to be achieved only by gravitational theories (described by many-body systems, in accordance with the holographic principle). However, from at least a heuristic perspective, this inequality is also sensible physically: at low temperature, the dynamics is restricted to increasingly few thermally activated degrees of freedom. The dynamics must necessarily slow down accordingly, and (1.2) (ignoring the precise 2π prefactor) is simply fixed by the Heisenberg uncertainty principle: This is one manifestation of a conjectured "Planckian" bound on quantum dynamics and thermalization, whereby the fastest time scale (at least, of thermalization) in a low temperature quantum system is 1/T . Heuristic evidence for quantum dynamics being limited by the time scale ħ h/k B T has arisen in many fields ranging from holographic field theories [3,4] to quantum critical theories [5,6], strongly correlated electrons [7][8][9][10] and quark-gluon plasma [11].
Of course, the argument (1.3) is far from rigorous, and strictly speaking there are plenty of counter-examples to (1.1), e.g. in free fermion models [1,12]. Is it possible, at least under certain circumstances, to prove that quantum dynamics truly must slow down at low energy? More broadly, can we show unambiguously that quantum dynamics has to slow down in any kind of constrained subspace? While this might seem intuitive, and there is certainly evidence for this [13][14][15], proving such a statement has been notoriously challenging, and very few rigorous results are known. The standard approach for constraining quantum dynamics is based on the Lieb-Robinson theorem [16], which applies to operator norms and holds in every state. By construction, therefore, Lieb-Robinson bounds are not useful at finding temperaturedependent bounds on quantum dynamics [17]. While recently these techniques have been improved to obtain temperature-dependent bounds on the velocity of information scrambling in one dimensional models [15], the resulting bounds depend on multiple microscopic model details.
It is almost certain that a rigorous derivation of (1.2), even in restricted models, is quite challenging without physical assumptions about thermalization. Clearly, a qualitatively important feature of low temperature dynamics is that it is restricted to low energy states in the Hilbert space. In this paper, we elect to study a simpler way to restrict dynamics to exponentially small parts of Hilbert space. Rather than cooling a system down to low temperature, we elect to study a system with a conserved U(1) charge, and in a highly polarized state with very low charge density n 1; here n denotes the probability that any lattice site is occupied. Elementary combinatorics demonstrates that exponentially small fractions of quantum states have n 1, even at infinite temperature. We will show explicitly how these constraints on accessible states qualitatively modify bounds on quantum dynamics and OTOC growth. In Section 3, we will show that in models where a Lyapunov exponent is well-defined, λ L ≤ λ * n γ , (1.4) where the exponent γ > 0 depends on basic details about the model (number of terms in interactions). There is a universal bound γ ≥ 1 2 , (1.5) valid for every theory; however, in certain cases, we can do parametrically better. (1.4) implies that in every theory with a U(1) conservation law (and a discrete Hilbert space), at least some kinds of quantum scrambling always become parametrically slow. The way which we derive this result is inspired by [12], which conjectured a similar phenomenon for energy conserving dynamics at finite temperature. However, our work will be more precise.
To understand the origin of the generic bound (1.5), let us recall that the growth of outof-time-ordered correlators comes from the "size" of operators increasing (see Section 2 for details). Consider for simplicity a model of fermions, with creation and annihilation operators c † i and c i . The simplest possible growth mechanism for a time evolving annihilation operator is (schematically) c i (t) = c i + j,k,l J i jkl c † j c k c l t + · · · . (1.6) Now, at low density, the second term above can only survive an expectation value if c † j acts on a state where site j is occupied. The fraction of states in the thermal ensemble where this site is occupied is n. So one might naively expect λ L ∼ n is the fastest possible operator growth, since each power of t will come with a factor of n (we can only add c † j and c k in "pairs", by charge conservation). However, since OTOCs such as (1.1) contain two commutators, the leading order density-dependent contribution to the OTOC will be 1 tr( ρc † i ρc i ) tr ρ j,k,l At the same time, we will see that the density dependence of (1.4) is qualitatively different in two models of quantum dynamics with time-dependent randomness: the Brownian SYK model [20,21] (Section 4) and a quantum automaton circuit (Section 5). In these models, the value of γ effectively doubles: γ → 2γ. As we will carefully explain, the discrepancy between the Hamiltonian quantum dynamics and the random time-dependent quantum dynamics is that the former relies on many-body quantum coherence effects, while the latter does not. This slowdown in effectively classical operator growth processes relative to quantum-coherent operator growth processes is reminiscent of the quadratic speed up of quantum walks over classical random walks [22,23].

Hilbert space
In this paper, we study quantum many-body systems with Hilbert space We interpret each copy of 2 = span(|0〉, |1〉) as consisting of either an empty site |0〉 or an occupied site |1〉. We define the density operators which measure whether the site i = 1, . . . , N is occupied, together with the total conserved charge The Hilbert space H can be written as the direct sum of subspaces with a fixed number of up spins: 11) with H N ↑ = span |n 1 n 2 · · · n N 〉 : Let U(t) be a time-dependent unitary transformation on H. In this paper, we are interested in studying the growth of operators when namely charge is conserved. The fact that charge is conserved means that the dynamics will separate the Hilbert space into N + 1 sectors corresponding to the allowed values of Q = 0, 1, . . . , N . In this paper, we will be interested in quantum dynamics in subspaces when Q and N are taken to be very large, while the ratio is held fixed. We will refer to n as the charge density, and focus on the limit n 1.

Operator dynamics and operator size
This paper is about operator growth: intuitively, given an operator such as n i which initially acts on only a finite number of sites (in this case 1), how much does time evolution scramble the information in n i ? Put another way, how complicated is the operator n i (t)? To answer this question carefully, we introduce a new formalism, following [24]. Let B denote the space of operators acting on H. For any operator A ∈ B, time evolution is defined by It is useful to write elements A ∈ B as "kets" |A) to emphasize linearity of quantum mechanics on operators, which will play a critical role in this paper. When t is a continuous parameter (i.e. we are not studying quantum circuit dynamics) we can also define the Liouvillian It is obvious that charge conservation places some constraints on operator growth. An operator which takes states of charge Q to states of charge Q will do so at all times. However, at finite n, there are O(exp(N )) such operators, so this constraint is not immediately useful or physically illuminating.
The purpose of this paper is to present a better way of thinking about operator growth in such systems, at a given density n. To do so, it is helpful to switch to a (grand) canonical perspective, and think about fixed chemical potential µ rather than fixed charge Q. Let denote the (grand) canonical density matrix at chemical potential µ, normalized so that tr(ρ) = 1.
Here we introduce the dimensionless chemical potential µ, which is the physically meaningful quantity in the infinite temperature limit. Note that (2.18) Given ρ, we now define the following inner product on B: For any value of µ, the length of an operator, which we define as (A|A), does not grow: Equipped with this inner product, we are now ready to define a physically useful notion of operator size and operator growth at finite µ. There are two possible interpretations of H, either in the language of spin models with a conserved z-spin, or in the language of fermion models with conserved charge. A non-local Jordan-Wigner-type transformation can convert between the two, but operator dynamics is not invariant under this transformation. For almost every quantum system [25], dynamics will only appear local in one language. So while there are clear similarities between how we talk about operator size and operator growth for a bosonic system vs. fermionic system, we must discuss each separately.
Let us first describe the physics when we interpret the Hilbert space in terms of bosonic degrees of freedom. First consider a single copy of 2 -i.e. a single site. There are 4 linearly independent operators acting on this two level system, forming the span of the operator vector space B i . With respect to the inner product (2.19), an orthogonal set of them is The lengths of these operators are For reasons that will become clear as we go through this paper, we define the size "superoperator" as a linear transformation on B i : Note also the useful identity e −µ/2 1 + e −µ = n(1 − n) . The generalization of (2.22) holds. The definition of size is now slightly more intuitive, as it counts the number of creation and annihilation operators: now denoting |n) = c † c − n: Thus far we have defined the size superoperator acting on a single Hilbert space, but it is straightforward to extend it to the N -body Hilbert space. Letting |T a ) (a = 1, . . . , 4) denote the four orthogonal operators above with length L a given in (2.22) and size S a given in (2.23) or (2.26) on a single two-level system, we observe that the following is an orthogonal basis for B: The length of each operator is Size is then defined as Let s denote a projector onto many-body operators of size s. Due to (2.20), we may define the probability that the operator A has size s at time t to be See [26] for a different interpretation of size at finite density or temperature. The probability that an operator has size s is related to the more convential out-of-timeordered correlation functions (OTOCs) which have been used to probe many-body chaos. As a simple example, let us consider a fermionic system, and ask for the typical magnitude of the normalized OTOC we conclude that this commutator will be non-vanishing whenever an operator string has either a c i or c † i on site i. Therefore, If the operator c j (t) did not have any strings with c † i c i on any site, then (2.32) would be an equality. We conclude that just as in the uncharged models [24,27], a typical OTOC C i j between a randomly chosen fermion i and our initial fermion j is non-vanishing only when the average operator size of c j (t) is large. However, crucially, this is when the operator size is measured with respect to the non-trivial inner product (2.19) at finite µ.

Bounds on dynamics
In the limit n 1, which corresponds to µ 1, we can estimate the canonical operators of size s (in our basis) as having length ∼ n s/2 . Recall that the "length" here refers to the Frobenius-like norm of the operator in the finite µ ensemble (2.19), whereas size counts the number of non-identity operators (with the inner product described above). As we now show, the fact that the canonical operators of size s have an exponentially small length leads to a significant slowdown in the dynamics of our operator size.

Lyapunov exponent
For illustrative purposes, we focus on Hamiltonian quantum dynamics generated by the fermionic q-body (also called q-local) Hamiltonian (note q must be even) If the J are all random, and are appropriately normalized, then this model is the complex SYK model [18,19,28] described in the next section. But we can also consider a more general class of models. Now let us start with the operator c j , as in our previous discussion. Our goal is to bound P s (t). In general, this is a challenging task [24], and requires finding the maximal eigenvalue of s L s . For illustrative purposes, it suffices to focus on what happens when s = 1 and s = q − 1. Without loss of generality, 1 consider the size 1 operator Observe that Now, the summation above does not depend on µ. The maximal eigenvalue of q−1 L 1 corresponds to maximizing the sum, which can be done independently of µ. Moreover, this argument did not depend on the choice of sizes s and s . We conclude that the maximal eigenvalue of s L s , denoted as s L s , obeys (3.37) Note the square root above, which arises due to the fact that there were two Ls in (3.36). Therefore, the growth of larger operators from smaller operators is parametrically slowed down at large |µ|, or when the system becomes low or high density. Note that (3.37) holds whether or not s < s or s > s . After all, for a charge conserving system, H commutes with ρ, and so (A|L|B) = −(B|L|A). s L s and s L s are transposes, and have the same maximal singular value (i.e. operator norm).
A quick route to justifying (3.37) is to simply observe from (2.22) that operators of size r always have their length reduced by a factor of (sech µ 2 ) r compared to what we might have naively expected based on the conventional Frobenius (µ = 0) inner product. Still, the reason that (3.37) is not trivial is that as we change the value of µ, the definition of |n) also changes, and so which operators have a given size must also change! After all, we know that the probability distribution P s (t) -which does have the µ-rescaled lengths built into it -is a well-defined probability distribution at any µ, and this would simply be impossible if our procedure was nothing more than re-scaling the lengths of strings of r c and c † operators by an r-dependent factor. The remarkable feature of all charge-conserving dynamics is that the operator evolution proceeds in just the right way so as to ensure the cancellation of two µ-dependent changes to our prescription: the change in the size-2 operator |n), and the change in the inner product (2.19).
Having now understood the physics behind (3.37), we can now immediately apply it to problems of interest. We begin by discussing the Lyapunov exponent of infinite temperature fermionic theories of the form (3.33) with a U(1) conservation law, defined by the exponential growth of (2.31): (3.38) for times smaller than the scrambling time t ∼ λ −1 L log N . As we described in (2.32), the growth of C i j (t) for generic i and j arises from the growth in effective µ-dependent size of the operator. Since the growth rates in size have been rescaled by (3.37), and due to charge conservation operator size can only grow by an even number, we can immediately find the universal bound Here λ * is a constant which comes from the µ = 0 bounds, following [24]. The n scaling above comes from applying (2.24) to (3.37). The point of this paper is not the evaluation of λ * , which can be quite challenging, but rather in the universal n dependence of (3.39).
A key ingredient in (3.39) is that operator size can only grow by an even number. In the fermion language, for example, we understand this as follows: the size 1 operators are c i and c † i . Suppose we have an operator A which transitions between the Hilbert spaces H Q 1 and H Q 2 at fixed charges Q 1 and Q 2 , as defined in (2.11). Each time we multiply by a product of creation/annihilation operators of size s, |Q 1 − Q 2 | (mod 2) changes by an amount s (mod 2). Hence, A(t) will only involve operators whose size is either all even or all odd. This result also holds in the spin language.
In certain models with all-to-all interactions, including the SYK model, we can parametrically improve upon (3.39). In the SYK model, operator growth in the large N limit is dominated by processes that grow operators by q − 2 c and c † at a time [24,27]. In this case, we can strengthen (3.39) to It is interesting that our approach readily leads to density-dependent bounds on Lyapunov exponents, which appear challenging to derive by other means [29]. We also emphasize that (3.40) does not depend on the precise choice of operators used in the OTOC.

Butterfly velocity
A more conjectural application of our rigorous result (3.37) is to constrain a suitably defined butterfly velocity. Consider a d-dimensional fermionic many-body system on a lattice of the form where the sum X runs over sufficiently local sets (e.g. no two sites in X 1,2 are farther than m sites apart, where m is some O(1) number). Charge conservation means that |X 1 | = |X 2 | in the sum above. Let us define v * B as follows: in a chaotic system, an operator c j (t) grows such that a typical OTOC C i j (t) is order 1 inside a ball of radius v * B t around site j. We propose for a generic system that there exists a The exponent q−2 4 above should be understood in the same context as (3.40): in the worst case scenario, we should set this exponent to 1 2 , however in certain models it may be possible to improve the exponent to q−2 4 . To (heuristically) obtain (3.42), we use a technique from [30]. For simplicity, assume that we have an operator supported at the origin of a d-dimensional lattice, O 0 , and that all terms in the Hamiltonian are either single-site fields, or nearest-neighbor interactions. We are interested in the weight of this operator on a site j a distance d j from the origin. So, let us define the superoperator F = j e d j j , (3.43) where j denotes the size of an operator on lattice site j. We now bound where R denotes a projection onto operators which have support on site i if and only if i ∈ R. Then, where we have used the Cauchy-Schwarz inequality and used the µ-dependent inner product similarly to (3.37). However, this step is not rigorous because we will not prove that nearly all weight corresponds to sets R and R that differ in q − 2 sites (or even just q − 2 fermions in the operator). Nevertheless proceeding with our argument, the key observation is that spatial locality on the lattice demands there exists a finite positive constant K such that for any two sets R and R contained in (3.44), (3.46) Moreover, the sum over sets R and R which share a union R ∩ R = i is finite. Therefore, we may replace the sum over R and R by a sum over lattice sites i: where z is another O(1) constant related to the number of sets X 1,2 in H containing site i, and K is yet another O(1) constant. Hence we conclude that However, comparing with (3.43), and using Markov's inequality [30], we conclude that which implies (3.42).
It is straightforward to generalize these results to spin models, rather than fermionic models. In models that are not of the form (3.33) and involve couplings with multiple different "q" (i.e. numbers of fermions), then the q in the above bounds should be replaced with the smallest value of q > 2 that appears in the Hamiltonian (since q = 2 terms do not grow operators).

Charged SYK model and its Brownian version
In this section, we consider two versions of the SYK model: one with Brownian motion couplings [20,21], and one with time-independent couplings. We will see that the behavior of operator growth in these two models qualitatively differs when n 1.

General methodology
The Brownian SYK and the regular SYK have the same form of the Hamiltonian (3.33) with a different random ensemble for the coupling constants (we upgrade the coupling J → J(t) to be generally time dependent): Note that in both cases above, J has the units of energy.
We now wish to calculate the Lyapunov exponent at fixed µ = βµ, as we take the infinite temperature limit β → 0. Unfortunately, neither of the analytically controlled limits of the SYK model -the strong coupling limit T J, or the large q limit -can be directly applied for our problem. After all, we are interested in β = 0 in this paper, invalidating the former approach. Moreover, if q → ∞, we can expect from (3.40) that for n < 1 2 , the calculation becomes trivial: λ L will vanish at leading order since for any 0 < c < 1, c q → 0 as q → ∞. This means the latter approach also is not directly useful.
Ultimately, we rely on an approximate method, which we expect will miss O(1) factors for the time-independent SYK model, but will otherwise be accurate. For the Brownian SYK model, however, our results will be exact. We use the Keldysh formalism with the assumption of a quasi-particle-like Green's function where Γ is the quasi-particle decay rate that will be self-consistently estimated. With the above approximated form of retarded Green function, we will find the Lyapunov exponent via the following kinetic equation [19,31,32]: where F R stands for the vertex function that contains the information of OTOC as a function of relative time (not the center of mass time which has been characterized by the λ L here). R(ω) is the rung function R = δΣ K /δG K obtained in the Keldysh formalism via the input G R we have in (4.51). We adopt a commonly used approximation [31] Therefore F R (ω) ≈ δ(ω + µ), and we have obtained here R(0) is the zero frequency component of the rung function.
Before applying the above formulas to the regular and Brownian SYK, let us clarify the validity of the approach used here. As is commented in Ref. [33], the above procedure is an approximation method for regular SYK because (1) in general the determination of the self energy at IR requires full knowledge of the Green function, not just the IR and UV limits. Therefore, the prefactor of the quasi-particle decay rate Γ is not expected to be accurate, while the scaling is still expected to be valid. (2) the approximation (4.53) also introduces inaccuracy for the prefactor of Γ . However, the above two sources of error will not occur for the Brownian SYK, because (1) the interaction is localized in time, so the self energy can be determined by the UV of the Green's function completely; (2) the relation (4.54) can be justified when R(ω) = R(0) is a constant in frequency. One way to see this is to rewrite (4.52) as follows 1 Integrating over ω for both sides, and eliminating the integral, we obtain (4.54).

Time-independent (regular) SYK
We first study the time-independent (regular) SYK model. The first step is to obtain Γ selfconsistently. By definition, we have iΓ = −Σ R (ω → −µ), and the retarded self energy can be obtained via Schwinger-Dyson equations (see Appendix A) and the the assumed form of G R .
In the limit β → 0 at fixed βµ = µ and Γ , we have whose ω → −µ component is Now, equating the above expression with −iΓ , we obtain However, as we commented above, we should not trust the constant prefactor in Γ for general q. 2 What is important is the dependence on J and µ, therefore for the rest of this subsection we will drop the unimportant prefactors. Next, we compute the rung function (4.59) Thus, the J and µ dependence for its zero frequency component is given as follows (4.60) Recall that within our approximation method, the Lyapunov exponent (4.54) is a linear combination of R(0) and Γ , so we conclude that

Brownian SYK
Next, we will move to the Brownian SYK where we will see a different scaling w.r.t cosh µ 2 . The computational logic for Brownian SYK is the same as for the regular SYK model; the only difference is that the interaction is uncorrelated in time. As a consequence, the two approximations in the above section become exact. For example, the self energy only relies on the UV behavior of the Green's function. Its Fourier transform 3 is a constant: (4.64) Thus, (4.65) Comparing with (4.58), we notice that the power law exponent of cosh µ 2 is twice that of the regular SYK model.
Similarly, the rung function (4.67) The Lyapunov exponent λ L = R(0) − 2Γ is therefore obtained as follows As commented before, this formula for the Brownian SYK is exact 4 , and we also note that the power is twice the result in the regular SYK.

Physical comparison between regular SYK and Brownian SYK
Let us now give a few physical arguments for the discrepancy between the Brownian/regular SYK models, as we believe this physics is somewhat universal (especially in models related to holographic gravity).
In the regular SYK model, we can loosely think of the density-dependence of λ L as follows. Consider a Taylor expansion of a time evolving operator, which looks schematically like In the first term of the Taylor series above, the operator has increased in size by q − 2 c and c † . By the formalism we developed above in (2.22), we know that each additional c and c † leads to an effective change in length of order n 1/4 . Recognizing that each subsequent commutator with H adds q − 2 more fermions, we can immediately see that the coefficient of c 1 (t) at order t k has length n k(q−2)/4 , which immediately implies (3.40).
Alternatively, if we are at low density n, then we can ask how many states there are which have a fermion on sites j 2 , . . . , j q 2 -the second term in (4.69) will annihilate any state where even one of those sites is unoccupied. At low density, the fraction of such states is n per site. So we might estimate the disorder-averaged average size to be (4.70) Again, the series above will be a function of t n (q−2)/4 . In the Brownian SYK, due to the time-dependent disorder average in (4.50b), we would instead find (4.71) Here, the series is a function of t n (q−2)/2 , which heuristically explains the doubling of the Lyapunov exponent. Ultimately, therefore, the difference between the Lyapunov exponents of the regular SYK model and the Brownian SYK model is the role of quantum coherence effects. Randomness in time, and not among the different coupling constants J, was responsible for the decoherence in the Brownian operator growth. This is analogous to the quadratic speed-up of coherent quantum walks over incohererent quantum walks, the latter of which behave identically to classical random walks [22,23]. Our universal bound (3.40) will be saturated by models, like SYK, with highly quantum coherent dynamics. It cannot be parametrically improved.

Butterfly velocity
We can generalize the discussions above to a spatially local version of the SYK model as introduced in [34,35]. Consider the Hamiltonian where the hopping matrix S x y = 0 only if x and y are nearest neighbors, or x = y. For example, in one dimension, we could take The coefficients J in (4.72) are defined so that H is Hermitian. On this simple one dimensional lattice, the eigenvectors of S x y are plane waves e ipx , with eigenvalues The growth of OTOCs in space can be characterized by the hopping matrix S x y above, which enters the kinetic equation (4.52) in the following way Note that the spatial and temporal dependence are factorized. Therefore, we can directly diagonalize the hopping S matrix using plane waves on the lattice. Within the approximation scheme we used before, we have the following p-dependent Lyapunov exponent where λ L (0) := λ L (p → 0) denotes the Lyapunov exponent we obtained in the case without spatial structure, while we remind that R(0) := R(ω → 0) is the zero frequency component not the momentum.
In the weak coupling, the butterfly velocity is determined by the saddle point of the following Fourier transform 5 from which we read out v 2 B = 4bλ L (0)R(0).

(4.78)
Regarding the dependence on the chemical potential/charge filling, we note that R(0) and λ L (0) have the same dependence as we demonstrated in previous sections, therefore we conclude that v B scales in the same way as λ L , namely .

(4.79)
This relation applies both to the regular and Brownian SYK. In particular, for the regular SYK, the above formula saturates the bound (3.42). It is easy to show that the discussion above for the nearest neighbor one dimensional lattice -in particular, the conclusion (4.78), generalizes to any other lattice.

Random automaton circuit
In this section, we discuss a random quantum automaton (QA) circuit, composed of N number of qubits (spin-1 2 degrees of freedom) with a global U(1) symmetry. Under QA dynamics, states expressed in the number basis (e.g. eigenstates of all Pauli Z operators) are sent to other eigenstates, without generating quantum superposition. Due to this special property, QA circuits can be simulated using the classical Monte Carlo algorithm. They have been extensively used to study quantum dynamics in both integrable and chaotic systems with local interaction [37][38][39][40][41].

Lyapunov exponent
Here, we construct a QA model consisting of k-qubit gates which acts on k qubits randomly selected in the system. This model has all-to-all interactions, and at each time step, we apply roughly N /k gates, to ensure extensive scaling of the dynamics in the large N limit. We expect that under time evolution, this QA model exhibits similar operator growth to a large class of other random circuit models with U(1) symmetry, including Haar random circuits without locality [42,43] and the Brownian SYK model above.
In the QA circuit, the k-qubit gate is randomly chosen to be U k with probability f or the identity with probability 1 − f . The U k gate is defined in the following way: For the k number of qubits, if the middle one has |1〉 with the rest k − 1 qubits having total 〈Z〉 = 0, U k will flip all these k − 1 qubits. It will leave other configurations invariant. The simplest case is k = 3, where we have Clearly, this circuit conserves total z-spin, and is U(1)-symmetric. Similarly, if k = 5, our QA circuit swaps between |00111〉 and |11100〉, |10110〉 and |01101〉, |01110〉 and |10101〉 and leaves other states invariant.
To understand the operator dynamics, we define the following operator basis for a single site: The space of many-body operators is a tensor product of this local basis. Any operator can be written as a superposition of these basis operators (Pauli string operators). This choice follows [41] and differs from the choice made in (2.21); however, due to the non-Hamiltonian nature of the QA circuit, this choice will prove a little more convenient here. Under the U k gate, a Pauli string operator maps to another Pauli string operator.
Let P N ↑ denote a projector onto the Hilbert space H N ↑ defined in (2.11). Consider the operator dynamics for X + x (t)P N ↑ in the limit n = N ↑ /N 1. In the operator basis defined in ((5.81)), X + x (t = 0)P N ↑ can be written as the superposition of Pauli strings with N ↑ P ↑ s and N − N ↑ − 1 P ↓ s. Under time evolution, the sum of the number of X + and P ↑ remains invariant, as does the number of X − and P ↓ together, due to charge conservation. Furthermore, the number of X + is always larger than X − by one. Operator growth can be characterized by counting the number of X + in X + x (t), which is 1 at t = 0 and eventually saturates to a value of order N ↑ .
Let us first assume k = 3 and define the number of X + as s. Under random dynamics governed by our QA circuit, at early time, the most important update rule for the growth of X + is P ↑ X + P ↓ → X + X + X − . Notice that the probability for P ↑ , X + and P ↓ are proportional to n, s and 1 − n respectively. Therefore in the continuous limit, we expect that ds dt ∼ ns, (5.82) which implies s ∼ exp(nt). The Lyapunov exponent λ is proportional to the ratio n. We can quickly generalize the above argument to any k. Since the probability to find (k − 1)/2 P ↑ s (which, at low density, is the limiting constraint) is proportional to n (k−1)/2 , the Lyapunov exponent obeys In order to compare (5.83) to our general bound (3.40), we observe that the Hamiltonian which would generate the U k gate (by applying it for a finite time) is schematically) We confirm this result numerically by computing OTOCs in the random QA circuit. For numerical ease, we study the following OTOC: where |s * 〉 = X i |s〉 flips a single spin/bit. We numerically computed C XZ (t) by averaging over the index i and j of C i j XZ (t). As shown in Fig. 1(a) for the case k = 3, C XZ (t) increases exponentially at early times. The Lyapunov exponent λ is linearly proportional to n when n 1. In Fig. 1(b), we show that the Lyapunov exponents of the k = 3, 5, 7 QA circuits are consistent with (5.83) for n 1.
We further computed the two point auto correlation function In terms of operator dynamics, this can be understood as the probability for the overlap between P N ↑ O i (t) and O i under time evolution, which should decay exponentially under the operator growth. As in the SYK models, we expect this decay rate is proportional to λ L . As shown in Fig. 1(c), we numerically computed the averaged C Z , and observed that Note that (1 − 2n) 2 is the saturation value C Z (∞). As shown in Fig. 1(d), we find that

Butterfly velocity
We have also studied the butterfly velocity v B in QA circuits where the degrees of freedom are arranged in a one-dimensional line [41]. The circuit for k = 5 is shown in Fig. 2; observe that the U k gates can now only act on a set of k adjacent degrees of freedom on the line: |s i+1 s i+2 · · · s i+k 〉. The QA circuits with k = 3 and k = 7 are constructed in an analogous way. In this case, there is no Lyapunov exponent due to the spatial locality. Nevertheless, we expect that v B ∼ n (k−1)/2 , (n 1).
The time-dependent randomness ensures that the n exponent "derived" in (3.42) must be multiplied by a factor of 2. Numerically, we computed v B by performing data collapse of the front of C XZ (r, t) ( See the example in Fig. 3(a)). We confirmed this prediction, as shown in Fig. 3(b).

Conclusions
We derived a new bound (3.40) on the growth of operators (as measured by OTOCs in a suitable (grand) canonical ensemble) in arbitrary many-body quantum systems. We studied several large N models with U(1) symmetry and showed that in the highly polarized sector with charge density n 1, the charged SYK model saturates our bound while the random dynamics including Brownian SYK model and random quantum automaton circuit do not. Due to the randomness in the time direction, the latter class of models lose the quantum coherence which allows the SYK model to saturate our bound. The Lyapunov exponents in these two classes of models satisfy the scaling relation and therefore classical systems are much less chaotic than quantum systems. Remarkably, a similar phenomenon to (6.90) arises in the study of systems with long-range interactions, where operator growth is much slower in effectively classical models [44,45] than in quantum coherent models [46][47][48].
There are a number of interesting applications and extensions of our work, which we briefly mention. Firstly, it is certainly interesting to try and generalize our results to other kinds of symmetry groups. An obvious candidate is SU(2) symmetry, which is easily realized in models of interacting qubits of the kind discussed in this paper. Such systems can approximately be realized in cold atomic gases [49], and our bounds may be relevant for designing models where highly entangled and metrologically useful states [50] exhibit very long lifetimes.
Secondly, we proposed a heuristic "bound" (3.42) on the butterfly velocity v B , which characterizes the growth of operators in a many-body model on the lattice. It would be interesting to make that argument more rigorous, if possible. More interestingly, it is worth investigating whether or not the density dependence of the butterfly velocity is captured by (3.42), or by the random unitary circuit models, which predicts (for a fermionic model such as SYK) v B ∼ n (q−2)/2 . (6.91) We postulate that, as in [41], the scaling (6.91) is more robust, as it incorporates destructive interference effects that seem natural for a typical chaotic system. Thirdly, recent work has used similar random circuits to model aspects of quantum gravity. Our work suggests that such analogies could be misleading for understanding short-time dynamics [21,27,35,51], because the mechanism for the exponential OTOC growth (1.2) is subtly different in a random circuit versus a holographic model. It would be interesting to understand better the crossover between the quantum coherent operator growth in the SYK model, and quantum incoherent operator growth in a random circuit, in particular to better understand quantum dynamics in chaotic lattice models. We also comment that random circuits have more recently been used to model holographic questions on much longer time scales, including the dynamics of a large and evaporating black hole [52,53]. Our work has no obvious relationship to this interesting problem.
Lastly, we note that other authors [54,55] has recently obtained the following bound for charged systems at chemical potential µ and temperature T : where µ c is a constant beyond which the (grand) canonical ensemble does not exist. We believe that this result, while it could be tight, is special to rotating black holes and their holographic duals. For example, the rotating three-dimensional black hole is dual to a twodimensional conformal field theory with holomorphic factorization, in which case T represents the harmonic mean of the left/right-moving temperatures (each of which controls a separate Lyapunov bound). In contrast, our result shows that (at least at infinite temperature) dynamics slows down by going to a constrained part of the Hilbert space. We expect that our results are much more universal, especially in non-holographic models. It would be interesting to generalize our result to finite temperature T , in which case a more detailed comparison with (6.92) could be made, along with other holographic results [56].
Re backwards, and also introduce α = 1, 2 . . . N for the contour indices. The interaction vertex is diagonal in (u, d) basis, so it will be convenient to first express the self energy in the (u, d) basis, and later make the basis change to the conventional (K, R, A) as follows: For complex fermions, we need to be careful about the arrows when drawing diagrams. G(t 1 , t 2 ) is represented by an arrow from t 2 to t 1 . We then find that Σ a b αβ (t 1 , t 2 ) = i t 1 , α, a t 2 , β, b = ±iJ 2 (iG a b αβ (t 1 , t 2 ))(G a b (t 1 , t 2 ) αβ G ba (t 2 , t 1 ) βα ) q−2 2 , +(−) for α = (=)β .
. (A.94) Here superscripts a, b ∈ {u, d} label the rail, the + sign is for a = b, and the − sign for a = b. Subscripts α, β = 1 . . . N label the contour index. The sign structure is due to the rule that each vertex is associated with a coupling constant: −iJ for u vertex, +iJ for d vertex.
Now we are ready to compute the self-energy which is diagonal in contour index. To proceed, we use the quasi-particle form (4.53) to obtain the following Green's functions, in the limit β → 0 with µ = βµ fixed: