A Universal Operator Growth Hypothesis

We present a hypothesis for the universal properties of operators evolving under Hamiltonian dynamics in many-body systems. The hypothesis states that successive Lanczos coefficients in the continued fraction expansion of the Green's functions grow linearly with rate $\alpha$ in generic systems. The rate $\alpha$ --- an experimental observable --- governs the exponential growth of operator complexity in a sense we make precise. This exponential growth even prevails beyond semiclassical or large-$N$ limits. Moreover, $\alpha$ upper bounds a large class of operator complexity measures, including the out-of-time-order correlator. As a result, we conjecture a sharp bound on Lyapunov exponents $\lambda_L \leq 2 \alpha$, which generalizes the known universal low-temperature bound $\lambda_L \leq 2 \pi T$. We illustrate our results in paradigmatic examples such as non-integrable spin chains, the Sachdev-Ye-Kitaev model, and classical models. Finally we use the hypothesis in conjunction with the recursion method to develop a technique for computing diffusion constants.


I. INTRODUCTION
The emergence of ergodic behavior in quantum systems is an old puzzle.Quantum mechanical time-evolution is local and unitary, but many quantum systems are effectively described by irreversible hydrodynamics, involving familiar quantities such as electrical conductivity.Understanding this emergent thermal behavior at both a conceptual and computational level is a central goal of theoretical research on quantum dynamics, of which a cornerstone is the eigenstate thermalization hypothesis [1][2][3][4][5].
Recent work has shifted focus from states to operator growth in many-body systems [6][7][8][9][10][11].Under Heisenbergpicture evolution, simple operators generically decay into an infinite "bath" of increasingly non-local operators.The emergence of this dissipative behavior from unitary dynamics is believed to be at the origin of thermalization, the decay of dynamical correlation functions, and the accuracy of hydrodynamics at large scales.This picture was recently confirmed in random unitary models of quantum dynamics [6,7], and extended to increasingly realistic systems involving conservation laws [8,9], Floquet dynamics [11], and even interacting integrable models [10].
While random unitary models are valuable proxies for studying operator growth, one would like to confirm this picture in genuine Hamiltonian systems.In semiclassical systems, a quantitative measure is provided by the outof-time-order correlation function (OTOC).The classical butterfly effect gives rise to an exponential growth of the OTOC, characterized by the Lyapunov exponent, which may be computed in a variety of models.It is conjectured that the Lyapunov exponent is bounded [12] and this bound is achieved in certain large-N strongly interacting models with a classical gravity dual, such as the Sachdev-Ye-Kitaev (SYK) model [13][14][15].Unfortunately, the OTOC does not necessarily exhibit exponential growth outside of semiclassical or large-N limits, rendering the Lyapunov exponent ill-defined [9,[16][17][18].A general theory of operator growth under generic, nonintegrable Hamiltonian dynamics is, therefore, still lacking.
The amount of information required to describe a growing operator increases exponentially in time.Computationally, this bars the exact calculation of operators at long times.Yet, the exponential size of the problem has a positive aspect: it acts as a thermodynamic bath, so a statistical description should emerge and become nigh-exact.This idea indicates operator growth should be governed by some form of universality.In this work we present a hypothesis specifying universal properties of growing operators in non-integrable quantum systems in any dimension.

II. SYNOPSIS
Our hypothesis has a simple formulation in the framework of the continued fraction expansion or recursion method, which we review in Section III.This is a wellunderstood technique, dating back to the 1980s [19], and has recently been used to compute conductivities in strongly-interacting systems [20][21][22].It is surveyed in great detail in Ref. [23].Essentially, it converts any linear-response calculation to a 1d quantum mechanics problem parameterized by the Lanczos coefficients.Section IV presents our hypothesis: operators in generic, non-integrable systems have Lanczos coefficients with asymptotically linear growth.The linear growth rate, denoted α, is the central quantity of this work.Although we are unable to prove the hypothesis rigorously, we shall support it with extensive numerical evidence, calculations in SYK models, and general physical arguments in Section IV.In particular, the hypothesis is equivalent to the exponential decay of the spectral function at high frequency, which can be (and has been) observed exper-imentally [24][25][26].
We explore several consequences of the hypothesis.In Section V, we develop a precise picture of the universal growth of operators.We show that under the hypothesis, the 1d quantum mechanics, governed by the Lanczos coefficients b n ∼ αn, captures the irreversible process of simple operators into complex ones.Furthermore, the 1d wavefunction delocalizes exponentially fast on the n axis, at a rate exactly given by α.Asymptotically, the expected position of the 1d wavefunction satisfies (n) t ∼ e 2αt . ( The expectation value (n) t has a succinct interpretation as an upper bound for a large class of operator complexity measures called "q-complexities", which we define in section V B. Crucially, this class includes out-of-time-order correlators.This allows us to establish a quantitative connection between α and the Lyapunov exponent, which will be the subject of Section VI.We conjecture that the growth rate gives an upper bound for the Lyapunov exponent whenever the latter is well-defined: This bound is equivalent to the universal one λ L ≤ 2πk B T / [12] in low temperature quantum systems, and remains non-trivial beyond this regime.We prove (2) partially, and confirm it in two examples: the SYK model and a classical tops model.The bound turns out to be tight in both cases.
A further application of the hypothesis, discussed in Section VII, is a semi-analytical technique to compute diffusion coefficients of conserved quantities.We leverage the hypothesis to extend classical methods of the continued fraction expansion to directly compute the pole structure of the Green's function, thus revealing the dispersion relation of the dynamics.We conclude in Section VIII by discussing conceptual implications of our results and perspectives for future work.

III. PRELIMINARIES: THE RECURSION METHOD
We briefly review the recursion method in order to state the hypothesis.A comprehensive treatment may be found in [23].Consider a local Hamiltonian H and fix a Hermitian operator O.We regard the operator as a state |O) in the Hilbert space of operators, endowed with the infinite-temperature inner product (O 1 |O 2 ) := Tr[O † 1 O 2 ]/ Tr [1].We write ||O|| := (O|O) for the norm.We will focus on systems in the thermodynamic limit.
Just as states evolve under the Hamiltonian operator, operators evolve under the Liouvillian superoperator L := [H, •].Our central object is the autocorrelation function C(t) = Tr[O(0)O(t)]/ Tr [1] = (O| exp (iLt) |O) , (3) where the second equality follows from Baker-Campbell-Hausdorff.One is interested in the long-time behavior of C(t), and also its frequency dependence (encoded in the spectral function).
Computing C(t) is inherently difficult.Suppose O(t = 0) is a relatively simple operator that can be written as the sum of a few basis vectors in any local basis [27].As the spatial support of O(t) grows, the number of nonzero coefficients of O(t) in any local basis can blow up exponentially.To make progress, one must compress this information.Intuitively, there are so many basis vectors at a given spatial size or "complexity" that we can think of them as a thermodynamic bath; no single basis vector has much individual relevance, only their statistical properties are important.In this interpretation, the operator flows though a series of "operator baths" of increasing size.The dynamics of an operator is then reduced to how the baths are connected -a much simpler problem.In particular, the second law then dictates that an operator eventually flows to the largest possible baths, running irreversibly away from small operators.This is shown schematically in Fig. 1.
Artist's impression of the space of operators and its relation to the 1d chain defined by the Lanczos algorithm starting from a simple operator O.The region of complex operators corresponds to that of large n on the 1d chain.Under our hypothesis, the hopping amplitudes bn on the chain grow linearly asymptotically in generic thermalizing systems.This implies an exponential spreading (n) t ∼ e 2αt of the wavefunction ϕn on the 1d chain, which reflects the exponential growth of operator complexity under Heisenberg evolution, in a sense we make precise in Section V.The form of the wavefunction ϕn is only a sketch; see Fig. 4 for a realistic picture.
We now quantify this idea precisely.This is done by applying the Lanczos algorithm, which iteratively computes a tridiagonal representation of a matrix.The idea is to find the sequence {L n |O)}, and then apply Gram-Schmidt to orthogonalize.Explicitly, start with a normalized vector |O 0 ) := |O).As a base case, let The output of the algorithm is a sequence of positive numbers, {b n }, called the Lanczos coefficients, and an orthonormal sequence of operators, {|O n )}, called the Krylov basis.(This is a bit of a misnomer, as the Krylov basis spans an operator space containing O(t) for any t, but does not usually span the full space of operators).
The Liouvillian is tridiagonal in this basis: We make four remarks.First, if the operator Hilbert space is d-dimensional with d finite (or if the subspace spanned by |O 0 ) , |O 1 ) , |O 2 ) , . . . is so), the algorithm will halt at n = d + 1: in this work, we work always in the thermodynamic limit and discard this trivial situation.Second, the Lanczos algorithm presented here is adapted to operator dynamics.Generally, a tridiagonal matrix will have non-zero diagonal entries, but they vanish in (5).This is because one can inductively show that i n O n is Hermitian for all n, hence (O n |L|O n ) = 0. Third, the knowledge of the Lanczos coefficients b 1 , . . ., b n is equivalent to that of the moments µ 2 , µ 4 , . . ., µ 2n , defined as The non-trivial transformation between them is reviewed in Appendix A. Fourth, the Lanczos coefficients have units of energy.
In the Krylov basis, the correlation function C(t) is: Hence the autocorrelation depends only on the Lanczos coefficients, and not on the Krylov basis.One way to interpret the Lanczos coefficients, which we will employ extensively below, is as the hopping amplitudes of a semi-infinite tight-binding model -see Fig. 1.
The wavefunction on the semi-infinite chain is defined as ϕ n (t) := i −n (O n |O(t)).Heisenberg evolution of O(t) becomes a discrete Schrödinger equation: where b 0 = ϕ −1 = 0 by convention.The autocorrelation is simply C(t) = ϕ 0 (t).Intuitively, we can think of the Krylov basis as stratifying operators by their 'complexity' (with respect to the initial operator O).In terms of the bath picture above, the "occupation number" of the nth bath is |ϕ n | 2 , and the b n 's describe the connections between adjacent baths.The goal of this work is to study aspects of operator growth that can be reduced to the quantum mechanics on this semi-infinite chain.

IV. THE HYPOTHESIS
We now state the hypothesis.Suppose that H describes an infinite, non-integrable [28], many-body system and O is a local operator having zero overlap with any conserved quantity (in particular, (O|H) = 0).Then the Lanczos coefficients are asymptotically linear: for some real constants α > 0 and γ.This linear growth is an example of universality.We will refer to α as the growth rate, and it will play a multitude of roles.In fact, it quantitatively captures the growth of "operator complexity" in a precise sense (Section V B).On the other hand, it is observable by standard linear response measures (Section IV C).Before presenting the evidence for the hypothesis, we note that the idea of classifying operator dynamics by Lanczos coefficients asymptotics is as old as the recursion method itself.Many examples have been explored, resulting in a broad zoology, as surveyed in [23].In particular, it is known that non-interacting models (such as lattice free fermions) give rise to a bounded sequence b n ∼ O(1).If we start with a two-body operator O in such free models, all O n 's will remain two-body.In this sense, the operator dynamics is simple.Models with obstructions to thermalization lead to more involved behavior.For example, a square root behavior b n ∼ √ n has been found in certain integrable models [23,29].To our knowledge, the ubiquity of asymptotically linear growth and its consequences have not been systematically studied in quantum systems.
A. Numerical Evidence: Spin Models Fig. 2(a) shows the Lanczos coefficients for a variety of spin models in the thermodynamic limit.(Numerical details are given in Appendix C.) One can clearly see that the asymptotic behavior is linear in all non-integrable cases.There is often an onset period before the universal behavior sets in; the first few Lanczos coefficients are highly model-dependent.We have observed that the more strongly-interacting the system, the sooner universal behavior appears [30].Fig. 2(b) shows the robustness of this linear behavior.The pure transverse field Ising model may be mapped to free fermions so, as expected, the Lanczos coefficients are bounded.But as soon as a small integrability-breaking interaction is added, the coefficients become asymptotically linear, and the linearity Lanczos coefficients in a variety of strongly interacting spin-half chains: H1 = i XiXi+1 + 0.709Zi + 0.9045Xi, H2 = H1 + i 0.2Yi, H3 = H1 + i 0.2ZiZi+1, and the SYK model (10) with q = 4 and J = 1.The initial operator O = √ 2γ1 for SYK; otherwise it is the energy density wave with momentum q = 0.1.(b) Cross-over to linear growth as interactions are added to a free model.Here H = i XiXi+1 − 1.05Zi + hX Xi, and O ∝ i 1.05XiXi+1 + Zi. bn is bounded when hX = 0 but asymptotically linear for any hX = 0. Numerical details are given in Appendix C. sets in at smaller n as the strength of the interaction increases.This is reminiscent of the crossover from Poisson to GOE distributed level statistics as integrability is broken [31,32].We have also checked a variety of other models believed to have chaotic behavior, and a number of operators in each.The hypothesis appears to hold in all cases.We may therefore conclude that the hypothesis is at least plausible.

B. Analytical Evidence: The SYK Model
It is an ironic point that the hypothesis (9) fails in virtually all exactly solvable models, as those are often integrable, even non-interacting.This explains why, to the best of our knowledge, the linear growth was not recognized in any of the extensive literature on the recursion method as a universal behavior (except for certain classical systems [33]).However, there is one solvable model where we can compute the linear behavior analytically: the SYK model (see, e.g.[13][14][15]).Its Hamiltonian is (10) where the γ i 's, with 1 ≤ i ≤ N , are Majorana fermions with anti-commutators {γ i , γ j } = δ ij , and the J i1...iq 's are disordered couplings drawn from a Gaussian distribution with mean zero and variance (q −1)!J 2 /N q−1 .We study the dynamics of a single Majorana O = √ 2γ 1 [34].To leverage the SYK solvability, we shall compute the moments µ 2n = O|L 2n |O , averaged over disorder in the large-N limit.For any finite q, the moments can be computed efficiently, thanks to the well-known large-N Schwinger-Dyson type equations satisfied by the correlation functions.The self-averaging properties of the SYK model allows the typical Lanczos coefficients to be computed from the averaged moments via a general numerical procedure [23].This is described in detail in Appendix B.
We find that the Lanczos coefficients follow the universal form (9) quite closely, as shown in Fig. 2(a).In the large-q limit, there is a closed form expression for the coefficients, computed in Appendix B: where J = √ q 2 (1−q)/2 J. Therefore in the large-q limit, the SYK model follows the universal form ( 9) with α = −2γ = J .We may conclude that our hypothesis is obeyed in a canonical model of quantum chaos.

C. Spectral Functions
Now that we have observed that the hypothesis holds in various models, we present a general argument supporting it.We consider here the hypothesis in light of the relation between the Lanczos coefficients and the spectral function.
Recall that the spectral function is the Fourier transform of the autocorrelation: In interacting many-body systems, the spectral function has a tail extending to arbitrarily high frequencies.The asymptotic behavior of the tail is directly related to the Lanczos coefficients, with faster growth of Lanczos coefficients corresponding to slower decay of Φ(ω).The precise asymptotic behavior is [35,36] b for any δ > 0 and some constant ω 0 .In particular, δ = 1 corresponds to asymptotically linear Lanczos coefficients and an exponentially decaying spectral function.
The decay of the spectral function is constrained by a powerful bound.A rigorous and general result of Refs [37,38] (see also [39,40]) is that, given an r-local lattice Hamiltonian H = i h i in any dimension, for some C > 0 and a known O(1) geometrical factor G r .We may conclude δ ≤ 1 in (13), so the Lanczos coefficients grow at most linearly.We provide a direct proof of this fact in Appendix F. In this sense, the hypothesis says that Lanzcos coefficients grow as fast as possible in non-integrable systems.When the hypothesis is satisfied, the growth rate α is quantitatively related to the exponential decay rate in the spectral function.In fact, Appendix A shows the following asymptotics are equivalent (see Fig. 3): We stress that this is a purely mathematical equivalence, which holds independently of physical considerations such as the dimension, the temperature, or even if the system is quantum or classical.However, this equivalence has a key physical consequence: it implies that α is observable in linear response measurements.In fact, high-frequency power spectra for quantum spin systems can be measured with nuclear magnetic resonance, and exponential decays were reported for CaF 2 [24][25][26].This experimental technique therefore provides a practical way of measuring α.Additionally, comparing (14) and (15) shows that α ≤ π/2κ, so the growth rate is limited by the local bandwidth of the model and the geometry.This inequality is a consequence of the natural energy scale for the Lanczos coefficients being set by the local bandwidth.However, we shall see that α itself is not merely the bandwidth, but contains a great deal of physical information about the system.
We find it useful to dispel a possible misconception.The hypothesis governs the high-frequency behavior of the spectral function Φ(ω).On dimensional grounds it is tempting -though ultimately erroneous -to interpret this as a statement about the short-time behavior of C(t).To see why this is wrong, notice that the shorttime behavior is captured by the first moment alone, as ).The high-frequency information instead governs the asymptotics of moments µ 2n as n → ∞ (which involve increasingly large operators) and the analytical structure of C(t) on the imaginary-t axis, as shown in Fig. 3.In particular, the exponential decay rate sets the location of the closest pole to the origin on the imaginary axis.The high-frequency information also does not control the large time limit t → +∞; we will come back to this point in Section VII B below.In brief, the hypothesis governs large ω behavior of Φ(ω) and, correspondingly, the behavior of C(t) on the imaginary axis.
To close this section, let us reiterate that the evidence presented constitutes no definitive proof that the hypothesis holds in all non-integrable systems.It is only certain that no faster growth is physically possible.Appendix F discusses the difficulty of showing the hypothesis from a microscopic point of view, and gives some partial results.It may be possible to put the hypothesis on a more solid footing using a type of random matrix model, a possibility we discuss in Section VIII.Nevertheless, we find the empirical evidence for the hypothesis compelling and expect it to be true.

V. EXPONENTIAL GROWTH OF COMPLEXITIES
Now that we have presented evidence in favor of the hypothesis, we shall turn to the analysis of its consequences.In this section we study the universal behavior of operators obeying the hypothesis.This is done in two steps.First, by studying the quantum mechanics problem (8) on the semi-infinite chain, we show that α measures the rate of exponential growth in operator complexity, in a sense we shall make precise below.Second, we prove that α gives an upper bound on a large class of operator complexity measures.
We remark that our notion of complexity is a priori distinct from other notions bearing the same name, such as circuit complexity (see the reviews [41,42] and references therein).Indeed, a satisfactory definition of operator complexity of any sort is an unresolved problem, and may not have a unique answer.

A. Exponential Growth in the Semi-infinite Chain
Recall that the Lanczos algorithm reduces the operator dynamics to a discrete Schrödinger equation ( 8),  (17) in the semiinfinite chain at various times.The wavefunction is defined only at n = 0, 1, 2 . . ., but has been extrapolated to intermediate values for display.
We shall analyze this quantum mechanics problem when the hypothesis is satisfied.
As a first step, we take the continuum limit, by linearizing around momenta 0 and π.This yields a Dirac equation ∂ t ϕ = ±2αx∂ x ϕ, whose characteristic curves x ∝ e ±2αt show the wavefunction spreads exponentially fast to the right in the semi-infinite chain with rate 2α.We remark that among all power-law Lanzcos coefficient asymptotics b n ∼ n δ , the linear growth δ = 1 is the only one which results in exponential spreading.When δ > 1, the characteristic curves reach x = ∞ at finite time [43].When δ < 1, the spreading follows a power law x ∼ t 1/(1−δ) .
To undertake a more careful analysis of the wavefunction on the semi-infinite chain, we employ a family of exact solutions.Suppose where η = 2γ/α+1.For any system when the hypothesis is satisfied, the b n 's will approach the b n 's asymptotically, so the properties of the exact solution using the b n 's are universal properties at large n.It is shown in Appendix D that the full wavefunction for the operator evolving under the b n 's is where (η is the Pochhammer symbol and |O n ) is the nth Krylov basis vector.Note that this example is not artificial but arises naturally in the SYK model, studied in Section VI B below.The exact solution (17) benefits from a detailed analysis.Recall that the component of the wavefunction at some fixed site n is ϕ n (t) = (−i) n (O n |O 0 (t)).For each n, ϕ n (t) is a purely real function which starts at 0 (for n > 1), increases monotonically until reaching a maximum at t ∼ log n, then decreases as ∼ e −αηt .The fact that exponential decay, reminiscent of dissipative dynamics, emerges under unitary evolution is quite remarkable, and is only possible in an infinite chain [44].Physically, the wavefunction is decaying by "escaping" off to n → ∞, which serves as a bath.Note, however, that the hypothesis is not sufficient to show that ϕ n (t) decays exponentially with time for small n, a fact whose consequences are studied in VII B below.
We now come to a central consequence of the linear growth hypothesis: the exponential spreading of the wavefunction.At any fixed time and large n, the wavefunction (17) has the form |ϕ n (t)| 2 ∼ e −n/ξ(t) , where ξ(t) is a "delocalization length" that grows exponentially in time: ξ(t) ∼ e 2αt for αt 1.This exponential spreading is reflected in the the expected position of the operator wavefunction ( 17) on the semi-infinite chain More generally, n k t ∼ e 2kαt for k ≥ 1.This result agrees, of course, with the one obtained in the simple continuum-limit above.We believe that the asymptotic growth in (18) holds whenever the Lanczos coefficients grow linearly.Although we have not proven this assertion, we have checked that it holds for many cases, such as artificially generated sequences of Lanczos coefficients b n = αn + γ n with various kinds of bounded "impurity" terms γ n ∼ O(1).We will consider (18) as a fact that follows directly from the hypothesis: the position of an operator in the abstract Krylov space grows exponentially in time.
We may interpret this exponential growth as a quantitative measure of the irreversible tendency of quantum operators to run away towards higher "complexity" [45].Indeed, we identify the position on the semi-infinite chain (n) t as a notion of operator complexity.We refer to (n) t as the "Krylov-complexity" (or "K-complexity" for short) of an operator After all, as n increases, the operators O n becomes more "complex", in the following sense: in the Heisenberg-picture, the equations of motions for O n 's form a hierarchy: that is, the dynamics of O n (t) depends on O n+1 (t).This is analogous to the BBGKY hierarchy in statistical mechanics, in which the evolution of the n-particle distribution depends on the (n + 1)-particle one.Similarly, as n increases, the O n 's becomes less local in real space, involve more basis vectors in any local basis, and are more difficult to compute.We remark that K-complexity is a distinct notion from precise terms such as circuit complexity and no relation should be inferred between the two.Closer precedents are the ideas of f-complexity and s-complexity [46].We know from Section IV C that linearly growing Lanczos coefficients are the maximal rate so, in turn, the wavefunction may not spread faster than exponentially.Thus the hypothesis implies that non-integrable systems have maximal growth of K-complexity: exponential, with rate 2α.

B. A bound on complexity growth
The physical meaning of K-complexity is far from transparent.After all, it depends on the rather abstract Krylov basis, the initial operator, and the choice of dynamics.To help pin down the idea of K-complexity, we study its relation to more familiar quantities.We shall consider a class of observables, "q-complexities" (q stands for quelconque), that includes familiar notions like out-of-time-order correlators and operator size.We will show that the growth of any q-complexity is bounded above by K-complexity.
We now define the q-complexity.Suppose Q is a superoperator that satisfies two properties: 1. Q is positive semidefinite.We denote its eigenbasis as |q a ), indexed by a, so that 2. There is a constant d > 0 such that Then q-complexity is defined to be the expectation value where O(t) is evolved under the Liouvillian dynamics of L. A q-complexity is, in principle, an observable, and requires Hamiltonian (or Liouvillian) dynamics.The rationale for the conditions is as follows: (20a) ensures the q-complexity is always non-negative, (20b) guarantees it cannot change too much under one application of the Liouvillian, and (20c) assigns a low complexity to the initial operator.To illustrate this concept, we now consider three examples: K-complexity, operator size, and out-of-time-order correlators.Example 2: operator size.A second example of a q-complexity is provided by operator size [34].For concreteness, we work in the framework of a spin-1/2 model (though Majorana fermions or higher spins work equally well).Consider the basis of Pauli strings, e.g.strings IXY ZII • • • with finitely many non-identity operators.Define Q to be diagonal in this basis, where the action of Q on a Pauli string is the number of non-identity Pauli's.So, for instance, Any choice of dynamics with at most d-body interactions (even long-ranged ones) will satisfy (20b), while (20c) requires simple that O is d-local.So, under these conditions, the q-complexity (Q) t becomes the average size of Pauli strings contained in O(t): Example 3: OTOCs.Our third -and most interesting -example of q-complexity is out-of-time-order commutators (OTOCs).Given O(t), there is an OTOC for each choice of a local operator For simplicity, we work with a many-body lattice model, and consider an on-site operator V i .We then define the OTOC superoperator by where the sum runs over all lattice sites i. Provided the Hamiltonian and initial operator are r-local, and that the dimension D of the on-site Hilbert space is finite, ( 23) is a q-complexity.
To see this, let us work in the eigenbasis of Q.For each site i, there is a basis Q i |q i,a ) = q i,a |q i,a ) with 1 ≤ a ≤ D 2 .We take |q i,0 ) to be the identity operator with eigenvalue 0, and note that 0 ≤ q i,a ≤ Q for some finite Q.Since [Q i , Q j ] = δ ij , the eigenbasis for the full operator space is the tensor product of the on-site bases.So for any sequence a = {a i }, |q a ) := ⊗ i |q i,ai ) is an eigenvector satisfying For the eigenvalue to be finite, a i must be zero for all but a finite number of i's and all eigenvalues are non-negative, so Q is positive semidefinite.Since the Hamiltonian is rlocal, the matrix element (q a |L|q b ) = 0 only if a and b differ on at most r sites.So by (24), we may bound the difference |q a − q b | ≤ d = rQ.Similarly, any r-local operator satisfies (20c).Having verified all the properties (20), we may conclude that OTOCs of this form are a qcomplexity.
OTOCs are known to be closely related to the operator size [12,34].It is usually possible to bound either quantity from the other, and to choose V i such that the OTOC reduces to the operator size.
We have now seen three examples of q-complexities, two of which are quantities that have been studied in recent times to understand the complexity of operators.A few remarks are in order.The q-complexities (including K-complexity) are quadratic in O(t) and not linear response quantities, although the growth rate α is, via the spectral function.The definition ( 20) is rather rigid and may be only sensible at infinite temperature; at finite temperatures it may be advisable to relax (20b), for instance.We will see in Section VI C that q-complexities may also be applied to classical systems, though they work somewhat differently there.
A rigorous argument in Appendix E proves that, for any q-complexity, The following section will focus on the application of this general bound in the specific case of OTOCs.

VI. GROWTH RATE AS A BOUND ON CHAOS
We showed in the preceding section that K-complexity provides an upper bound for any q-complexity whatsoever, which includes certain types of OTOCs.Combining ( 25) and ( 18), we see that q-complexities grow at most exponentially in time.If that is the case, with (Q) t ∼ e λ Q t , then the exponent is bounded above by 2α: In the rest of this section we focus on the case where the q-complexity is an OTOC.When the OTOC grows exponentially at late times, its growth rate λ L is called the Lyapunov exponent, since in the classical limit it reduces to the Lyapunov exponent characterizing the butterfly effect in classical deterministic chaos [47].The bound (25) then suggests the following conjectural bound on the Lyapunov exponent: We insist on calling Eq. ( 28) a conjecture because we shall boldly extrapolate (28) beyond where (25) and the arguments of the previous section apply directly; we have essentially proven (28) in quantum systems at infinite temperature already.Remarkably, (28) appears to be widely valid -at any temperature and in either classical or quantum systems.

A. Equivalence with the universal bound
We first examine finite temperature and show that our conjectural bound (28) is equivalent to the universal bound on the Lyapunov exponent of Ref. [12].The universal bound was derived for quantum field theories at finite temperature T = β −1 , and reads as follows in the natural units = k B = 1.
To compare (29) with our bound, we must adapt our formalism to finite temperature, which we discuss briefly here.The key novelty at finite temperature is that one must make a physical choice of operator norm from the infinity of mathematically-valid options.Any positive, even function g(λ) on the thermal circle defines a valid norm [23]: where y := e −H and Z := Tr[y β ] [48].Various choices of g(λ) have been made in the literature.For instance, linear response calculations have often dealt with the Kubo inner product g(λ) = 1; a simpler choice is afforded by . At infinite temperature they all reduce to the one (A|B) = Tr[A † B]/ Tr [1] standard in this work.
Here, we shall choose g(λ) = βδ(λ − β/2), which gives: This choice corresponds to inserting the operators A and B in the thermal circle [0, β) with even spacing, which is the same regularization scheme used for four-point OTOCs in [12] to argue for the universal bound.The autocorrelation function is modified accordingly In field theories, C β (t) is generally analytic in the strip |Im(t)| < β/2, and has contact singularities at t = ±iβ/2.This analytic structure determines the exponential decay rate of the spectral function Φ β (ω) ∼ e −|ω|/ω0 to be ω 0 = 2T , which implies 2α = 2πT (33) by the relation (15).Equation (33) implies the equivalence of the two bounds ( 28) and ( 29) at low temperature.However, (28) has a wider realm of validity.For instance, at infinite temperature or in classical systems, the universal bound becomes irrelevant, yet our conjecture remains nontrivial.

B. SYK Model
The clearest evidence for the conjecture (28) comes from the SYK model (10), where it holds at both infinite and low temperatures.
At low temperatures T = 1/β J, it is well-known that λ L = 2πT [14] saturates the universal quantum bound (29).The finite-T autocorrelation function (32) with O = √ 2γ 1 may be computed exactly by conformal invariance [13]: This is the autocorrelation function of the exact solution (17), and corresponds to Lanczos coefficients b n = πT n(n − 1 + η) with η = 2/q and α = πT , in agreement with (33).Therefore the low-temperature SYK model saturates our bound (28).At finite temperatures, using analytic results in the large-q limit [13], it is not hard to check (see Appendix B) that our bound is also saturated, whereas the universal bound ( 29) is not.
Returning to infinite temperature, no analytic formula for the Lyapunov exponent is available, but it has been computed numerically in, e.g.[13,34].Table I shows that not only does (28) hold for the whole range of q-SYK models, but α is almost equal to λ L /2, with exact agreement in the limit q → ∞ [49].These results show that the bound λ L ≤ 2α is tight: the prefactor cannot be improved in general.Moreover, in the large q limit, the probability distribution |ϕ n (t)| 2 on the semi-infinite line is identical to the operator size distribution of γ 1 (t) [34].(See (B17) in Appendix B for the precise statement.)So the large-q SYK model is an instance where the quantum mechanics problem on the semi-infinite chain can be concretely interpreted.
We remark that in models with all-to-all interactions like SYK and its variants may be the only circumstances where the bound (28) can be nearly saturated.For spatially extended quantum systems with finitely many local degrees of freedom, Lieb-Robinson bounds [50] and its long-range generalizations [51] guarantee that the OTOC has slower-than-exponential growth in most physical systems at infinite temperature [52].
Such a difference can be understood as follows.Due to the lack of spatial structure in the SYK model, we expect operator complexity (by any reasonable definition) is almost completely captured by operator size which, in turn, is directly probed by OTOCs.In finite-dimensional systems, complexity should be a distinct concept from operator size.For instance, long Pauli strings generated in the non-interacting Ising models have nonetheless low complexity, since they can be transformed to simple fewbody operators under the Jordan-Wigner transform.In non-integrable systems, by contrast, operator size growth is limited by Lieb-Robinson, while complexity can grow exponentially in the bulk of an operator's support.

C. Classical Chaos
We now transition to the classical setting.After briefly explaining how the recursion method carries over almost q 2 3 4 7 10 ∞ α/J 0 0.461 0.623 0.800 0.863 1 λL/(2J ) 0 0.454 0.620 0.799 0.863 1 TABLE I.The growth rate α versus half the OTOC-Lyapunov exponent λL/2 in the q-SYK model (10) in units of J = √ q2 (1−q)/2 J.Here α is obtained by exact numerical methods discussed in Appendix B, while λL is taken from the Appendix of [34].The q-SYK is physical only for even integers q, large-N methods allow an extrapolation to any q ≥ 2.
verbatim to classical systems, we shall examine the classical form of the conjecture (28).However, the arguments of Section V B do not carry over in full, and we are only able to prove a weaker bound.We close with a numerical case-study that suggests the stronger conjectural bound may well be true (and tight).

A (Weaker) Bound on Classical Chaos
The recursion method applies to classical and quantum systems in exactly the same manner [23].Classically, operator space is the (Hilbert) space of functions on classical phase space and the Liouvillian L = i{H , •} is defined by the Poisson bracket against the classical Hamiltonian H (we take = 1).The appropriate classical inner product at infinite temperature is (f |g) = f * g dΩ, where dΩ is the symplectic volume form on the phase space [53].The Liouvillian L is a self-adjoint operator, and the entire framework of the Lanczos coefficients carries over wholesale.
Indeed, the Lanczos coefficients have been studied more in the classical context.It is known [23,33] that linear growth of the Lanczos coefficients appears in general finite-dimensional, non-linear systems, to which we restrict ourselves [54].The growth rate α is well-defined in such systems, as is the (classical) Lyapunov exponent λ L , and the conjecture (28) takes on the same form as before: λ L ≤ 2α.In short, the similarity of classical and quantum Liouvillian evolution means that the recursion method -and its consequences -carry over unchanged.
There is, however, one important caveat: a classical OTOC does not generally qualify as a q-complexity.We will demonstrate this through an explicit, and instructive, example.Let us consider a single classical SU (2) spin.Its classical phase space is the two-sphere, and classical operator space is spanned by the basis of spherical harmonics |Y m ), = 0, 1, 2 . . ., m = − , . . ., .
A typical Hamiltonian is a polynomial of the classical spin operators S x , S y , S z with Poisson bracket {S a , S b } = −ε abc S c .We consider the simple nonlinear example Using Clebsch-Gordon coefficients one can show that the classical Liouvillian is quite sparse, and only the following matrix elements are non-zero: whenever the states in question exist.We now examine the classical OTOC for the local operator S z , given by matrix elements of a super-operator Q z .This operator is diagonal in the basis of spherical harmonics and we may immediately read off the eigenvalues as m 2 .When m changes by 1 upon application of the Liouvillian, the eigenvalue m 2 changes by 1 ± 2m, which can be arbitrarily large.Hence the condition (20b) cannot be satisfied for any finite constant d.It is helpful to recall that Section V B showed the quantum OTOC is a q-complexity whenever the on-site Hilbert space is finitedimensional.This fails in the case of a spin s, whose on-site dimension 2s + 1, in the classical limit s → ∞.
We have therefore seen that classical OTOCs are not qcomplexities and, hence, the conjectural bound (28) does not follow from the reasoning of Section V B in the classical case.Nonetheless, for any Hamiltonian and initial operators that are polynomials of the spin variables S a , we can show the following general bound which is weaker than the conjecture (28), λ L ≤ 2α.
To show (38), observe that by (37), the superopera- z satisfies (20b), since its has eigenvalue m for Y m , which can change only by δ upon one Liouvillian application, where δ is the polynomial degree of the Hamiltonian.Other conditions in (20) are satisfied straightforwardly.We then have which implies (38).Here the first ∼ is by definition, the the inequality is a straightforward generalization of the bound on q-complexity, Eq. (E8) of Appendix E, and the last ∼ is a generalization of (18) (see below that equation).This argument carries over to the OTOC with spin variable of any direction by spherical symmetry, and applies almost verbatim to systems with a few spins, S x,y,z i , i = 1, . . ., N .A Lyapunov exponent associated with a finite sum such as satisfies the same bound since every term does so.In summary, (38) is established in general classical few-spin models.We expect it is possible to show (38) rigorously.
An interesting corollary of ( 38) is a relation between chaos and the decay rate of the spectral function.Recall from Section IV C that the hypothesis is equivalent to the exponential decay of the spectral function Φ(ω) ∼ exp(−|ω|/ω 0 ) at high frequency, where ω 0 = 2 π α.Then (38) is equivalent to (The conjectured bound would instead imply λ L ≤ πω 0 .)In numerous classical systems, the power spectrum decay of time series has been used as an empirical probe of deterministic chaos [55][56][57][58][59][60][61].To the best of our knowledge, the bound (41) provides the first quantitative justification for this usage.
We mention that the relation between chaos and longtime decay of correlation functions has also been studied: long-time relaxation to equilibrium was shown to be controlled by Ruelle resonances in specific chaotic models [62,63].However, the long-time and high-frequency behaviors are a priori unrelated, as we discuss further in Section VII.
We stress that the growth rate is an upper bound on chaos, but not a diagnostic of it.Indeed, our bound is correct but not tight for most classical integrable systems which, generically, have non-zero growth rate but no chaos [33].
Unfortunately, we are not able to improve the argument and prove the stronger conjectured bound.Instead, we resort to testing the validity of the conjectured bound (28) in a canonical example of classical chaos.

Numerical study
The Feingold-Peres model of coupled tops [64] is a wellstudied model of few-body chaos, both classically and at the quantum level [65,66].The quantum model is a system of two spin-s particles, 1 and 2, with Hamiltonian where c ∈ [−1, 1] is a parameter and S α i satisfy the SU (2) algebra [S α i , S β j ] = i δ ij ε αβγ S γ i and act on a spins Hilbert space.This is non-interacting when c = ±1 and chaotic in the intermediate region.Correspondingly, the Lanczos coefficients are asymptotic to a constant at c = ±1 and increase linearly in intermediate regions.However, since the operator space dimension is finite (equal to (2s + 1) 4 ), the sequence of Lanczos coefficients is finite; in fact, the Lanczos coefficients saturate.The classical limit is obtained by taking s to infinity.There the Hamiltonian becomes where S α i , i = 1, 2 are two sets of classical SU (2) spins.As an SU (2) representation, the classical operator space contains all integer spins, whereas the quantum operator space has only integer spins up to 2s.We compute the classical Lanczos coefficients for the operator O ∝ S z 1 S z 2 (S z 1 S z 2 in the classical case).As shown in Fig. 5(b), the quantum Lanczos coefficients converge to the classical ones as s → ∞, as expected, and they increase linearly near c = 0. We have checked that α does not depend on the choice of initial operator O, so long as O does not overlap with any conserved quantity.
To test the conjectured bound (28), we compare the growth rate α with the classical Lyapunov exponent (λ L /2 in our notation), which can be calculated by the standard variational equation method [67].Remarkably, the data shown in Fig. 5(a) corroborates the conjectured bound α ≥ λ L /2 in the parameter region explored, with equality up to numerical accuracy in the regime c ≈ 0, where the model is known to be maximally chaotic, with almost no regular orbits [64,65].Enlarging the parameter space, for instance by adding terms such as S z i to the Hamiltonian, give further results consistent with the bound.It is thus possible that the conjectured bound is valid in classical systems and becomes tight in highly chaotic ones.

VII. APPLICATION TO HYDRODYNAMICS
Structural information about quantum systems can enable numerical algorithms.As an example, the success of the density matrix renormalization group algorithm is a consequence of the area law of entanglement entropy [68,69].We now apply the hypothesis to develop a semi-analytical technique to calculate decay rates and autocorrelation functions of operators and, in particular, compute diffusion coefficients of conserved charges.The key idea is to use the hypothesis to make a meromorphic approximation to the Green's function.This section introduces the continued fraction expansion of the Green's function, describes the zoology of operator decay, and finally presents the semi-analytical method.

A. Continued Fraction Expansion: Brief review
We briefly review the continued fraction expansion of the Green's function [23].The Green's function is defined as the expectation value and is related to the autocorrelation C(t) by the following transform: where the integration contour is taken to be the shifted real axis shifted down by −i for some small > 0. Since C(t) is bounded on the real axis, G(z) is analytic in the lower half-plane, but may contain singularities on the upper half plane.We shall refer to (45) as the Laplace transform, despite the fact that it differs from the usual definition by a factor of i.
In the Krylov basis, 00 corresponds to all paths that start on the first site, propagate through the chain, and return.We can divide all paths into those that stay on the first site, and those that first hop to the second site, propagate on sites n ≥ 2, and then return.More formally, for each n ≥ 0, let L (n) := L p≥n,q≥n be the hopping matrix on the semi-infinite chain restricted to sites n and above, and let nn be the corresponding Green function.(Note that G (0) (z) = G(z).)We then have the following recursion relationhence the name "recursion method" - (For a quick derivation [22], consider the polynomial P n (z) := det(z − L (n) ).By Cramer's rule we have G (n) (z) = P n+1 (z)/P n (z); a cofactor expansion gives P n (z) = zP n+1 (z)−b 2 n+1 P n+2 (z).Then (46) follows from the two preceding equations.) Applying Eq. ( 46) recursively yields the continued fraction expansion: To save space, we denote the recursion 46 by , where M n is the Möbius transform w → 1/(z − b 2 n w) and "•" denotes function composition.It is crucial that the convergence of the continued fraction expansions is quite subtle and quite different from the convergence of, say, Taylor series.Practically speaking, one can compute only a finite number of the b n 's in most situations.Truncating the expansion by taking the rest of the b n 's to be zero (or any constant) rarely provides a good approximation to the whole function [23].

B. Hydrodynamical Phenomenology
Long-time and large-wavelength properties of correlation functions are governed by emergent hydrodynamics.For each conserved charge (e.g.energy, spin), the density field should relax to equilibrium in a manner prescribed by a classical partial differential equation.Often this is a diffusion equation, though more exotic possibilities such as anomalous diffusion and ballistic transport (infinite conductivity) can also appear.
A numerical (and sometimes experimental) protocol to probe the emergent hydrodynamics is to study the autocorrelation function of the density wave operator O q = x e iqx Q x (here Q x is the operator of the conserved charge at x) at a range of momenta q.The behavior at large time is of especial interest, and can, in turn, be read off from the singularity structure of the Green's function.Let us give a few examples.If the closest pole to the origin is at z = iγ, then the autocorrelation function will decay exponentially as e −γt , while if the location of the closest pole varies quadratically as z = iDq 2 /2, then the dynamics are diffusive.However, the presence of non-linear terms in addition to the linear diffusive ones can give rise to exotic behavior where the diffusion constant itself becomes a function of frequency.An example of this is G(z) = z − iD(z)q 2 /2 −1 , where At any fixed q, G(z) has a branch cut in addition the diffusive pole, so although the diffusion constant D 0 is still well-defined, autocorrelation functions decay as a power law in time [70][71].Regardless, the full singularity structure of the Green's function determines the long-time behavior.Of course, computing the singularity structure of the Green's function is a demanding task.Even in integrable models, determining if the correct hydrodynamics is, say, diffusion or anomalous diffusion is non-trivial -let alone computing diffusion coefficients (see Refs [72][73][74] for recent developments).Indeed, accurately computing diffusion coefficients has been the goal of much recent numerical work [75][76][77].This difficulty reflected in the continued fraction expansion (47): the location of the poles change with each new fraction, so the full analytic structure of G(z) depends on all of the b n 's.
Knowing that the coefficients obey the universal form ( 9) is not enough, because even though the wavefunction is spreading out into the semi-infinite chain exponentially fast, we are given no guarantee about the wavefunction at the origin n = 0.For instance, the correlation functions C 1 (t) = sech(αt) and C 2 (t) = 1 + t 2 −γ [23] both correspond to Lanczos coefficients that obey the hypothesis But C 1 (t) decays exponentially while C 2 (t) decays as a power law, so clearly the hypothesis alone is insufficient to establish long-time behavior.The power law decay is nonetheless reflected in the Lanczos coefficients for C 2 (t), which have an alterating subleading tail.Precisely, they have the form b n = αn + γ + (−1) n f n where the f n 's are positive and decay to zero.Therefore determining the long-time tail of C(t) probably requires additional information about the subleading corrections to the hypothesis.

C. Numerical Diffusion Coefficients
Despite the complex behavior of autocorrelation functions in the time domain, there are situations where the hypothesis alone suffices to compute diffusion coefficients.In the case where the b n 's approach the universal form (9) especially quickly and regularly [78], we are able to make a meromorphic approximation to G(z).The idea is as follows.In the semi-infinite chain picture, we may hope to calculate the first few Lanczos coefficients exactly, so we may describe behavior near the origin n = 0 exactly.For large n, on the other hand, the hypothesis gives the coefficients almost exactly, so we can describes the dynamics by some exact solution.By stitching the dynamics at large and small n together, we can hope to find the dynamics on the whole chain.This allows us to recover a diffusive dispersion relation and numerically extract the diffusion constant in specific models.
We remark that there are a number of existent extrapolation schemes to determine the Green's function from the first few Lanczos coefficients [22,23].The new ingredient here is the hypothesis, which controls the approximation.
To make this idea into a precise numerical technique, we need three ingredients: a way to compute the Lanczos coefficients at small n, an exact solution at large n, and a robust way to meld them together.For a 1D spin chain in the thermodynamic limit of large system size, it is straightforward to compute the first few dozen Lanczos coefficients exactly through repeated matrix multiplication.Details are given in Appendix C.
To find the large n-behavior, we employ an exact solution for the quantum mechanics problem on the semiinfinite chain.If the hypothesis is obeyed, then the b n 's also asymptotically approach the form where η = 2γ/α + 1.The agreement is better, of course, at large n.The coefficients b n have the virtue that the quantum mechanics problem they describe on the semiinfinite chain is exactly solvable.Appendix D applies the theory of Meixner orthogonal polynomials of the second kind to determine the autocorrelation analytically: C(t) = sech(αt) η .(This is the same exact solution used in Section V above.)By Laplace transform, the corresponding Green's function is Here 1 F 2 is the hypergeometric function and M k depends on b k .It is crucial that G (n) (z) is known analytically, so that (49) provides the asymptotically exact large nbehavior.Now we stitch the small and large n information together.The true Green's function G (N ) (z) only depends on the coefficients b n with n ≥ N .So for sufficiently large N , where the b n 's are approximately the same as the b n 's, we may approximate an approximation that becomes better at large N .Equation ( 50) is our semi-analytical approximation to the Green's function.One can check that this is a meromorphic approximation for G(z), whose poles lie only in the upper half plane.In practice, one must calculate the b n 's until the universal behavior appears and fit α and η.Then the approximate G(z) can be calculated from (49) and a sequence of two-by-two matrix multiplications.One can then find the location of the first pole on the imaginary axis for a range of wavevectors q and fit z = iDq 2 /2+O(q 4 ) to extract the diffusion coefficient D. The diffusion constant for the chaotic Ising model is shown in Fig. 6 as a proof of concept.Almost all the computational effort goes into in computing the first few b n 's exactly.
In short, the hypothesis is sometimes sufficient to describe the emergent hydrodynamic behavior of operators.We reiterate that the hypothesis governs the leading order asymptotics of the Lanzcos coefficients only, while the autocorrelation depends on further corrections, so there is no a priori reason it should be computable just from the hypothesis.Subsequent work will provide further examples of this algorithm and discuss its theoretical and practical accuracy.

A. Discussion
We have presented a hypothesis on the universal growth of operators: the Lanczos coefficients follow the The Lanczos coefficients for q = 0.15 are fit to (48) with α = 0.35 and η = 1.74.We found it actually better not to approximate G (N ) (z) by G (N ) (z), but instead by G (N +δ) (z) for some integer offset δ so that η ≈ 1 (in the example shown, δ = 12).Large η or negative values lead to numerical pathologies.(b) The approximate Green's function (50) at q = 0.15.The arrow shows the "leading" pole that governs diffusion.(c) The locations of the leading poles for a range of q.One can clearly see the diffusive dispersion relation z = iDq 2 /2 + O(q 4 ).Fitting yields a diffusion coefficient D = 3.3 (5).
asymptotically linear form b n = αn + γ + o(1) in nonintegrable systems.We have seen copious evidence that the hypothesis is satisfied in a wide variety of nonintegrable models.Over the course of this work, the growth rate α has emerged as a quantity of prime importance, tying a diverse array of seemingly-disparate ideas together.Let us recount them now: • α > 0 is the slope of asymptotically linear growth of the Lanczos coefficients.
• ±iπ/(2α) are the locations of singularities closest to the origin in the (analytic continuation) of the autocorrelation C(t), see Appendix A.
• 2α is the exponential growth rate of Krylovcomplexity.
• 2α is an upper bound for the growth of all qcomplexities.
• 2α is an upper bound for the Lyapunov exponent (whenever the latter is well-defined), since quantum OTOCs are an example of q-complexities.
We have, of course, put aside the precise conditions and qualifiers of each statement.In light of these results, α plays a central role in operator growth and dynamics of complex systems.
Complexity -especially the Krylov-complexityarose as a key concept in this work.We would like to highlight its temporal nature which, as we now argue, makes it a more general notion than chaos.Chaos essentially tracks the development of structures at ever-smaller scales in phase space.In classical systems, of course, this may proceed indefinitely, while in quantum systems, features smaller than are ruled out and the process saturates.Chaos therefore cannot carry over straightforwardly to systems deep in the quantum regime, where the phase space volume is comparable to and saturation occurs almost immediately.The K-complexity, in sharp contrast, measures structures at ever-smaller scales in the time domain.We believe this is a fundamental difference; as we have seen, the K-complexity can grow exponentially in quantum systems beyond semiclassical or large-N limits.Operator complexity may well supersede the notion of chaos in quantum dynamics.

B. Outlook
We would like to understand how our hypothesis can be affected by obstructions to thermalization.In free and integrable models, there are an extensive number of conserved local or quasi-local quantities.The behavior of the Lanczos coefficients in integrable models is likely non-universal, and depends strongly on the model and operator in question [23].In the integrable case, it may be appropriate to modify the Lanczos algorithm to promote the semi-infinite line to a lattice where the perpendicular direction is generated by commutators against quasi-local conserved charges.Another exceptional case is quantum scar states [79][80][81], isolated states that fail to thermalize in otherwise chaotic systems, possibly due to emergent or approximately conserved charges.It would be revealing to see how scars are reflected in the Lanczos coefficients.Finally it would be of great interest to understand the interplay of the hypothesis with many-body localized systems (see [82] and references therein) where thermalization fails.
An open question is how the hypothesis may be extended to general finite-temperature systems.The key physical question is what operator norm (30) to choose.Should the linear growth prevail when T < ∞, the value of α would depend on the choice of norm.This opens the possibility that different norms should be chosen for different questions.For instance, the choice (31) for comparing with the universal bound on Lyapunov exponent may not be appropriate when calculating dynamical response with the Kubo formula.Any choice of norm at finite temperature presents a numerical challenge due to the presence of the thermal density matrix [22,83,84].Quantum Monte Carlo seems promising for this problem, as the Lanczos coefficients can be computed without analytic continuation.In low dimensions, DMRG can be also useful: matrix product operators can be used to approximate the thermal state, and the operators in the Lanczos algorithm.
One would like to put the hypothesis on more solid mathematical footing.We speculate that this may be achieved within an extended random matrix theory.Standard proofs of the Wigner semicircle law exploit the connections between the moments of a distribution, the combinatorics of Dyck paths, Catalan numbers, and the Stieltjes transform of a distribution [85].These are directly analogous to the moments µ 2n , the combinatorics of Motzkin paths, secant numbers, and the continued fraction expansion for G(z) -all of which arose in the calculation of our exact wavefunction in Appendix D).The non-trivial appearance of the same type of objects in both contexts suggests a strong analogy.We thus conjecture that the hypothesis can be derived analytically by introducing new type of random matrix ensemble that incorporates locality and translation invariance.(This is similar to the framework of [86].)In this case, a Hamiltonian such as H = i h i,i+1 , where h i,i+1 is a random matrix, should obey the hypothesis (9) in expectation.Therefore generic, 2-local Hamiltonians would also be expected to obey the hypothesis by concentration of measure.It may well be that showing the hypothesis holds for a specific Hamiltonian is of comparable difficulty to showing the ergodic hypothesis applies to specific classical systems.
To close, we wish to point out that the territory of qcomplexities beyond K-complexity and OTOCs is completely unexplored.In generic many-body systems (i.e.not semiclassical) at infinite temperature, these two examples represent two extremes, showing maximal and non-existent exponential growth rates, respectively.The significant gap between them should be filled with potentially more meaningful measures of complexity.These complexities could be entirely new concepts or disguised forms of existing notions such as entanglement or circuit complexity.Hopefully, charting this terra incognita will continue to shed new light on the complex nature of many-body quantum dynamics.
The first two moments µ 2 and µ 4 are (twice) the norm squared of the Lγ 1 and L 2 γ 1 , respectively.Under disorder averaging, the terms on the right-hand side are orthogonal, and each corresponds to a different diagram: The combinatorial factor is trivial for each of the above graphs.The SYK diagrams encode the Schwinger-Dyson equations governing the autocorrelation and Green's function which are, up to trivial transformations, the exponential and ordinary generating functions of the moments, respectively: that is, Σ(z) and Σ(t) are related by (non-standard) Laplace transform (45) just as G(z) and C(t) are.Equation (B7) can be represented diagrammatically (here for the case q = 4) by The dot represents a general SYK diagram (a fullydressed Green's function).This is the sum of the bare Green's function, or the time-domain product of (q − 1) dressed Green's functions.Note that both exponential and ordinary generating functions are needed to take the combinatorial factors into account: a serial (respectively, parallel) composition of diagrams correspond to product of ordinary (resp.exponential) generating function.Equation (B7) has no closed form solution for general q.However, working with the power series representations, it enables the numerical calculation of µ 2 , . . ., µ 2n in polynomial time and space complexity in n.Concretely, the following iteration algorithm can be easily implemented in a computer algebra system: 1. Set g 0 (z) := z −1 , and let j = 0.
6. Increment j by 1 and repeat from step 2.
When the above procedure is stopped at j = n, the result g n (z) will be a polynomial truncation of the Green function: g n (z) = n j=0 µ 2j z −2j−1 , which contains the correct moments up to µ 2n .They can be then used to compute Lanczos coefficients b 2 1 , . . ., b 2 n by the recipe (A3).Arbitrary-precision rational number arithmetic is necessary for n ∼ 10 2 , since the moments grow very fast.We calculated b n for a few different values of q up to n = 100, and extracted the linear slope by a linear fit.The results are reported in Table I and Fig. 2 (a).
In the large-q limit, (B7) can be solved analytically.It is convenient to define the coupling constant [13,34] J 2 := 2 1−q q J 2 . (B9) It is then known [13,34] that C(t) admits a 1/q expansion where the leading non-trivial term satisfies the following differential equation: whose solution is The corresponding moments where (T n ) ∞ n=0 = (1, 2, 16, 272, 7936, . . . ) are the tangent numbers [89].The generating function of T n admits a and complex coefficients.The Liouvillian is applied by combining (C2), (C3), and (C5).Of course, it is not necessary to take O to be translation invariant.One could equally well take a small single-site operator and apply the same technique without the sum over all sites.We note that Lanczos algorithm (4) only requires the storage of three operators at any time.In practice the method described here allows a few dozen Lanczos coefficients to be computed in a few minutes on a modern laptop and is generally memory-limited by the exponential increase in the number of Pauli strings required.
Once the Lanczos coefficients and Krylov vectors have been computed, it is possible to understand how the operators O n grow in physical space.One way to characterize this is in terms of the distribution of string lengths in each O n .If O n = a c a σ a , where the sum runs over all Pauli string a, then the distribution is defined by P n (s) = a:|a|=s |c a | 2 .This distribution is shown for the Hamiltonian H 1 with the parameters given in Fig. 2. The mean and variance of the distribution grow with n.We have observed that the distribution P n (s) appears to be highly model-dependent.This makes it difficult to translate information about the exponential spreading of the wavefunction in the semi-infinite chain back to physical space.

Appendix D: A Family of Exact Solution with Linear Growth
This section will provide a derivation for the exact solution (17) of the 1d quantum mechanics problem with Lanczos coefficients To solve this problem, notice that our infinite, tridiagonal matrix is actually quite a familiar setup.If instead we had b n = √ n, then L would be the matrix representing the Hamiltonian for the quantum harmonic oscillator in the basis of raising and lowering operators.So really this is just a 1d quantum mechanics problem, albeit not a standard one.In particular, it is known that system described by L has very high symmetry, due to an infinite-dimensional representation of the Lie algebra su(1, 1), enabling us to find an exact solution [91,92].Indeed, there is a rich mathematical literature on the close connections between representations of su(1, 1), the combinatorics of Motzkin paths, and Meixner orthogonal polynomials [93,94].Our solution will be a simple application of these mathematical results.
We start with some generalities on orthogonal polynomials.Define L (n) = L 0≤i<n,0≤j<n to be the n × n matrix in the upper-left block of L. For example, We then define polynomials for each n via By performing a cofactor expansion for the determinant on the nth row, the Q's admit a three-term recursion relation together with initial conditions Q 0 (z) = 1 and Q −1 (z) = 0. Eq. (D4) should be compared with where {e n } is the natural orthonormal basis of L. In fact, (D4) and (D5) are equivalent, under the identification: Therefore, the polynomials Q n (z) are orthogonal, but not normalized.Instead they are monic, i.e., the highest order coefficient is 1: . By construction, both {Q k (z)} and {z n } are a basis of C[z] and can be related by a triangular linear transform with matrix elements µ n,k : Combined with (D6), and by orthonormality of {e n }, we have and therefore The statements so far are general and apply to any set of Lanczos coefficients.
In the specific case b n = n(n − 1 + η); the extra overall factor α in (D1) can be later recovered by a simple time rescaling;, one may recognize from the recursion relation (D4) that Q n 's are a special case of the Meixner polynomials of the second kind [95].They are a nonclassical family of orthogonal polynomials defined by the following three-term recursion: [96,97] In particular, Q n (z) = M n (z; δ = 0, η).For these polynomials, the matrix elements µ n,d have been exactly calculated, in terms of the following generating function [94]: As a side note, we mention that the above generating function, referred to as that of the "inverse polynomials" in the theory of orthogonal polynomial, is closely related to the generating function of Meixner polynomials themselves.The latter has also a closed form expression, known to be of Sheffer type [93,96]: Now, taking δ = 0 and the series coefficient of w d in (D11), we have Applying this to (D9), and recalling b n = n(n − 1 + η), we obtain the wavefunction solution where (η is the Pochhammer symbol.The general solution for b n = α n(n − 1 + η) can be obtained by a simple rescaling t → αt, and is precisely Eq. ( 17) of the main text where, of course, O n |e iLt |O 0 = e n |e iLt |e 0 .The special case η = 1 of this family of solutions is well-known [23,29].To the best of our knowledge, the general solution (D13) has not been applied to the recursion method.
Appendix E: Derivation of the q-Complexity Bound This Appendix will derive Eq. ( 25), (Q) t ≤ C (n) t for C = 2d.The main idea of is that the definition of Q guarantees that the eigenbasis of Q is dilated by a factor of at most C compared to the Krylov basis.
We first show that the Krylov basis vectors have a bounded number of components in the Q basis due to the dilation property.For any operator Φ where there is an M > 0 such that (q a |Φ) = 0 for q a > M , the hypothesis (20b) implies that (q a |L|Φ) = 0 for q a > M + d.Using (20c), as a base case for induction, we have (q a |L n |O) = 0 for q a > d(n + 1) and, in particular, for q a > Cn.By the construction of the Krylov basis, (q a |O n ) = 0 if q a > Cn. (E1) We claim that (E1) implies for any operator wavefunction Φ; taking Φ = O(t), we obtain (25).
To show (E2), we introduce projectors to large spectral values in the Krylov and Q bases, respectively: Then, we have for n = q/C, because m < n = q/C ≤ q a /C, (q a |O m ) = 0 by (E1).Equivalently, Applying this equation and its Hermitian conjugate, we have where the inequality follows from the fact that P Q q is a projector.Finally we need a standard integration-byparts identity that converts the expectation value to an integral over the projectors: for any k = 1, 2, 3, . . . .Combining the case k = 1 and (E5), we obtain which finishes the proof.More generally, for any k, we have This is useful as a bound on the growth rate of higher moments of the q-complexity super-operator.See Section VI C for an application.
Appendix F: Microscopic Origin of the Hypothesis While the evidence presented in Section IV suffices alone to support the hypothesis, we find it instructive to discuss its microscopic origin.We will first re-iterate why linear growth is maximal, comment on the importance of interference, and finally show how our hypothesis can be understood as an assumption on the form of the spectral function which arises in ETH.Though our arguments here will not be truly rigorous, we aspire to higher precision in this section.
To this end, we consider the specific setting of a translation-invariant spin chain.The restriction to a single dimension is a wholly artificial one to simplify notation.The arguments in this section extend to any dimension at the cost of greater book-keeping.Suppose we have an r-local Hamiltonian H = x h x and an r-local operator O.The Liouvillian becomes a sum of terms L = x x with x = [h x , •] and we note that O k is r + k(r − 1)-local.Finally, we suppose that the local bandwidth ||h x || ≤ E is finite.
The main tool in our analysis will be the growth of the moments µ 2n = O|L Hence linearly increasing Lanczos coefficients are equivalent to factorial-squared growth of the moments.Exponential corrections to the moments amount, therefore, to a change in the value of α, but are negligible when checking if the hypothesis holds.We have already seen through two indirect arguments that (9) is the fastest possible growth rate of the Lanczos coefficients.Now we show this fact directly using the moments (this is essentially the same method used in [37,38]).The moment µ 2n is the norm-squared of the sum This sum is highly constrained by the spatial structure of the spin chain.The operator O is supported only on a few sites, and each application of the Liouvillian grows that support at the edges.Each term in (F2) can be visualized as a discrete quantum circuit, where each gate Hence the moments grow no faster than (F1), so we have (again) shown that ( 9) is the maximal growth rate.It remains to be seen why linear growth ( 9) is achieved in generic non-integrable systems.It is important to keep the quantum nature of the problem in mind; each term in (F2) has both a value and a complex phase, so terms interfere either constructively or destructively.This interference is paramount to achieving linear growth.To see why, let us examine a fallacious hypothetical: each term in the sum (F2) is uncorrelated with the rest and their phases are random and independent.In this case, a simple estimate shows µ 2n ∼ n! and, hence, b n ∼ √ n.However, this assumption is mistaken, because large classes of terms interfere constructively.For instance, whenever x k and x k+1 are separated by more than a distance r they commute and thus xn Such identical terms are commonplace and occur whenever two of the x 's are outside each other's "lightcone".The presence of such classes amounts to a conspiracy of correlations in the sum (F2), which enables the moments to exceed the uncorrelated estimate n! qualitatively.Indeed, even if the value of the sum is reduced exponentially by destructive interference, the moments will realize factorial-squared growth µ 2n ∼ (n!) 2 .(This also depends on having sufficiently many non-zero terms in the sum, an assumption that breaks down in generic non-integrable systems.)In other words, any small bias in favor of constructive rather than destructive interference is sufficient to realize the hypothesis.
While it is not easy to precisely estimate the small bias, we know that it is reflected in the extensiveness of energy in short-range systems.To see explicitly how extensiveness is related to the hypothesis, let us "derive" it again, assuming the off-diagonal eigenstate thermalization hypothesis (ETH) [3][4][5].We stress that our hypothesis does not require ETH, nor is it implied by ETH alone; what follows is equivalent to saying the exponential decay of the spectral function implies the hypothesis as (15).In the energy eigenbasis, applying the Liouvillian is easy, so We wish to show this has factorial-squared (n 2n ) dependence.The extensiveness of energy says that E a −E b ∝ n for some proportion of the eigenstates, so it is sufficient to show those contributions to the sum are suppressed no more than exponentially by the other terms.
Since n is large, we may approximate the sum by an integral The integral may be carried out in the saddle-point approximation to find µ 2n ∼ n 2n e −nSmin = n 2n e O(n) , in agreement with (F1).Note that due to the logarithmic term, the saddle point must be away from the diagonal = .In other words, the dominant contribution to the linear growth of the b n 's comes from matrix elements O ab with extensive energy differences E a − E b ∝ n.

FIG. 4 .
FIG.4.The exact solution wavefunction(17) in the semiinfinite chain at various times.The wavefunction is defined only at n = 0, 1, 2 . . ., but has been extrapolated to intermediate values for display.

Example 1 :
K-complexity.The K-complexity is always a q-complexity, with Q = n.The basis |q a ) is just the Krylov basis |O n ) and the conditions (20b) and (20c) are satisfied by construction of the Krylov basis with d = 1.

2 FIG. 5 .
FIG. 5. (a) The growth rate α versus the classical Lyapunov exponent λL/2 in the classical Feingold-Peres model of coupled tops, (43).α ≥ in general, with equality around the c = 0 where the model is the most chaotic.The growth rate appears to be discontinuous at the non-interacting points c = ±1, similarly to Fig. 2-(b).(b) The first 40 Lanczos coefficients of quantum s = 2, . . ., 32 and classical (s = ∞) FP model, with c = 0.

FIG. 6 .
FIG. 6. Numerical computation of the diffusion coefficient for the energy density operatorO = Eq in H = i XiXi+1 − 1.05Zi + 0.5Xi.(a)The Lanczos coefficients for q = 0.15 are fit to(48) with α = 0.35 and η = 1.74.We found it actually better not to approximate G (N ) (z) by G (N ) (z), but instead by G (N +δ) (z) for some integer offset δ so that η ≈ 1 (in the example shown, δ = 12).Large η or negative values lead to numerical pathologies.(b) The approximate Green's function(50) at q = 0.15.The arrow shows the "leading" pole that governs diffusion.(c) The locations of the leading poles for a range of q.One can clearly see the diffusive dispersion relation z = iDq 2 /2 + O(q 4 ).Fitting yields a diffusion coefficient D = 3.3(5).
C G = 1.The first non-trivial combinatorial factor C G = 6 for the diagram , which contributes to µ 6 .The six vertex orderings are

FIG. 7 .
FIG.7.The size distribution of the Pauli strings in the Krylov vectors On for the Hamiltonian H1 with parameters and initial operator as in Fig.2.Though the distribution drops quickly after its peak, Pn(s) is supported on [0, n/2 + 2].

x
k+1 must act on at least one site that is already in the support of x k • • • x1 O -otherwise the term vanishes due to the commutator.This condition is satisfied by at most r + k(r − 1) positions x k , so the total number of non-zero terms in (F2) is at most n!r n for large n.The value of each non-zero term is itself bounded due to the finite bandwidth E, so|| xn • • • x1 O|| 2 ≤ (2E) 2n.By the triangle inequality,µ 2n = ||L n O|| 2 ≤ (n!) 2 r 2n (2E) 2n .(F3) Fix a large n and consider a finite subsystem S of size O(n), large enough to contain the support of each O k for 0 ≤ k ≤ n.We work in the basis of energy eigenstates of H (restricted to S), |E a for 1 ≤ a ≤ D, where D ∼ exp(O(n)) is the Hilbert space dimension.The Liouvillian has normalized energy eigenoperatorsO ab = D 1/2 |E a E b | (F4)with eigenvalue E a − E b .The off-diagonal ETH states that we can write the local operator O as(O ab |O) = D −1/2 E b |O|E a = D −1 f (E, ω)R ab , a = b(F5) where R ab are independent random variables with mean zero and unit variance,E = (E a + E b )/2, ω = E a − E b ,and f (E, ω) is the ETH spectral function.The spectral function is a smooth function of its arguments and is exponentially suppressed at large |ω|: |f |2 ∼ e −|ω|/ω0 .In fact, |f (E = 0, ω)| 2 is equal to the spectral function Φ(ω) defined in Eq. (12).