A Universal Operator Growth Hypothesis

TL;DR

This paper proposes a universal hypothesis linking operator growth rates in many-body systems to exponential complexity growth, providing bounds on chaos measures and new methods for computing diffusion constants.

Contribution

It introduces a universal operator growth hypothesis with a growth rate  that bounds chaos and complexity measures, applicable across various models and systems.

Findings

Lanczos coefficients grow linearly with rate  in generic systems

The growth rate  bounds Lyapunov exponents as   ; a universal bound

The hypothesis is validated in spin chains, SYK model, and classical models.

Abstract

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, with an extra logarithmic correction in 1d. 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 obtain a sharp bound on Lyapunov exponents $\lambda_L \leq 2 \alpha$, which complements and improves 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.

Figures

Tables

Table 1. Table 1: The growth rate α 𝛼 \alpha versus half the OTOC-Lyapunov exponent λ L / 2 subscript 𝜆 𝐿 2 \lambda_{L}/2 in the q 𝑞 q -SYK model ( 18 ) in units of 𝒥 = q ​ 2 ( 1 − q ) / 2 ​ J 𝒥 𝑞 superscript 2 1 𝑞 2 𝐽 \mathcal{J}=\sqrt{q}2^{(1-q)/2}J . Here α 𝛼 \alpha is obtained by exact numerical methods discussed in Appendix  B , while λ L subscript 𝜆 𝐿 \lambda_{L} is taken from the Appendix of  Roberts  et al. ( 2018 ) . The q 𝑞 q -SYK is physical only for even integers q 𝑞 q , but large- N 𝑁 N methods allow an extrapolation to any q ≥ 2 𝑞 2 q\geq 2 .

$q$	2	3	4	7	10	$\infty$
$\alpha /𝒥$	0	0.461	0.623	0.800	0.863	1
${\lambda }_L/\left(2𝒥\right)$	0	0.454	0.620	0.799	0.863	1

Table 2. Table 2: Correspondence between finite and infinite temperature definitions and results.

		$T=\infty$
Inner Product		$\left(A|B\right)\propto \mathrm{Tr}\left[A^{†}B\right]$	Eq. (53)
Lanczos Algorithm		Eq. (4)	Eq. (55)
$C\left(t\right),G\left(z\right),\Phi \left(\omega \right),{\mu }_{2n}$		Section III	Eq. (56)
$b_n\leftrightarrow C\leftrightarrow G\leftrightarrow \Phi \leftrightarrow \mu$		App. A	App. A
Hypothesis		Eq. (22)	Eq. (58)
$b_n\sim \alpha n$ for SYK		Eq. (93)	Eq. (99)
Bound ${\lambda }_L\le 2\alpha$		Proven	Conjectured

Equations

$$(n{)}_t\sim e^{2\alpha t}.$$

$${\lambda }_L\le 2\alpha .$$

$$C(t)=\mathrm{Tr}[\mathcal{O}(0)\mathcal{O}(t)]/\mathrm{Tr}[1]=\left(\mathcal{O}|\exp \left(i\mathcal{L}t\right)|\mathcal{O}\right),$$

$$|A_n)$$

$$b_n$$

$$|{\mathcal{O}}_n)$$

$$L_{nm}:=\left({\mathcal{O}}_n|\mathcal{L}|{\mathcal{O}}_m\right)=\left(\begin{array}{ccccc}0& b_1& 0& 0& \cdots \\ b_1& 0& b_2& 0& \cdots \\ 0& b_2& 0& b_3& \cdots \\ 0& 0& b_3& 0& \ddots \\ ⋮& ⋮& ⋮& \ddots & \ddots \end{array}\right).$$

$$\mu_{2n}:=\mathinner{\left(\mathcal{O}|\mathcal{L}^{2n}|\mathcal{O}\right)}=\frac{d^{2n}}{dt^{2n}}C(t)\big{|}_{t=0}$$

$$C(t)={\left(e^{iLt}\right)}_{00}.$$

$${\partial }_t{\phi }_n=-b_{n+1}{\phi }_{n+1}+b_n{\phi }_{n-1},{\phi }_n(0)={\delta }_{n0}.$$

$$G(z)=\left(\mathcal{O}\left|\frac{1}{z-\mathcal{L}}\right|\mathcal{O}\right)=i{\int }_0^{\infty }{e}^{-izt}C(t)dt.$$

$$\Phi (\omega )$$

$$C(t)\leftrightarrow G(z)\leftrightarrow \Phi (\omega )\leftrightarrow \{{\mu }_{2n}\}\leftrightarrow \{b_n\}$$

$$b_n=\alpha n+\gamma +\mathrm{o}(1),$$

$$b_n\sim n^{\delta }⟺\Phi (\omega )\sim \exp (-|\omega /{\omega }_0{|}^{1/\delta })$$

$$\Phi (\omega )\le C{e}^{-\kappa |\omega |},\kappa =\frac{1}{2e{G}_r||h_i||}$$

$$b_n$$

$$\Phi (\omega )$$

$$\alpha \le \pi e{G}_r||h_i||,$$

$$t=\pm \frac{i\pi }{2\alpha }.$$

$$H_{\text{SYK}}^{(q)}=i^{q/2}\sum _{1\le i_1<i_2<\cdots <i_q\le N}{J}_{i_1\dots i_q}{\gamma }_{i_1}{\gamma }_{i_2}\cdots {\gamma }_{i_q}$$

$$b_n^{\text{SYK}}=\{\begin{array}{ll}\mathcal{J}\sqrt{2/q}+\mathrm{O}(1/q)& n=1\\ \mathcal{J}\sqrt{n(n-1)}+\mathrm{O}(1/q)& n>1,\end{array}$$

$$H=\sum _{x,y}{X}_{x,y}{Z}_{x+1,y}+Z_{x,y}{X}_{x,y+1}$$

$${\mu }_{2n}=n^{2n}{e}^{O(n)}$$

$$b_n=\{\begin{array}{ll}A\frac{n}{W(n)}+\mathrm{o}(n/\ln n)\sim A\frac{n}{\ln n}+\mathrm{o}(n/\ln n)& d=1\\ \alpha n+\gamma +\mathrm{o}(1)& d>1\end{array}$$

$$A\le \frac{4}{e}||h_x||.$$

$${\partial }_t{\phi }_n=-b_{n+1}{\phi }_{n+1}+b_n{\phi }_{n-1},{\phi }_n(0)={\delta }_{n0}.$$

$${\tilde{b}}_n:=\alpha \sqrt{n(n-1+\eta )}\stackrel{n\gg 1}{\to }\alpha n+\gamma ,$$

$$|\mathcal{O}(t))=\sum _{n=0}^{\infty }\sqrt{\frac{{\left(\eta \right)}_n}{n!}}\tanh (\alpha t{)}^n\mathrm{sech}(\alpha t{)}^{\eta }{i}^n|{\mathcal{O}}_n)$$

$${\left(n\right)}_t:=\left(\mathcal{O}(t)|n|\mathcal{O}(t)\right)=\eta \sinh (\alpha t{)}^2\sim e^{2\alpha t},$$

$$-i{\dot{\mathcal{O}}}_0(t)$$

$$-i{\dot{\mathcal{O}}}_1(t)$$

Full text

1. I Introduction

2. II Synopsis

3. III Preliminaries: The Recursion Method

4. IV The Hypothesis

5. IV.1 Upper Bounds

6. IV.2 Analytical Evidence

7. IV.3 The Special Case $d=1$

8. V Exponential Growth of Complexities

9. V.1 Exponential Growth in the Semi-infinite Chain

10. V.2 A Bound on Complexity Growth

11. VI Growth Rate as a Bound on Chaos

12. VI.1 SYK Model

13. VI.2 Classical Chaos

14. VI.2.1 A (Weaker) Bound on Classical Chaos

15. VI.2.2 Numerical Case Study

16. VII Application to Hydrodynamics

17. VII.1 Continued Fraction Expansion: Brief Review

18. VII.2 Hydrodynamical Phenomenology

19. VII.3 Numerical Diffusion Coefficients

20. VIII Finite Temperature

21. VIII.1 Choice of Inner Product

22. VIII.2 Bound on Chaos

23. VIII.3 SYK Model

24. IX Conclusions

25. IX.1 Discussion

26. IX.2 Outlook

27. A Brief Review of the Recursion Method

28. A.1 From Moments to Lanczos Coefficients

29. A.2 From Lanczos Coefficients to Moments

30. B Moments and Lanczos Coefficients in the SYK Model

31. B.1 General Method

32. B.2 Large-$q$ limit

33. C Numerical Details for 1d Spin Chains

34. D A Family of Exact Solution with Linear Growth

35. E Derivation of the q-Complexity Bound

36. F Geometric Origin of the Upper Bounds

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 Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016); Deutsch (2018).

Recent work has shifted focus from states to operator growth in many-body systems von Keyserlingk et al. (2018); Nahum et al. (2018); Rakovszky et al. (2017); Khemani et al. (2018); Gopalakrishnan et al. (2018); Chan et al. (2018). Under Heisenberg-picture 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 von Keyserlingk et al. (2018); Nahum et al. (2018), and extended to increasingly realistic systems involving conservation laws Rakovszky et al. (2017); Khemani et al. (2018), Floquet dynamics Chan et al. (2018), and even interacting integrable models Gopalakrishnan et al. (2018).

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 out-of-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 Maldacena et al. (2016) 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 Maldacena and Stanford (2016); Kitaev (2015); Sachdev and Ye (1993). Unfortunately, the OTOC does not necessarily exhibit exponential growth outside of semiclassical or large-$N$ limits, rendering the Lyapunov exponent ill-defined Fine et al. (2014); Xu and Swingle (2018a); Khemani et al. (2018); Xu and Swingle (2018b). A general theory of operator growth under generic, non-integrable 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 well-understood technique, dating back to the 1980s Mattis (1981), and has recently been used to compute conductivities in strongly-interacting systems Lindner and Auerbach (2010); Khait et al. (2016); Auerbach (2018). It is surveyed in great detail in Ref. Viswanath and Müller (2008). Essentially, it converts any linear-response calculation to the problem of a quantum particle on a half chain, with the hopping matrix elements given by the Lanczos coefficients $b_{n}$. Section IV presents our hypothesis: operators in generic, non-integrable systems have Lanczos coefficients with asymptotically linear growth with $n$, suppressed by a logarithmic correction in one dimension. The linear growth rate, denoted $\alpha$, is the central quantity of this work. It has dimensions of energy and can be bounded from above by the local bandwidth [see (16) and (23)]. In light of this, the hypothesis essentially asserts that the Lanczos coefficients grow as fast as possible in non-integrable systems. 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 experimentally McArthur et al. (1969); Engelsberg and Lowe (1974); Lundin et al. (1990).

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}\sim\alpha n$, captures the irreversible process of simple operators evolving into complex ones. Furthermore, the 1d wavefunction delocalizes exponentially fast on the $n$ axis, at a rate exactly given by $\alpha$. Asymptotically, the expected position of the 1d wavefunction satisfies

[TABLE]

The expectation value $\mathinner{\left(n\right)}_{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.2. Crucially, this class includes out-of-time-order correlators. This allows us to establish a quantitative connection between $\alpha$ and the Lyapunov exponent, which will be the subject of Section VI. We show for quantum systems at infinite temperature that the growth rate gives an upper bound for the Lyapunov exponent whenever the latter is well-defined:

[TABLE]

For classical systems, this statement is a conjecture but it is posible to prove a somewhat weaker bound. We check (2) in the SYK model and a classical tops model, and find it 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.

Section VIII discusses the generalization to finite temperatures, which involves many open questions. Nevertheless, we show that the universal bound on chaos $\lambda_{L}\leq 2\pi k_{B}T/\hbar$ Maldacena et al. (2016) can be implied and improved by a proper finite-temperature extension of the bound (2), and provide evidence supporting this conjecture. We conclude in Section IX 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 Viswanath and Müller (2008). Consider a local Hamiltonian $H$ and fix a Hermitian operator $\mathcal{O}$. We regard the operator as a state $\mathinner{|\mathcal{O})}$ in the Hilbert space of operators, endowed with the infinite-temperature inner product $\mathinner{\left(\mathcal{O}_{1}|\mathcal{O}_{2}\right)}:=\operatorname{Tr}[\mathcal{O}_{1}^{\dagger}\mathcal{O}_{2}]/\operatorname{Tr}[1]$. We write $\left|\left|\mathcal{O}\right|\right|:=\mathinner{\left(\mathcal{O}|\mathcal{O}\right)}^{1/2}$ 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 $\mathcal{L}:=[H,\cdot]$. Our central object is the autocorrelation function

[TABLE]

where the second equality follows from Baker-Campbell-Hausdorff.

Computing $C(t)$ is inherently difficult. Suppose $\mathcal{O}(t=0)$ is a relatively simple operator that can be written as the sum of a few basis vectors in any local basis 111A local basis in, say, a spin chain is any basis related to the basis of Pauli strings by a finite-depth local unitary circuit.. As the spatial support of $\mathcal{O}(t)$ grows, the number of non-zero coefficients of $\mathcal{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.

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 $\left\{\mathcal{L}^{n}\mathinner{|\mathcal{O})}\right\}$, and then apply Gram-Schmidt to orthogonalize. Explicitly, start with a normalized vector $\mathinner{|\mathcal{O}_{0})}:=\mathinner{|\mathcal{O})}$. As a base case, let $\mathinner{|\mathcal{O}_{1})}:=b_{1}^{-1}\mathcal{L}\mathinner{|O_{0})}$ where $b_{1}:=\mathinner{\left(\mathcal{O}_{0}\mathcal{L}|\mathcal{L}\mathcal{O}_{0}\right)}^{1/2}$. Then inductively define

[TABLE]

The output of the algorithm is a sequence of positive numbers, $\left\{b_{n}\right\}$, called the Lanczos coefficients, and an orthonormal sequence of operators, $\left\{\mathinner{|\mathcal{O}_{n})}\right\}$, called the Krylov basis. (This is a bit of a misnomer, as the Krylov basis spans an operator space containing $\mathcal{O}(t)$ for any $t$, but does not usually span the full space of operators). The Liouvillian is tridiagonal in this basis:

[TABLE]

We make four remarks. First, if the operator Hilbert space is $d$-dimensional with $d$ finite (or if the subspace spanned by $\mathinner{|\mathcal{O}_{0})},\mathinner{|\mathcal{O}_{1})},\mathinner{|\mathcal{O}_{2})},\dots$ is so), the algorithm will halt at $n=d+1$: in this work, we work always in the thermodynamic limit and discard this non-generic 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}\mathcal{O}_{n}$ is Hermitian for all $n$, hence $\mathinner{\left(\mathcal{O}_{n}|\mathcal{L}|\mathcal{O}_{n}\right)}=0$. Third, the knowledge of the Lanczos coefficients $b_{1},\dots,b_{n}$ is equivalent to that of the moments $\mu_{2},\mu_{4},\dots,\mu_{2n}$, defined as the Taylor series coefficients of the correlation function

[TABLE]

The non-trivial transformation between the Lanczos coefficients and the moments is reviewed in Appendix A. Fourth, the Lanczos coefficients have units of energy.

In the Krylov basis, the correlation function $C(t)$ is:

[TABLE]

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 $\varphi_{n}(t):=i^{-n}\mathinner{\left(\mathcal{O}_{n}|\mathcal{O}(t)\right)}$. Heisenberg evolution of $\mathcal{O}(t)$ becomes a discrete Schrödinger equation:

[TABLE]

where $b_{0}=\varphi_{-1}=0$ by convention. The autocorrelation is simply $C(t)=\varphi_{0}(t)$, so the Lanczos coefficients are completely equivalent to the autocorrelation function.

Just as different bases are well-suited for particular computations, a number of equivalent representions of the autocorrelation function appear in this work, namely the Green’s function

[TABLE]

and the spectral function

[TABLE]

In summary, we have reviewed five equivalent ways to describe the dynamics

[TABLE]

Just as with a choice of basis, we shall use the most convenient representation for the task at hand and translate freely between them. We note that $\{b_{n}\}$ is special in the sense that it is a non-linear representation of the autocorrelation while all other representations are linearly related. We provide the details on the mapping to $b_{n}$ in Appendix A, with a particular focus on asymptotic properties.

The nonlinearity involved in $\{b_{n}\}$ also makes them more abstract. Intuitively, we can think of the Krylov basis $\{\mathcal{O}_{n}\}$ as stratifying operators by their ‘complexity’ (with respect to the initial operator $\mathcal{O}$), and $b_{n}$’s describe how operators of different complexities transform into one another. 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. Informally, in a chaotic quantum system, the Lanczos coefficients $\{b_{n}\}$ should grow as fast as possible. The maximal possible growth rate turns out to be linear (with logarithm corrections in 1D). Our precise statement is therefore as follows. Suppose that $H$ describes an infinite, non-integrable 222As a working definition, we say that a system is integrable if it has an extensive number of quasi-local conserved quantities., many-body system in dimension $d>1$ and $\mathcal{O}$ is a local operator having zero overlap with any conserved quantity (in particular, $\mathinner{\left(\mathcal{O}|H\right)}=0$). Then the Lanczos coefficients are asymptotically linear:

[TABLE]

for some real constants $\alpha>0$ and $\gamma$. This linear growth is an example of universality. We will refer to $\alpha$ 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.2). On the other hand, it is observable by standard linear response measures (Section IV.1). This section first describes why linear growth is maximal, amasses a weight of evidence in favor of the hypothesis, and finally discusses the special case of one dimension.

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 Viswanath and Müller (2008). In particular, it is known that non-interacting models (such as lattice free fermions) give rise to a bounded sequence $b_{n}\sim{}\mathrm{O}(1)$. If we start with a two-body operator $\mathcal{O}$ in such free models, all $\mathcal{O}_{n}$’s will remain two-body. In this sense, the operator dynamics is simple. In this work, we focus on the opposite extreme of generic chaotic dynamics. To our knowledge, the ubiquity of asymptotically linear growth in these systems and its consequences have not been systematically studied in quantum systems. Interacting models with obstructions to thermalization (e.g., integrable systems) lead to more involved behaviors, which have not been thoroughly explored. Nevertheless, a square root behavior $b_{n}\sim\sqrt{n}$ has been observed in a few examples  (Viswanath and Müller (2008); Lee (2001), see also Fig. 2).

IV.1 Upper Bounds

We start by showing that linear growth is the maximal possible growth of the Lanczos coefficients. This is most easily done starting with the spectral function. 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 $\Phi(\omega)$. The precise asymptotic behavior is Lubinsky (1987); Magnus (2012)

[TABLE]

for any $\delta>0$ and some constant $\omega_{0}$. In particular, $\delta=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 Abanin et al. (2015) (see also Strohmaier et al. (2010); Arad et al. (2016); Abanin et al. (2017), and Appendix F for a self-contained proof) is that, given an $r$-local lattice Hamiltonian $H=\sum_{i}h_{i}$ in any dimension,

[TABLE]

for some $C>0$ and a known $\mathrm{O}(1)$ geometrical factor $G_{r}$. We may conclude $\delta\leq 1$ in (13), so the Lanczos coefficients grow at most linearly.

When linear growth of the $b_{n}$’s is achieved, the growth rate $\alpha$ 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):

[TABLE]

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 $\alpha$ 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 CaF2 McArthur et al. (1969); Engelsberg and Lowe (1974); Lundin et al. (1990). This experimental technique therefore provides a practical way of measuring $\alpha$. On a theoretical note, the spectral function also appears in the off-diagonal Eigenstate Thermalization Hypothesis, which is therefore related to our hypothesis.

Additionally, comparing (14) and (15) shows that $\alpha\leq\pi/2\kappa$, so the growth rate is limited by the local bandwidth of the model and the geometry:

[TABLE]

c.f. (14). This inequality is the consequence of the natural energy scale for the Lanczos coefficients being set by the local bandwidth. However, we shall see that $\alpha$ 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 related to the high-frequency tail of the spectral function $\Phi(\omega)$. On dimensional grounds it is tempting — though ultimately erroneous — to interpret (15) as a statement about the short-time behavior of $C(t)$. To see why this is wrong, notice that the short-time behavior is captured by the first moment alone, as $C(t)=1-\mu_{2}\,t^{2}/2+O(t^{4})$. The high-frequency information instead governs the asymptotics of moments $\mu_{2n}$ as $n\to\infty$ (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\to+\infty$; we will come back to this point in Section VII.2 below. In brief, the hypothesis governs large $\omega$ behavior of $\Phi(\omega)$ and, correspondingly, the behavior of $C(t)$ on the imaginary axis. Explicitly, a growth rate of $\alpha$ gives rise to a singularity at

[TABLE]

IV.2 Analytical Evidence

The upper bounds of the previous section show that the Lanczos coefficients cannot grow faster than linearly. We now show that this bound is tight through two analytic examples.

It is an ironic point that the assumptions for the hypothesis (12) fail in virtually all known solvable models, as those are often integrable, or even non-interacting. This explains why, to the best of our knowledge, linear growth was not recognized in any of the extensive literature on the recursion method as a universal behavior (except for certain classical systems Liu and Müller (1990)). However, there is one solvable model where we can compute the linear behavior analytically: the SYK model (see, e.g. Kitaev (2015); Sachdev and Ye (1993); Maldacena and Stanford (2016)). Its Hamiltonian is

[TABLE]

where the $\gamma_{i}$’s, with $1\leq i\leq N$, are Majorana fermions with anti-commutators $\{\gamma_{i},\gamma_{j}\}=\delta_{ij}$, and the $J_{i_{1}\dots i_{q}}$’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 $\mathcal{O}=\sqrt{2}\gamma_{1}$ Roberts et al. (2018). To leverage the SYK solvability, we shall compute the moments $\mu_{2n}=\mathinner{\left(\mathcal{O}|\mathcal{L}^{2n}|\mathcal{O}\right)}$, 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 Viswanath and Müller (2008). This is described in detail in Appendix B.

We find that the Lanczos coefficients follow the universal form (12) quite closely, as shown in Fig. 4(a). In the large-$q$ limit, there is a closed form expression for the coefficients, computed in Appendix B:

[TABLE]

where $\mathcal{J}=\sqrt{q}\,2^{(1-q)/2}J$. Therefore in the large-$q$ limit, the SYK model follows the universal form (12) with $\alpha=\mathcal{J}$. We may conclude that our hypothesis is obeyed in a canonical model of quantum chaos and that the upper bound of linear growth of the Lanczos coefficients is tight.

The SYK model is quite unusual in several respects: it is a disordered, large-$N$ model in zero dimensions. However, none of these special features are required to achieve linear growth. To demonstrate this we turn to a model studied in the mathematical literature, defined on the 2d square lattice Bouch (2015):

[TABLE]

where $X$ and $Z$ are the normal Pauli matrices. A theorem Bouch (2015) states that the moments of the operator $X_{0,0}$ grow as

[TABLE]

which implies that the Lanczos coefficients grow linearly (see Appendix A for translation between asymptotics). Thus linear growth (12) is a tight-upper bound for the growth of the Lanczos coefficients in dimensions greater than one for “realistic” spin models. The content of our hypothesis is that achieving this upper bound is generic in chaotic systems.

IV.3 The Special Case $d=1$

We now turn to the special case of one dimensional systems. Let us first present some numerical evidence. Fig. 4(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 still appears linear whenever the model is non-integrable. 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 333This is quite fortuitous, computationally: as a general rule, in more strongly interacting systems, exponentially more parameters are required to compute a given $b_{n}$, so fewer $b_{n}$’s may be computed overall.. Fig. 4(b) shows the robustness of this asymptotic 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 appear to become asymptotically linear, and the asymptotic behavior 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 Bohigas et al. (1984); Ullmo (2016). Observe also that the slope of the asymptotic growth depends only weakly on the (integrability breaking) interaction strength. This seems to be a general phenomenon, as it occurs also in the SYK model plus two body interactions, see Fig. 9 for details.

The numerical evidence is apparently compatible with linear growth of the Lanczos coefficients in 1d — but only apparently. We can see this by considering the singularity structure of the correlation function. When the Lanczos coefficients achieve linear growth, there is a singularity in $C(t)$ on the imaginary axis, given by Eq. (3). However, there is a classical theorem Araki (1969) which says, roughly, that $C(t)$, $t\in\mathbb{C}$, is entire for any local system in one dimension. Lanczos coefficients, therefore, must have strictly sublinear growth in one dimension. We note that this is an entirely geometric constraint, and has been previously noted by several works in a variety of contexts Abanin et al. (2015, 2017), and derive it from first principles in Appendix F.

To formulate the hypothesis in one-dimension, we return to the informal version: the Lanczos coefficients should grow as fast as possible. More concretely, the Lanczos coefficients should achieve the upper-bound imposed by the geometry. Following Bouch (2015), we compute this bound in Appendix F and can therefore formulate the hypothesis as follows. Suppose $H$ describes an infinite, non-integrable, many-body system in dimension $d$ and $\mathcal{O}$ is a local operator having zero overlap with any conserved quantity. Then the asymptotic behavior of the Lanczos coefficients is

[TABLE]

for some constants $\alpha,\gamma,A$ and $W$ is the Lambert $W$-function which is defined by the implicit equation $z=W(ze^{z})$ and has the asymptotic $W(n)=\ln n-\ln\ln n+o(1)$. In other words, the hypothesis acquires a logarithmic correction in one dimension. The coefficient $A$, like the growth rate $\alpha$, has dimensions of energy and can be bounded above by the local bandwidth; for Hamiltonians with nearest-neighbor local term $h_{x}$, we have (see Appendix F)

[TABLE]

We note that, unlike in higher dimensions, we are not aware of any analytic examples which achieve the maximal growth rate in 1D, leaving open the possibility that the first line of (22) is an over-estimate.

In some sense, the linear growth “barely breaks” in one dimension; the Lanczos coefficients can still grow faster than $b_{n}\sim n^{\delta}$ for any $\delta<1$. The phenomenological difference between linear growth in all dimensions and (22) is often slight — such as in Fig. 4. Indeed, resolving logarithmic corrections in numerical data is a hard problem that often requires several decades of scaling. Altogether, we see that there is a subtle logarithmic correction to the operator growth hypothesis in one dimension.

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 which have linear growth of Lanczos coefficients with rate $\alpha$. This is done in two steps. First, by studying the quantum mechanics problem (8) on the semi-infinite chain, we show that $\alpha$ measures the rate of exponential growth in operator complexity, in a sense we shall make precise below. Second, we prove that $\alpha$ gives an upper bound on a large class of operator complexity measures. Finally we shall remark on the case of linear growth with log-corrections.

We remark that our notion of complexity is prima facie distinct from other notions bearing the same name, such as circuit complexity (see the reviews Aaronson (2016); Susskind (2018a) and references therein). Indeed, a satisfactory definition of operator complexity of any sort is an unresolved problem, and may not have a unique answer.

V.1 Exponential Growth in the Semi-infinite Chain

Recall that the Lanczos algorithm reduces the operator dynamics to a discrete Schrödinger equation (8),

[TABLE]

We shall analyze this quantum mechanics problem when the hypothesis is satisfied in $d>1$, i.e. $b_{n}=\alpha n+\gamma+o(1)$.

As a first step, we take the continuum limit, by linearizing around momenta [math] and $\pi$. This yields a Dirac equation $\partial_{t}\varphi=\pm 2\alpha x\partial_{x}\varphi$, whose characteristic curves $x\propto e^{\pm 2\alpha t}$ show the wavefunction spreads exponentially fast to the right in the semi-infinite chain with rate $2\alpha$. We remark that among all power-law Lanzcos coefficient asymptotics $b_{n}\sim n^{\delta}$, the linear growth $\delta=1$ is the only one which results in exponential spreading. When $\delta>1$, the characteristic curves reach $x=\infty$ at finite time 444This seems non-physical and indeed, has only been observed in exotic classical systems Liu and Müller (1990). It is ruled out whenever the dynamics are local by Eq. (14).. When $\delta<1$, the spreading follows a power law $x\sim t^{1/(1-\delta)}$. In the case of $d=1$, with the logarithmic correction, the wavefunction spreads as a stretched exponential — faster than any power law, but still slower than exponential.

To undertake a more careful analysis of the wavefunction on the semi-infinite chain, we employ a family of exact solutions. Suppose

[TABLE]

where $\eta=2\gamma/\alpha+1$. For any system when the hypothesis is satisfied, the $b_{n}$’s will approach the $\widetilde{b}_{n}$’s asymptotically, so the properties of the exact solution using the $\widetilde{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 $\widetilde{b}_{n}$’s is

[TABLE]

where $\left(\eta\right)_{n}=\eta(\eta+1)\cdots(\eta+n-1)$ is the Pochhammer symbol and $\mathinner{|\mathcal{O}_{n})}$ is the $n$th Krylov basis vector. Note that this example is not artificial but arises naturally in the SYK model, studied in Section VI.1 below.

The exact solution (25) benefits from a detailed analysis. Recall that the component of the wavefunction at some fixed site $n$ is $\varphi_{n}(t)=(-i)^{n}\mathinner{\left(\mathcal{O}_{n}|\mathcal{O}_{0}(t)\right)}$. For each $n$, $\varphi_{n}(t)$ is a purely real function which starts at [math] (for $n>1$), increases monotonically until reaching a maximum at $t\sim{}\ln n$, then decreases as $\sim{}e^{-\alpha\eta 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 555This shows the importance of the thermodynamic limit. With any finite-dimensional Hilbert space, the chain would be finite, and the results in this section would be affected.. Physically, the wavefunction is decaying by “escaping” off to $n\to\infty$, which serves as a bath. Note, however, that the hypothesis is not sufficient to show that $\varphi_{n}(t)$ decays exponentially with time for small $n$, a fact whose consequences are studied in VII.2 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 $\eqref{eq:universal_wavefunction_main_text}$ has the form $\left|\varphi_{n}(t)\right|^{2}\sim{}e^{-n/\xi(t)}$, where $\xi(t)$ is a “delocalization length” that grows exponentially in time: $\xi(t)\sim e^{2\alpha t}$ for $\alpha t\gg 1$. This exponential spreading is reflected in the expected position of the operator wavefunction (25) on the semi-infinite chain

[TABLE]

More generally, $\mathinner{\left(n^{k}\right)}_{t}\sim e^{2k\alpha t}$ for $k\geq 1$. This result agrees, of course, with the one obtained in the simple continuum-limit above. We believe that the asymptotic growth in (26) 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}=\alpha n+\gamma_{n}$ with various kinds of bounded “impurity” terms $\gamma_{n}\sim\mathrm{O}(1)$. We will consider (26) 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” Susskind (2018b). Indeed, we identify the position on the semi-infinite chain $\mathinner{\left(n\right)}_{t}$ as a notion of operator complexity. We refer to $\mathinner{\left(n\right)}_{t}$ as the “Krylov-complexity” (or “K-complexity” for short) of an operator. After all, as $n$ increases, the operators $\mathcal{O}_{n}$ becomes more “complex”, in the following sense: in the Heisenberg-picture, the equations of motions for $\mathcal{O}_{n}$’s form a hierarchy:

[TABLE]

that is, the dynamics of $\mathcal{O}_{n}(t)$ depends on $\mathcal{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 $\mathcal{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 Prosen (2011).

We know from Section IV.1 that linearly growing Lanczos coefficients are the maximal rate so, in turn, the wavefunction may not spread faster than exponentially. Thus the hypothesis in $d>1$ implies that non-integrable systems have maximal growth of K-complexity: exponential, with rate $2\alpha$.

V.2 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 $\mathcal{Q}$ is a superoperator that satisfies two properties:

[TABLE]

$\mathcal{Q}$ is positive semidefinite. We denote its eigenbasis as $\mathinner{|q_{a})}$, indexed by $a$, so that

[TABLE] 2. 2.

There is a constant $M>0$ such that

[TABLE]

Then q-complexity is defined to be the expectation value

[TABLE]

where $\mathcal{O}(t)$ is evolved under the Liouvillian dynamics of $\mathcal{L}$. A q-complexity is, in principle, an observable, and requires Hamiltonian (or Liouvillian) dynamics. The rationale for the conditions is as follows:  (28a) ensures the q-complexity is always non-negative, (28b) guarantees it cannot change too much under one application of the Liouvillian, and (28c) 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 1: K-complexity.The K-complexity is always a q-complexity, with $\mathcal{Q}=\sum_{n}n\mathinner{|\mathcal{O}_{n})}\mathinner{(\mathcal{O}_{n}|}$. The basis $\mathinner{|q_{a})}$ is just the Krylov basis $\mathinner{|\mathcal{O}_{n})}$ and the conditions (28b) and (28c) are satisfied by construction of the Krylov basis with $M=1$.

Example 2: operator size. A second example of a q-complexity is provided by operator size Roberts et al. (2018). 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 $IXYZII\cdots$ with finitely many non-identity operators. Define $\mathcal{Q}$ to be diagonal in this basis, where the action of $\mathcal{Q}$ on a Pauli string is the number of non-identity Pauli’s. So, for instance, $\mathcal{Q}\mathinner{|IXYZI\cdots)}=3\mathinner{|IXYZI\cdots)}$. The eigenvectors of $\mathcal{Q}$ have non-negative eigenvalues, so $\mathcal{Q}$ is positive semi-definite.

Any choice of dynamics with at most $M$-body interactions (even long-ranged ones) will satisfy (28b), while (28c) requires simple that $\mathcal{O}$ is $d$-local. So, under these conditions, the q-complexity $\mathinner{\left(\mathcal{Q}\right)}_{t}$ becomes the average size of Pauli strings contained in $\mathcal{O}(t)$:

[TABLE]

Example 3: OTOCs. Our third — and most interesting — example of q-complexity is out-of-time-order commutators (OTOCs). Given $\mathcal{O}(t)$, each choice of local operator $V$ defines an OTOC $\mathinner{\left([V,\mathcal{O}(t)]\;|\;[V,\mathcal{O}(t)]\right)}$. 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

[TABLE]

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, (31) is a q-complexity.

To see this, let us work in the eigenbasis of $\mathcal{Q}$. For each site $i$, there is a basis $\mathcal{Q}_{i}\mathinner{|q_{i,a})}=q_{i,a}\mathinner{|q_{i,a})}$ with $1\leq a\leq D^{2}$. We take $\mathinner{|q_{i,0})}$ to be the identity operator with eigenvalue [math], and note that $0\leq q_{i,a}\leq Q$ for some finite $Q$. Since $[\mathcal{Q}_{i},\mathcal{Q}_{j}]=\delta_{ij}$, the eigenbasis for the full operator space is the tensor product of the on-site bases. So for any sequence $\boldsymbol{a}=\left\{a_{i}\right\}$, $\mathinner{|q_{\boldsymbol{a}})}:=\otimes_{i}\mathinner{|q_{i,a_{i}})}$ is an eigenvector satisfying

[TABLE]

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 $\mathcal{Q}$ is positive semidefinite. Since the Hamiltonian is $r$-local, the matrix element $\mathinner{\left(q_{\boldsymbol{a}}|\mathcal{L}|q_{\boldsymbol{b}}\right)}\neq 0$ only if $\boldsymbol{a}$ and $\boldsymbol{b}$ differ on at most $r$ sites. So by (32), we may bound the difference $\left|q_{\boldsymbol{a}}-q_{\boldsymbol{b}}\right|\leq M=rQ$. Similarly, any $r$-local operator satisfies (28c). Having verified all the properties (28), we may conclude that OTOCs of this form are a q-complexity.

OTOCs are known to be closely related to the operator size Roberts et al. (2018); Maldacena et al. (2016). 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. We remark that q-complexities (including K-complexity) are quadratic in $O(t)$ and not linear response quantities, although the growth rate $\alpha$ is, via the spectral function. We will see in Section VI.2 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,

[TABLE]

The following section will focus on the application of this general bound in the specific case of OTOCs.

To close this section, we show how the above results are affected by the log-correction to linear growth in 1d from Eq. (22): $b_{n}\sim An/\ln n$. The continuum Dirac equation analysis yields a stretched exponential growth of K-complexity:

[TABLE]

which is slower than any exponential growth but faster than any power law. Combined with (33), we conclude that all q-complexities have at most stretched exponential growth in 1d.

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 (33) and (26), we see that q-complexities grow at most exponentially in time, at least when the hypothesis holds for $d>1$. If that is the case, with $\mathinner{\left(Q\right)}_{t}\sim{}e^{\lambda_{Q}t}$, then the exponent is bounded above by $2\alpha$:

[TABLE]

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,

[TABLE]

its growth rate $\lambda_{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 666To be precise, the OTOC measures a generalized Lyapunov exponent with $q=2$, which is greater or equal to the typical one Politi (2013). We can then state following bound on Lyapunov exponents: for any system at infinite temperature where the operator growth hypothesis holds, then

[TABLE]

where we put $\lambda_{L}=0$ whenever the OTOC grows slower than exponentially, and similarly for $\alpha$. This follows directly from (33) and (26), so we have essentially proven (37) as a mathematical proposition.

It is interesting to note that, as $\lambda_{L}$ is defined via a four-point correlation function (the OTOC), while $\alpha$ depends on a two-point correlation function ($C(t)$), the bound (37) can be interpreted as a relation between correlation functions of distinct nature. Such a relation is, to our knowledge, rather unusual (see Murthy and Srednicki (2019) for a recent result). However, this point of view is not how we derived (37); an alternative proof working directly with the correlation functions would be illuminating.

Remarkably, the bound (37) appears to be valid under much less restrictive assumptions — at any temperature and in either classical or quantum systems. In this section, we examine the cases of quantum and classical systems at infinite temperature, and leave that of finite temperatures to Section VIII below.

VI.1 SYK Model

We illustrate the bound (37) for the SYK model (18). At infinite temperature, no analytic formula for the Lyapunov exponent is available, but it has been computed numerically in, e.g. Roberts et al. (2018); Maldacena and Stanford (2016). Table 1 shows that not only does (37) hold for the whole range of $q$-SYK models, but $\alpha$ is almost equal to $\lambda_{L}/2$, with exact agreement in the limit $q\to\infty$ 777Indeed, the difference may well be a numerical effect, see Roberts et al. (2018).. These results show that the bound $\lambda_{L}\leq 2\alpha$ is tight: the prefactor cannot be improved in general. Moreover, in the large $q$ limit, the probability distribution $\left|\varphi_{n}(t)\right|^{2}$ on the semi-infinite line is identical to the operator size distribution of $\gamma_{1}(t)$ Roberts et al. (2018). (See (95) 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 (37) can be nearly saturated. For spatially extended quantum systems with finitely many local degrees of freedom, Lieb-Robinson bounds Lieb and Robinson (1972) and its long-range generalizations Else et al. (2018) guarantee that the OTOC has slower-than-exponential growth in most physical systems at infinite temperature 888Indeed, generalized Lieb-Robinson bounds state that the OTOC between $\mathcal{O}(t)$ and $V_{i}$ is exponentially small if the site $i$ lies out of some volume which grows sub-exponentially. Then, a sum like (31) is essentially that volume..

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 few-body 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.

VI.2 Classical Chaos

We now transition to the classical setting. After briefly explaining how the recursion method carries over almost verbatim to classical systems, we shall examine the classical form of the bound (37). However, the arguments of Section V.2 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).

VI.2.1 A (Weaker) Bound on Classical Chaos

The recursion method applies to classical and quantum systems in exactly the same manner Viswanath and Müller (2008). Classically, operator space is the space of functions on classical phase space and the Liouvillian $\mathscr{L}=i\{\mathscr{H},\cdot\}$ is defined by the Poisson bracket against the classical Hamiltonian $\mathscr{H}$ (we take $\hbar=1$). The appropriate classical inner product at infinite temperature is $(f|g)=\int f^{*}g\,d\Omega$, where $d\Omega$ is the symplectic volume form on the phase space 999We therefore require a compact phase space, such as in a classical spin model.. The Liouvillian $\mathscr{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 Liu and Müller (1990); Viswanath and Müller (2008) that linear growth of the Lanczos coefficients appears in general finite-dimensional, non-linear systems, to which we restrict ourselves 101010Note that even if the phase space is finite-dimensional, the operator space is infinite-dimensional, allowing an infinite sequence of Lanzcos coefficients.. The growth rate $\alpha$ is well-defined in such systems, as is the (classical) Lyapunov exponent $\lambda_{L}$, and the bound (37) takes on the same form as before: $\lambda_{L}\leq 2\alpha$. 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 $\mathinner{|Y^{m}_{\ell})}$, $\ell=0,1,2\dots$, $m=-\ell,\dots,\ell$.

A typical Hamiltonian is a polynomial of the classical spin operators $\mathscr{S}^{x},\mathscr{S}^{y},\mathscr{S}^{z}$ with Poisson bracket $\{\mathscr{S}^{a},\mathscr{S}^{b}\}=-\varepsilon^{abc}\mathscr{S}^{c}$. We consider the simple non-linear example

[TABLE]

Using Clebsch-Gordon coefficients one can show that the classical Liouvillian is quite sparse, and only the following matrix elements are non-zero:

[TABLE]

whenever the states in question exist.

We now examine the classical OTOC for the local operator $\mathscr{S}^{z}$, given by matrix elements of a super-operator $\mathscr{Q}^{z}$. This operator is diagonal in the basis of spherical harmonics

[TABLE]

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\pm 2m$, which can be arbitrarily large. Hence the condition (28b) cannot be satisfied for any finite constant $d$. It is helpful to recall that Section V.2 showed the quantum OTOC is a q-complexity whenever the on-site Hilbert space is finite-dimensional. This fails in the case of a spin $s$, whose on-site dimension $2s+1$, in the classical limit $s\to\infty$. We have therefore seen that classical OTOCs are not q-complexities and, hence, the bound (37) does not follow from the reasoning of Section V.2 in the classical case, and remains a conjecture.

Nonetheless, for any Hamiltonian and initial operators that are polynomials of the spin variables $\mathscr{S}^{a}$, we can show the following general bound

[TABLE]

which is weaker than the conjectured $\lambda_{L}\leq 2\alpha$.

To show (41), observe that by (40), the superoperator $\mathcal{R}_{z}:=\mathcal{Q}_{z}^{\frac{1}{2}}$ satisfies (28b), since its has eigenvalue $m$ for $Y^{\ell}_{m}$, which can change only by $\delta$ upon one Liouvillian application, where $\delta$ is the polynomial degree of the Hamiltonian. Other conditions in (28) are satisfied straightforwardly. We then have

[TABLE]

which implies (41). Here the first $\sim$ is by definition, the the inequality is a straightforward generalization of the bound on q-complexity, Eq. (127) of Appendix E, and the last $\sim$ is a generalization of (26) (see below that equation).

This argument carries over to the OTOC with spin variables in any direction by spherical symmetry, and applies almost verbatim to systems with a few spins, $\mathscr{S}_{i}^{x,y,z},i=1,\dots,N$. A Lyapunov exponent associated with a finite sum such as

[TABLE]

satisfies the same bound since every term does so. In summary, (41) is established in general classical few-spin models. We expect it is possible to show (41) rigorously.

An interesting corollary of (41) is a relation between chaos and the decay rate of the spectral function. Recall that the linear growth of Lanczos coefficients is equivalent to the exponential decay of the spectral function $\Phi(\omega)\sim\exp(-|\omega|/\omega_{0})$ at high frequency, where $\omega_{0}=\tfrac{2}{\pi}\alpha$. Then (41) is equivalent to

[TABLE]

(The conjectured bound would instead imply $\lambda_{L}\leq\pi\omega_{0}$.) In numerous classical systems, the power spectrum decay of time series has been used as an empirical probe of deterministic chaos Frisch and Morf (1981); Greenside et al. (1982); Sigeti (1995a, b); Cheskidov et al. (2008); Elsayed et al. (2014); Maggs and Morales (2011). To the best of our knowledge, the bound (44) provides the first quantitative justification for this usage.

We mention that the relation between chaos and long-time 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 Ruelle (1986); Turiaci and Verlinde (2016). 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, but not a diagnostic of classical chaos. Indeed, our bound is correct but not tight for most classical integrable systems which, generically, have non-zero growth rate but no chaos  Liu and Müller (1990).

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 (37) in a canonical example of classical chaos.

VI.2.2 Numerical Case Study

The Feingold-Peres model of coupled tops Feingold and Peres (1983) is a well-studied model of few-body chaos, both classically and at the quantum level Feingold et al. (1984); Fan et al. (2017). The quantum model is a system of two spin-$s$ particles, $1$ and $2$, with Hamiltonian

[TABLE]

where $c\in[-1,1]$ is a parameter and $S_{i}^{\alpha}$ satisfy the $SU(2)$ algebra $[S_{i}^{\alpha},S_{j}^{\beta}]=i\hbar\delta_{ij}\varepsilon^{\alpha\beta\gamma}S_{i}^{\gamma}$ and act on a spin-$s$ Hilbert space. This is non-interacting when $c=\pm 1$ and chaotic in the intermediate region. Correspondingly, the Lanczos coefficients are asymptotic to a constant at $c=\pm 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

[TABLE]

where $\mathscr{S}_{i}^{\alpha},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 $\mathcal{O}\propto S_{1}^{z}S_{2}^{z}$ ($\mathscr{S}_{1}^{z}\mathscr{S}_{2}^{z}$ in the classical case). As shown in Fig. 6(b), the quantum Lanczos coefficients converge to the classical ones as $s\to\infty$, as expected, and they increase linearly near $c=0$. We have checked that $\alpha$ does not depend on the choice of initial operator $\mathcal{O}$, so long as $\mathcal{O}$ does not overlap with any conserved quantity.

To test the conjectured bound (37), we compare the growth rate $\alpha$ with the classical Lyapunov exponent ($\lambda_{L}/2$ in our notation), which can be calculated by the standard variational equation method Alligood et al. (1996). Remarkably, the data shown in Fig. 6(a) corroborates the conjectured bound $\alpha\geq\lambda_{L}/2$ in the parameter region explored, with equality up to numerical accuracy in the regime $c\approx 0$, where the model is known to be maximally chaotic, with almost no regular orbits Feingold et al. (1984); Feingold and Peres (1983). Enlarging the parameter space, for instance by adding terms such as $\mathscr{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 White (1992); Osborne and Nielsen (2002). 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.

VII.1 Continued Fraction Expansion: Brief Review

We briefly review the continued fraction expansion of the Green’s function Viswanath and Müller (2008). The Green’s function (9) is related to the autocorrelation $C(t)$ by the following transform:

[TABLE]

where the integration contour is taken to be the shifted real axis shifted down by $-i\epsilon$ for some small $\epsilon>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 (47) as the Laplace transform, despite the fact that it differs from the usual definition by a factor of $i$.

In the Krylov basis, $G(z)=\left[z-L\right]^{-1}_{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\geq 2$, and then return. More formally, for each $n\geq 0$, let $L^{(n)}:=L_{p\geq n,q\geq n}$ be the hopping matrix on the semi-infinite chain restricted to sites $n$ and above, and let $G^{(n)}(z):=\left[z-L^{(n)}\right]^{-1}_{nn}$ be the corresponding Green function. (Note that $G^{(0)}(z)=G(z)$.) We then have the following recursion relation — hence the name “recursion method” —

[TABLE]

(For a quick derivation Auerbach (2018), 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_{n+1}^{2}P_{n+2}(z)$. Then (48) follows from the two preceding equations.)

Applying Eq. (48) recursively yields the continued fraction expansion:

[TABLE]

To save space, we denote the recursion 48 by $G^{(n)}=M_{n+1}\circ G^{(n+1)}$, where $M_{n}$ is the Möbius transform $w\mapsto 1/(z-b_{n}^{2}w)$ and “$\circ$” 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 Viswanath and Müller (2008).

VII.2 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 $\mathcal{O}_{q}=\sum_{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\gamma$, then the autocorrelation function will decay exponentially as $e^{-\gamma 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)=\left[z-iD(z)q^{2}/2\right]^{-1}$, where $D(z)=D_{0}+D_{1}\sqrt{z}$. 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 Mukerjee et al. (2006) as a power law in time 111111We thank Achim Rosch for pointing out this possibility.. 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 Medenjak et al. (2017); Ljubotina et al. (2017); Nardis et al. (2019); Gopalakrishnan and Vasseur (2019); De Nardis et al. (2018a, b) for recent developments). Indeed, accurately computing diffusion coefficients has been the goal of much recent numerical work Leviatan et al. (2017); White et al. (2018); Hauschild et al. (2018). This difficulty is reflected in the continued fraction expansion (49): 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 $\eqref{eq:linear_lanczos}$ 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)=\operatorname{sech}(\alpha t)$ and $C_{2}(t)=\left(1+t^{2}\right)^{-\gamma}$ Viswanath and Müller (2008) both correspond to Lanczos coefficients that grow linearly But $C_{1}(t)$ decays exponentially while $C_{2}(t)$ decays as a power law, so clearly the asymptotics of $b_{n}$ 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}=\alpha n+\gamma+(-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. In particular, the results in this work are prima facie unrelated to a bound on transport Hartman et al. (2017).

VII.3 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 (12) especially quickly and regularly, 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 Viswanath and Müller (2008); Auerbach (2018). 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 semi-infinite chain. If the hypothesis is obeyed, then the $b_{n}$’s also asymptotically approach the form

[TABLE]

where $\eta=2\gamma/\alpha+1$. The agreement is better, of course, at large $n$. The coefficients $\widetilde{b}_{n}$ have the virtue that the quantum mechanics problem they describe on the semi-infinite chain is exactly solvable. Appendix D applies the theory of Meixner orthogonal polynomials of the second kind to determine the autocorrelation analytically: $C(t)=\operatorname{sech}(\alpha t)^{\eta}$. (This is the same exact solution used in Section V above.) By Laplace transform, the corresponding Green’s function is

[TABLE]

Here ${}_{1}F_{2}$ is the hypergeometric function and $\widetilde{M_{k}}$ depends on $\widetilde{b}_{k}$. It is crucial that $\widetilde{G}^{(n)}(z)$ is known analytically, so that (51) provides the asymptotically exact large $n$-behavior.

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\geq N$. So for sufficiently large $N$, where the $b_{n}$’s are approximately the same as the $\widetilde{b}_{n}$’s, we may approximate

[TABLE]

an approximation that becomes better at large $N$. Equation (52) 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 $\alpha$ and $\eta$. Then the approximate $G(z)$ can be calculated from (51) 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$. This procedure is illustrated for the energy diffusion in chaotic Ising model in Fig. 7. Almost all the computational effort goes into in computing the first few $b_{n}$’s exactly. We also note that the extrapolation is carried out with a linear fit to the Lanczos coefficients which is not strictly appropriate to $d=1$ (the log-correction is missing). Nevertheless, the numerical value of the diffusion coefficient appears to match other methods to within a few percent.121212We are greatful to Francisco Machado and Biantian Ye for sharing their density matrix truncation (DMT) results with us. Further numerical tests on this example indicate the the exact asymptotics of Lanczos coefficients may not be necessary to compute $D$ to a decent precision.

In short, the hypothesis is sometimes sufficient to describe the emergent hydrodynamic behavior of operators, even if we ignore the log correction in 1d. 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. On the other hand, in the better scenarios, less knowledge on the Lanczos coefficients is required to capture the hydrodynamic coefficients. We will provide further examples of this algorithm and discuss its theoretical and practical accuracy in subsequent work.

VIII Finite Temperature

So far our discussion has been confined to infinite temperature. In this section we generalize to finite temperature. Only a minor modification is required to carry out the Lanczos algorithm at finite temperature so many of our results carry over unaffected. A summary is provided in Table 2 for the reader’s convenience.

VIII.1 Choice of Inner Product

A single modification is required to adapt the formalism of recursion method to finite temperature: an operator inner product which incorporates the thermal density matrix. At temperature $T=1/\beta$ (we set $k_{\text{B}}=1$), a general operator scalar product is defined by the integral Viswanath and Müller (2008):

[TABLE]

where $g(\lambda)$ is some even function on the thermal circle $[0,\beta]$, $y:=e^{-H}$, and $Z:=\operatorname{Tr}[y^{\beta}]$ 131313Precisely, $g$ must satisfy $g(\lambda)\geq 0$, $g(\beta-\lambda)=g(\lambda)$, and $\beta^{-1}\int_{0}^{\beta}d\lambda g(\lambda)=1$. We also restrict to the subspace of operators with zero thermal expectation value, and omit the disconnected term in (53).. The choice of the inner product is not arbitrary, but is equivalent to the choice of the correlation function

[TABLE]

(where $\rho_{\beta}=e^{-\beta H}/Z$), which is in turn determined by the physical context; in fact, only a few choices of $g$ are physically relevant, such as (59) and (60) below.

Once the inner product is chosen, the Lanczos coefficients are defined by the same Lanczos algorithm with the new norm. Quite explicitly, the recursion is:

[TABLE]

for $n=1,2,3,\dots,$ starting from $\mathinner{|\mathcal{O}_{0})}_{\beta}^{g}:=\mathinner{|\mathcal{O})}$, $\mathinner{|\mathcal{O}_{-1})}_{\beta}^{g}:=0$ and $b_{0,T}^{(g)}:=0$. We emphasize that only the inner product has been changed compared to the infinite-$T$ version. In fact, the Krylov subspaces $\mathrm{span}\{\mathinner{|\mathcal{O})},\mathcal{L}\mathinner{|\mathcal{O})},\dots,\mathcal{L}^{n}\mathinner{|\mathcal{O})}\}$ are unchanged at finite temperature, and only the notion of orthogonality is different, giving us a new orthogonal basis for those spaces. Also, we have the same relationships between the Lanczos coefficients and the correlation function (54), as well as its linear transforms, the Green’s function and spectral function

[TABLE]

where the superscript $g$ is not an exponent. For example, the Green function (56a) admits the continuous fraction expansion

[TABLE]

which is identical to (49), except that $b_{n}$ are replaced by the finite-$T$ Lanczos coefficients. Similarly, the results of Appendix A carry over directly.

The statement of the hypothesis at finite temperature is also directly analogous. We hypothesize that a chaotic system should have maximal growth of the Lanczos coefficients,

[TABLE]

under the same conditions as before. Here $\alpha_{T}^{(g)}\geq 0$ depends on the inner product. Evidence for the hypothesis at finite $T$ will be provided in Section VIII.3.

Though the Lanczos algorithm proceeds in the same way for any choice of inner product, this choice will determine what physical correlation function we end up computing. There are two prominent choices of inner products:

• •

In linear response theory, we use the “standard” inner product given by $g(\lambda)=[\delta(\lambda)+\delta(\lambda-\beta)]/2$:

[TABLE]

that leads to the usual thermal correlation function.

• •

In quantum field theory, it is often natural to consider the Wightman inner product, which corresponds to $g(\lambda)=\delta(\lambda-\beta/2)$:

[TABLE]

In particular, this inner product allows us to relate our bound on chaos (37) and the finite-temperature bound of Ref. Maldacena et al. (2016).

In equations (59) and (60) and below, we replace the $g$ by $S$ or $W$ to indicate the choice of standard and Wightman inner product, respectively. At infinite temperature, both inner products reduce to the one $\mathinner{\left(A|B\right)}=\operatorname{Tr}[A^{\dagger}B]/\operatorname{Tr}[1]$ considered previously.

The spectral functions of the two choices are related by a well-known identity:

[TABLE]

which follows directly from the definition (10). The Wightman inner product therefore imposes an extra temperature-dependent exponential decay to the spectral function, due to the suppression of high energy excitation by the two $e^{-\beta H/2}$ factors in (60). This observation will be crucial in the following section. On the other hand, it would be very interesting to understand how the high-frequency tail of $\Phi(\omega)_{\beta}^{S}$ depends on the temperature.

VIII.2 Bound on Chaos

A key result on quantum chaos at finite temperature is the bound on chaos of Ref.  Maldacena et al. (2016). This universal bound was derived for quantum field theories at finite temperature $T=\beta^{-1}$, and reads as follows

[TABLE]

in natural units $\hbar=k_{\mathrm{B}}=1$. It is nontrivial in finite-temperature quantum systems, and is therefore complemented by our bound $\lambda_{L}\leq 2\alpha$ (37) which applies to infinite temperature quantum and classical system. This leads to two natural questions. Can our bound be extended to finite temperature? How does it compare to the universal one?

Since $\alpha_{T}^{(g)}$ depends on the inner product, and the finite-$T$ OTOC admits various regularizations, it is already a nontrivial task to find the correct formulation of the extension. To make progress we consider the regularization scheme used for four-point OTOCs in Maldacena et al. (2016) to derive the universal bound. This scheme inserts the operators in the thermal circle $[0,\beta)$ with even spacing, as does the Wightman inner product (60). This suggests that an extension of the bound $\lambda_{L}\leq 2\alpha$ to finite temperature can be obtained by comparing the finite-T Lyapunov exponent (as defined in Maldacena et al. (2016)) and the finite-T growth rate defined with the Wightman inner product:

[TABLE]

We stress that this is a conjecture below infinite temperature. Nevertheless, as we show in Section VIII.3 below, exact results in the $q$-SYK model suggest that (63) is plausible and tight.

We now turn to the relation between the conjecture (63) and the universal bound, and show that the former infers the latter. By (61), the Wightman spectral function decays at least as fast as $e^{-\beta\omega/2}$ at high frequency (because $\Phi_{\beta}^{S}(\omega)\leq 1$). By (15), this is equivalent to the following upper bound on the Lanczos coefficients growth rate:

[TABLE]

where $\alpha_{T}^{(W)}$ denotes the growth rate with Wightman inner product. Therefore, the conjecture (63), if true, would be tighter than the universal one $\lambda_{L,T}\leq 2\pi T$ (62). At low temperature ($\beta\to\infty$ limit), the decay of $\Phi_{\beta}^{W}(\omega)$ is dominated by the factor $e^{-\beta\omega/2}$, so $\alpha_{T}^{(W)}/(\pi T)\to 1$ and the conjectural bound (63) becomes equivalent to the universal one (62). This equivalence suggests further the plausibility of the conjecture (63).

VIII.3 SYK Model

To illustrate the foregoing discussion, and provide some evidence for the hypothesis at finite-$T$ (58) and the conjectural bound on chaos (63), let us consider again the example of SYK model.

At low temperatures $T=1/\beta\ll J$, it is well-known that $\lambda_{L,T}=2\pi T$ Kitaev (2015) saturates the universal quantum bound (62). In this limit, the finite-$T$ autocorrelation function of $\mathcal{O}=\sqrt{2}\gamma_{1}$ may be computed exactly by conformal invariance Maldacena and Stanford (2016). Choosing the Wightman inner product, we have

[TABLE]

This is the autocorrelation function of the exact solution (25), and corresponds to Lanczos coefficients $b_{n,T}^{(W)}=\pi T\sqrt{n(n-1+\eta)}$. They satisfy the hypothesis (58) with $\alpha_{T}^{(W)}=\pi T$ (64). Therefore the low-temperature SYK model saturates also our conjectural bound (63).

At finite (but not necessarily low) temperatures, using analytic results in the large-$q$ limit Maldacena and Stanford (2016), it is not hard to check (see Appendix B) that our conjectured bound (63) is saturated, whereas the universal bound (62) is not, see Fig. 8. This result indicates that an extension of our bound on chaos to finite temperature is at least plausible. The exact agreement between $\alpha_{T}^{(W)}$ and $\lambda_{L,T}$ is notable given that the former is defined solely from 2-point correlators whereas the latter requires 4-point functions.

We reiterate that the above SYK results depend crucially on the Wightman inner product. If the “standard” inner product (59) is chosen instead, the Lanczos coefficients $b_{n,T}^{(S)}$ cannot be extracted from the conformal solution, since that would require the Taylor expansion of $C_{\beta}^{S}(t)$ around $t=0$, at which the conformal solution is non-analytic. A numerical high-temperature expansion (extending the method of Appendix B) and an exact calculation in the large-$q$ limit both indicate that the Lanczos coefficients still grow linearly, but the growth rate increases as the temperature decreases.

To summarize, exact calculations in the SYK model support the universal operator growth hypothesis at finite temperature, and the conjectural bound on chaos.

IX Conclusions

IX.1 Discussion

We have presented a hypothesis on the universal growth of operators: the Lanczos coefficients follow the asymptotically linear form $b_{n}=\alpha n+\gamma+\mathrm{o}(1)$ in non-integrable systems, with a logarithmic correction in 1d. We have seen copious evidence that the hypothesis is satisfied in a wide variety of non-integrable models. Over the course of this work, the growth rate $\alpha$ has emerged as a quantity of prime importance, tying a diverse array of seemingly-disparate ideas together. Let us recount them now:

• •

$\alpha>0$ is the slope of asymptotically linear growth of the Lanczos coefficients.

• •

$\tfrac{2}{\pi}\alpha=\omega_{0}$ is the exponential decay rate of the spectral function $\Phi(\omega)\sim e^{-|\omega|/\omega_{0}}$, which can be (and has been) measured experimentally McArthur et al. (1969); Engelsberg and Lowe (1974); Lundin et al. (1990).

• •

$\pm i\pi/(2\alpha)$ are the locations of the singularities closest to the origin in the (analytic continuation) of the autocorrelation $C(t)$, see Appendix A.

• •

$2\alpha$ is the exponential growth rate of Krylov-complexity.

• •

$2\alpha$ is an upper bound for the growth of all q-complexities.

• •

$2\alpha$ 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, $\alpha$ plays a central role in operator growth and dynamics of complex systems.

Complexity — especially the Krylov-complexity — arose 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 $\hbar$ 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 $\hbar$ 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.

IX.2 Outlook

We would like to understand how our hypothesis can be affected by obstructions to thermalization. Based on evidence available to us, it is tempting to conjecture that they lead to a qualitative slower growth for quantum systems. Confirming this in general would be a remarkable result. However, given the diversity of non-thermalizing situations, it may be more reasonable to explore them on a case by case basis. In free and integrable models, there are an extensive number of conserved local or quasi-local charges. The behavior of the Lanczos coefficients in integrable models is likely non-universal, and depends strongly on the model and operator in question Viswanath and Müller (2008). We wish to gain general analytical insights in this direction (especially for interacting models), by leveraging the knowledge available on the quantum inverse scattering method Faddeev (1999); Kitanine et al. (1999); Maillet and Terras (2000). Also, it may be desirable 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 Bernien et al. (2017); Turner et al. (2018); Choi et al. (2018), 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 Abanin et al. (2018) and references therein for a review, and Khait et al. (2016) for numerical calculations of Lanczos coefficients in disordered spin chains) where thermalization fails.

Our treatment at finite temperatures is far from complete and leaves numerous open questions, especially those concerning the “standard” inner product: How do the Lanczos coefficients grow? If linearly, how does the growth rate depend on the temperature? How can we extend our bound on chaos to finite $T$? Numerical investigations into these questions are challenging due to the presence of the thermal density matrix Lucas (2018); Auerbach (2018); Martyn and Swingle (2018). 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, especially in 1d. Finding analytically tractable models far from the large-$N$ limit that achieve the maximal Lanczos coefficient growth seems a formidable problem, which is made even harder by the restriction to time-independent Hamiltonian systems; the only result in this direction is that of Bouch (2015) in 2d. Many solvable models of quantum chaos (see Refs Chan et al. (2018); Bertini et al. (2018) for notable recent progress) are only defined as unitary maps or Floquet systems. To this respect, a meaningful extension of the hypothesis to such contexts would be a highly rewarding advance.

An alternative route would be to develop an extended (Hermitian) 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 Tao (2012). These are directly analogous to the moments $\mu_{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 speculate that the hypothesis can be derived analytically by introducing a new type of random matrix ensemble that incorporates locality and translation invariance. (This is similar to the framework of Movassagh and Edelman (2011).) In this case, a Hamiltonian such as $H=\sum_{<x,y>}h_{x,y}$, where $h_{x,y}$ is a random matrix acting on neighboring sites $x$ and $y$, should obey the hypothesis (12) 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.

Coming back to physics, we argue that there should be a general principle, analogous of the second law of thermodynamics, that governs the operator growth in generic systems. Indeed, the latter is irreversible, in the same sense as the dynamics of an isolated gas is so in the thermodynamic limit. We cannot help but wonder what entropy is maximized by the operator growth process, and whether any notion of (quantum) dynamical entropy (see e.g. Connes et al. (1987); Alicki and Fannes (1994); Benatti (2012); Prosen (2007) is relevant in describing the process. Elusive as it seems, such a thermodynamic principle might be the ultimate explanation of our empirical observations of ubiquitous maximal operator growth.

To close, we wish to point out that the territory of q-complexities 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 circuit complexity and entanglement entropy. Hopefully, charting this terra incognita will continue to shed new light on the complex nature of many-body quantum dynamics.

Acknowledgments. We thank Assa Auerbach, Erez Berg, David Huse, Joel Moore, Steve Shenker, and Achim Rosch for insightful discussions. Part of the numerical computations are performed using the cluster of Laboratoire de Physique Théorique et Modèles Statistiques (CNRS, Université Paris-Sud). We acknowledge support from NSF Graduate Research Fellowship Program NSF DGE 1752814 (DP), the Emergent Phenomena in Quantum Systems initiative of the Gordon and Betty Moore Foundation (TS), ERC synergy Grant UQUAM (EA, XC and TS) and DOE grant DE-SC0019380 (EA and XC).

Appendix A Brief Review of the Recursion Method

In this appendix we recall the relations between Lanczos coefficients, correlation function, Green function, spectral function, and moments. These relations are mathematical in nature, and apply to any inner product on the operator space, and thereby to finite as well as infinite temperature. For simplicity, we will omit the sub- and superscripts indicating the inner product.

Let us recall the five equivalent representations of the dynamics of an operator:

[TABLE]

The first four are related by linear transformations given in the text. For instance, the moments $\mu_{2n}$ are the Taylor expansion coefficients of autocorrelation around $t=0$:

[TABLE]

where the odd terms vanish provided $\mathcal{O}$ is Hermitian. The moments can also be extracted from the spectral function via

[TABLE]

All the transformations between the first four quantities are similarly straightforward.

The Lanczos coefficients, on the other hand, are related to the others via a non-linear transformation. The rest of this Appendix discusses how to perform the non-trivial translation between the Lanczos coefficients and the moments both asymptotically and numerically.

A.1 From Moments to Lanczos Coefficients

Cumulative products of the first $n$ Lanczos coefficients are given by determinants of the Hankel matrix of moments Viswanath and Müller (2008)

[TABLE]

If the moments are known, the determinant can be computed efficiently by transforming the Hankel matrix into diagonal form. Doing this iteratively for $k\in[1,n]$ provides a fast algorithm that computes $b_{1},\dots,b_{n}$ from $\mu_{2},\mu_{4},\dots,\mu_{2n}$. The algorithm may be expressed concisely as a recursion relation (see Eq. 3.33 of Ref. Viswanath and Müller (2008)) as follows:

[TABLE]

If an analytic expression for $C(t)$ is known, then an arbitrary number of the Lanczos coefficients may be computed numerically via (70). We remark that this algorithm suffers from large numerical instabilities due to repeated floating-point divisions.

A.2 From Lanczos Coefficients to Moments

It follows from the tridiagonal form of $L$ that the moments may be expressed in terms of the Lanczos coefficients as

[TABLE]

If the Lanczos coefficients are known, this is a completely combinatorial object. In particular, the moments are given by a sum over Dyck paths. Formally, a Dyck path of length $2n$ can be defined as a sequence $(h_{0},h_{1},\dots,h_{2}n)$ such that: $h_{0}=h_{2n}=1/2$; $h_{k}\geq\frac{1}{2}$ and $|h_{k}-h_{k+1}|=1$ for any $k$. These are often visualized as paths starting at height zero where each segment either increases or decreases the height by one unit, with the constraint that the height is always non-negative and returns to zero at the end. Denoting the set of such paths by $\mathcal{D}_{n}$, we have

[TABLE]

For example, $\mu_{2}=b_{1}^{2}$ and $\mu_{4}=b_{1}^{4}+b_{1}^{2}b_{2}^{2}$. The number of Dyck paths of length $2n$ is given by the Catalan numbers $C_{n}=\frac{(2n)!}{(n+1)!n!}$. A consequence of (72) is the following lower bound:

[TABLE]

On the other hand, we have the upper bound $\mu_{2n}\leq\max_{k=1}^{n}\left(b_{k}^{2}\right)C_{n}$. Applying the upper and lower bounds, linear growth of the Lanczos coefficients $b_{n}$ corresponds to the following growth rate of moments:

[TABLE]

This equation is a useful reformulation of the linear growth hypothesis.

If the growth rate is known as well, $b_{n}=\alpha n+\mathrm{O}(1)$, it is possible to refine the asymptotic by specifying the next order exponential term:

[TABLE]

Combining this equation with the Stirling formula, the correlation function $C(t)=\sum_{n}\mu_{2n}(it)^{2n}/(2n)!$ has convergence radius $r=\pi/(2\alpha)$, due to singularities at $t=\pm ir$; in fact, $C(t)$ is analytical in the strip $-r<\operatorname{Im}(t)<r$, see Fig. 3. Therefore, the Fourier transform of $C(t)$, which is the spectral density $\Phi(\omega)$, has a exponential decay

[TABLE]

We illustrate the above results by a simple example: when $b_{n}=\alpha\,n$, then $C(t)=\operatorname{sech}(\alpha t)$ and $\Phi(\omega)=\frac{\alpha}{\pi}\operatorname{sech}\left(\frac{\pi\omega}{2\alpha}\right)$. The moments $\mu_{2n}=1,1,5,61,1385,\dots$ are known as Euler or secant numbers and have the asymptotic behavior $\mu_{2n}=4\sqrt{\frac{4n}{\pi}}\left(\frac{4n}{\pi e}\right)^{2n}(1+\mathrm{o}(1))$ Borwein et al. (1989). We checked that (76) and  (75) hold in all analytic examples we are aware of in the literature and believe them to hold in general.

Appendix B Moments and Lanczos Coefficients in the SYK Model

In this section we compute the Lanczos coefficients in the large-$N$ SYK model at infinite temperature with the initial operator $\mathcal{O}=\sqrt{2}\gamma_{1}$. Most often, this is done by computing the moments and applying the mapping described in Section A.

For convenience, we recall the SYK Hamiltonian and disorder normalization:

[TABLE]

where the line denotes disorder averages. We shall extend $J_{i_{1}\dots i_{q}}$ to all $i_{1},\dots,i_{q}$ by anti-symmetry. As discussed in the main text, disorder-averaging will be assumed throughout. We first describe the general method, and then discuss the large-$q$ limit.

B.1 General Method

Since the moments are closely related to the Green function, they can be calculated by the diagrammatic technique commonly used in the SYK literature. Indeed, $\mu_{2n}$ can be represented as a sum over diagrams $G$ diagrams with $2n$ vertices:

[TABLE]

where $C_{G}$ is the combinatorial factor of the diagram, which counts the number of labellings of the vertices by $1,\dots,2n$ such that the labels are increasing from left to right.

Let us illustrate the diagrams with some examples with $q=4$ and $n=1,2$. Direct calculation yields:

[TABLE]

The first two moments $\mu_{2}$ and $\mu_{4}$ are (twice) the norm squared of the $\mathcal{L}\gamma_{1}$ and $\mathcal{L}^{2}\gamma_{1}$, respectively. Under disorder averaging, the terms on the right-hand side are orthogonal, and each corresponds to a different diagram:

[TABLE]

The combinatorial factor is $C_{G}=1$ for each of the above graphs. The first non-trivial combinatorial factor is$C_{G}=6$ for the diagram , which contributes to $\mu_{6}$. The six vertex orderings are $1\;\begin{matrix}2&3\\ 4&5\end{matrix}\;6$, $1\;\begin{matrix}4&5\\ 2&3\end{matrix}\;6$, $1\;\begin{matrix}2&4\\ 3&5\end{matrix}\;6$, $1\;\begin{matrix}3&4\\ 2&5\end{matrix}\;6$, $1\;\begin{matrix}2&5\\ 3&4\end{matrix}\;6$, and $1\;\begin{matrix}3&5\\ 2&4\end{matrix}\;6$.

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:

[TABLE]

that is, $\widetilde{\Sigma}(z)$ and $\Sigma(t)$ are related by (non-standard) Laplace transform (47) just as $G(z)$ and $C(t)$ are. Equation (83) can be represented diagrammatically (here for the case $q=4$) by

[TABLE]

The dot represents a general SYK diagram (a fully-dressed 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 (83) has no closed form solution for general $q$. However, working with the power series representations, it enables the numerical calculation of $\mu_{2},\dots,\mu_{2n}$ in polynomial time and space complexity in $n$. Concretely, the following iteration algorithm can be easily implemented in a computer algebra system:

Set $g_{0}(z):=z^{-1}$, and let $j=0$. 2. 2.

Compute $c_{j}(t)$ from $g_{j}(z)$ by replacing $z^{-2n-1}$ with $(it)^{2n}/(2n)!$. 3. 3.

Set $\sigma_{j}(t):=c_{j}(t)^{q-1}$ up to order $t^{j}$. 4. 4.

Compute $\widetilde{\sigma}_{j}(z)$ from $\sigma_{j}(t)$ by replacing $(it)^{2n}$ with $z^{-2n-1}(2n)!$. 5. 5.

Set $g_{j+1}(z):=(1+J^{2}2^{2-q}g(z)\widetilde{\sigma}_{j}(z))/z$ up to order $t^{j+1}$. 6. 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)=\sum_{j=0}^{n}\mu_{2j}z^{-2j-1}$, which contains the correct moments up to $\mu_{2n}$. They can be then used to compute Lanczos coefficients $b_{1}^{2},\dots,b_{n}^{2}$ by the recipe (70). Arbitrary-precision rational number arithmetic is necessary for $n\sim 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 1 and Fig. 4 (a).

The above method can be readily adopted to variants of SYK where two-body and four-body interactions coexist:

[TABLE]

One only needs to replace the last term in (83b) by a sum over $q=2$ and $q=4$ with the corresponding coupling constants. Since the $q=2$ model is non-interacting, eq. (85) can be another model to study the effect of weak thermalizing interaction on the Lanczos coefficients. The results, shown in Fig. 9, are qualitatively consistent with those from the Ising model (Fig. 4): the linear growth rate depends only weakly on the interaction strength $J$ as it goes to zero. Quantitative, a logarithmic dependence

[TABLE]

describes the numerical data well for vanishing $J$.

B.2 Large-$q$ limit

In the large-$q$ limit, (83) can be solved analytically. It is convenient to define the coupling constant Maldacena and Stanford (2016); Roberts et al. (2018)

[TABLE]

It is then known Maldacena and Stanford (2016); Roberts et al. (2018) that $C(t)$ admits a $1/q$ expansion

[TABLE]

where the leading non-trivial term satisfies the following differential equation:

[TABLE]

whose solution is

[TABLE]

The corresponding moments

[TABLE]

where $(T_{n})_{n=0}^{\infty}=(1,2,16,272,7936,\dots)$ are the tangent numbers Sloane (2018a). The generating function of $T_{n}$ admits a continued fraction expansion Sloane (2018a):

[TABLE]

Using this, one can obtain the following Lanczos coefficients for the large-$q$ SYK model

[TABLE]

It is not hard to check using (8) that the wavefunction on the semi-infinite chain is

[TABLE]

The corresponding probability distribution is identical to the operator size distribution (see Eq. (5.11) of Ref. Roberts et al. (2018)):

[TABLE]

The large-$q$ SYK model is also solvable at any finite temperature Maldacena and Stanford (2016). The temperature $T$ is parametrized by $v\in(0,1)$ via

[TABLE]

The limits $T\to\infty$ and $T\to 0$ correspond to $v\to 0$ and $v\to 1$, respectively. The Lyapunov exponent is then

[TABLE]

and the autocorrelation under the Wightman inner product (60) is

[TABLE]

Comparing to (90), we see immediately that

[TABLE]

Therefore the finite-$T$ growth rate with the Wightman inner product is

[TABLE]

at any temperature. Thus, the bound $\lambda_{L,T}\leq 2\alpha_{T}^{(W)}$ is saturated at all temperature in the SYK model, whereas the bound $\lambda_{L,T}\leq 2\pi T$ is only so in the zero-temperature limit (see Fig. 8). Using the relation between growth rate and spectral function decay rate (15) and the relation (61) between spectral functions of different inner products, it is not hard to obtain the growth rate with the standard inner product from (100):

[TABLE]

Using (96), we obtain the limits $\alpha(T)_{S}\to\mathcal{J}\pi/2$ as $T\to 0$ and $\alpha_{T}^{(S)}\to\mathcal{J}$ as $T\to\infty$. We notice that $\alpha(T)_{S}$ increases at low temperatures while, in contrast, $\alpha_{T}^{(W)}$ decreases.

Appendix C Numerical Details for 1d Spin Chains

This section discusses the numerical details involved in computing the Lanczos coefficients and Krylov basis vectors in 1D spin chains. We work directly in the thermodynamic limit of a chain with $N\to\infty$ sites. However, bookkeeping will reduce this to finite-dimensional matrix multiplication.

Suppose we have a translation-invariant $k$-local Hamiltonian $H=\sum_{n}h_{n}$ and an $\ell$-local operator $\mathcal{O}=\sum_{n}\mathcal{O}_{n}$. Here $h_{n}$ and $\mathcal{O}_{m}$ are operators starting on sites $n$ or $m$ respectively. (For instance, we might have $\mathcal{O}_{2}=\cdots\otimes I_{1}\otimes X_{2}\otimes Z_{3}\otimes I_{4}\otimes\cdots$.) We normalize the operators so that $\mathinner{\left(h_{n}|h_{n}\right)}=1=\mathinner{\left(\mathcal{O}_{m}|\mathcal{O}_{m}\right)}$. At minor additional computational cost, we can work with an operator at a finite wavevector $q$:

[TABLE]

The crucial point is that applying the Liouvillian to $\mathcal{O}_{q}$ is another operator at wavevector $q$ by using translation-invariance to re-index the sum at the cost of phase factors. Explicitly,

[TABLE]

where

[TABLE]

where the shift is $s_{nm}$ is the index of the first non-identity site of $[h_{n},\mathcal{O}_{m}]$ minus $m$, which is needed to keep track of how much the support of the operator shifted due to the commutator. One can check that $\mathcal{O}_{m}^{\prime}$ starts on site $m$.

Therefore we only need to keep track of operators starting on a single site, say site [math]. We adopt the basis of Pauli strings and, following, e.g. Dehaene and De Moor (2003), we adopt a representation which minimizes the computational cost of taking commutators. Since $iY=ZX$, we may adopt a representation

[TABLE]

where $\delta,\epsilon,v_{k},w_{k}\in\left\{0,1\right\}$, i.e. a Pauli string of length $n$ may be represented by two binary vectors $\boldsymbol{v}$ and $\boldsymbol{w}$ of length $n$ and two binary digits. So if $\tau_{1}=i^{\delta_{1}}(-1)^{\epsilon_{1}}Z^{\boldsymbol{v}_{1}}X^{\boldsymbol{w}_{2}}$ and $\tau_{2}=i^{\delta_{2}}(-1)^{\epsilon_{2}}Z^{\boldsymbol{v}_{2}}X^{\boldsymbol{w}_{2}}$, then their commutator is a string $\tau^{\prime}=[\tau_{1},\tau_{2}]$ with

[TABLE]

All additions are performed over $\mathbb{Z}_{2}$.

With this setup, the Lanczos coefficients can be computed in a similar way to matrix-free exact diagonalization codes. A translation-invariant operator can be stored as a hash map of Pauli strings starting on site zero with complex coefficients. The Liouvillian is applied by combining (103), (104), and (106). Of course, it is not necessary to take $\mathcal{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 the 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 $\mathcal{O}_{n}$ grow in physical space. One way to characterize this is in terms of the distribution of string lengths in each $\mathcal{O}_{n}$. If $\mathcal{O}_{n}=\sum_{a}c_{a}\sigma^{a}$, where the sum runs over all Pauli strings $a$, then the distribution is defined by $P_{n}(s)=\sum_{a:\left|a\right|=s}\left|c_{a}\right|^{2}$. This distribution is shown for the Hamiltonian $H_{1}$ with the parameters given in Fig. 4. 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 (25) of the 1d quantum mechanics problem with Lanczos coefficients

[TABLE]

To solve this problem, notice that our infinite, tri-diagonal matrix is actually quite a familiar setup. If instead we had $b_{n}=\sqrt{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 the system described by $L$ has very high symmetry, due to an infinite-dimensional representation of the Lie algebra $\mathfrak{su}(1,1)$, enabling us to find an exact solution Hodges and Sukumar (2007); Sukumar and Hodges (2007). Indeed, there is a rich mathematical literature on the close connections between representations of $\mathfrak{su}(1,1)$, the combinatorics of Motzkin paths, and Meixner orthogonal polynomials Hetyei (2010); Viennot (1983). Our solution will be a simple application of these mathematical results.

We start with some generalities on orthogonal polynomials. Define $L_{(n)}=L_{0\leq i<n,0\leq j<n}$ to be the $n\times n$ matrix in the upper-left block of $L$. For example,

[TABLE]

We then define polynomials for each $n$ via

[TABLE]

By performing a cofactor expansion for the determinant on the $n$th row, the $Q$’s admit a three-term recursion relation

[TABLE]

together with initial conditions $Q_{0}(z)=1$ and $Q_{-1}(z)=0$. Eq. (110) should be compared with

[TABLE]

where $\left\{e_{n}\right\}$ is the natural orthonormal basis of $L$. In fact, (110) and (111) are equivalent, under the identification:

[TABLE]

Therefore, the polynomials $Q_{n}(z)$ are orthogonal, but not normalized. Instead they are monic, i.e., the highest order coefficient is unity: $Q_{n}(z)=z^{n}+\mathrm{O}(z^{n-1})$.

By construction, both $\left\{Q_{k}(z)\right\}$ and $\left\{z^{n}\right\}$ are a basis of $\mathbb{C}[z]$ and can be related by a triangular linear transform with matrix elements $\mu_{n,k}$:

[TABLE]

Combined with (112), and by orthonormality of $\left\{e_{n}\right\}$, we have

[TABLE]

and therefore

[TABLE]

The statements so far are general and apply to any set of Lanczos coefficients.

In the specific case $b_{n}=\sqrt{n(n-1+\eta)}$ (the extra overall factor $\alpha$ in (107) can be recovered by a simple time rescaling), one may recognize from the recursion relation (110) that $Q_{n}$’s are a special case of the Meixner polynomials of the second kind Sloane (2018b). They are a non-classical family of orthogonal polynomials defined by the following three-term recursion: Ismail (2009); Koekoek et al. (2010)

[TABLE]

In particular, $Q_{n}(z)=M_{n}(z;\delta=0,\eta)$. For these polynomials, the matrix elements $\mu_{n,d}$ have been exactly calculated, in terms of the following generating function Viennot (1983):

[TABLE]

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 Hetyei (2010); Ismail (2009):

[TABLE]

Now, taking $\delta=0$ and the series coefficient of $w^{d}$ in (117), we have

[TABLE]

Applying this to (115), and recalling $b_{n}=\sqrt{n(n-1+\eta)}$, we obtain the wavefunction solution

[TABLE]

where $(\eta)_{n}=\eta(\eta+1)\cdots(\eta+n-1)$ is the Pochhammer symbol. The general solution for $b_{n}=\alpha\sqrt{n(n-1+\eta)}$ can be obtained by a simple rescaling $t\mapsto\alpha t$, and is precisely Eq. (25) of the main text where, of course, $\mathinner{\left(\mathcal{O}_{n}|e^{i\mathcal{L}t}|\mathcal{O}_{0}\right)}=\mathinner{\left(e_{n}|e^{iLt}|e_{0}\right)}$. The special case $\eta=1$ of this family of solutions is well-known Viswanath and Müller (2008); Lee (2001). To the best of our knowledge, the general solution (119) has not been applied to the recursion method.

Appendix E Derivation of the q-Complexity Bound

This Appendix will derive Eq. (33), $\mathinner{\left(\mathcal{Q}\right)}_{t}\leq C\mathinner{\left(n\right)}_{t}$ for $C=2M$. The main idea of is that the definition of $\mathcal{Q}$ guarantees that the eigenbasis of $\mathcal{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 $\mathcal{Q}$ basis due to the dilation property. For any operator $\Phi$ where there is an $R>0$ such that $\mathinner{\left(q_{a}|\Phi\right)}=0$ for $q_{a}>R$, the hypothesis (28b) implies that $\mathinner{\left(q_{a}|\mathcal{L}|\Phi\right)}=0$ for $q_{a}>R+M$. Using (28c), as a base case for induction, we have $\mathinner{\left(q_{a}|\mathcal{L}^{n}|\mathcal{O}\right)}=0$ for $q_{a}>M(n+1)$ and, in particular, for $q_{a}>Cn$. By the construction of the Krylov basis,

[TABLE]

We claim that (120) implies

[TABLE]

for any operator wavefunction $\Phi$; taking $\Phi=\mathcal{O}(t)$, we obtain (33).

To show (121), we introduce projectors to large spectral values in the Krylov and $\mathcal{Q}$ bases, respectively:

[TABLE]

Then, we have for $n=q/C$,

[TABLE]

because $m<n=q/C\leq q_{a}/C$, $\mathinner{\left(q_{a}|O_{m}\right)}=0$ by (120). Equivalently,

[TABLE]

Applying this equation and its Hermitian conjugate, we have

[TABLE]

where the inequality follows from the fact that $\mathcal{P}^{Q}_{q}$ is a projector. Finally we need a standard integration-by-parts identity that converts the expectation value to an integral over the projectors:

[TABLE]

for any $k=1,2,3,\dots$. Combining the case $k=1$ and (124), we obtain

[TABLE]

which finishes the proof. More generally, for any $k$, we have

[TABLE]

This is useful as a bound on the growth rate of higher moments of the q-complexity super-operator. See Section VI.2 for an application.

Appendix F Geometric Origin of the Upper Bounds

In this appendix we derive the geometric upper bound for the Lanczos coefficients in one-dimensional quantum systems. The main object of our analysis will be the growth of the moments $\mu_{2n}=\mathinner{\left(\mathcal{O}|\mathcal{L}^{2n}|\mathcal{O}\right)}=\left|\left|\mathcal{L}^{n}\mathcal{O}\right|\right|^{2}$. Moments and Lanczos coefficients are equivalent, and Appendix A details how to translate between them.

To warm up, we first show a bound corresponding to linear growth (using essentially the same argument as in  Abanin et al. (2015, 2017)). This is asymptotically tight in $d>1$. Suppose we have a $2$-local Hamiltonian $H=\sum_{x}h_{x}$ and a $1$-local operator $\mathcal{O}$ (the general case of $r$-local $h_{x}$ and $r$-local $\mathcal{O}$ can be reduced to the previous case by a block renormalization step that groups consecutive sites into renormalized sites). The Liouvillian becomes a sum of terms $\mathcal{L}=\sum_{x}\ell_{x}$ with $\ell_{x}=[h_{x},\cdot]$. We suppose that the local terms are uniformly bounded, i.e., for all $x$, $\left|\left|h_{x}\right|\right|\leq\mathcal{E}$. Now, the moment $\mu_{2n}$ is the norm-squared of the sum

[TABLE]

This sum is highly constrained by the spatial structure of the spin chain. The operator $\mathcal{O}$ is supported only on one site, and the applications of the Liouvillian grow that support at the edges. Each term in (128) can be visualized as a discrete quantum circuit, where each gate $\ell_{x_{k+1}}$ must act on at least one site that is already in the support of $\ell_{x_{k}}\cdots\ell_{x_{1}}\mathcal{O}$ — otherwise the term vanishes due to the commutator. This condition is satisfied by at most $(k+1)\leq 2k$ positions $x_{k}$, so the total number of non-zero terms in (128) is at most $2^{n}n!$ for large $n$. The value of each non-zero term is itself bounded due to the finite local bandwidth $\mathcal{E}$, so $\left|\left|\ell_{x_{n}}\cdots\ell_{x_{1}}\mathcal{O}\right|\right|^{2}\leq(2\mathcal{E})^{2n}$. By the triangle inequality, we have

[TABLE]

By Stirling’s formula, the right hand side has the same asymptotics as (21), which corresponds to linear growth of the $b_{n}$’s. Hence (129) implies that the Lanczos coefficients can grow at most linearly in any dimension.

Notice that, the bound comes essentially from counting the number of sequences $x_{1},\dots,x_{n}$ that give rise to a nonzero contribution to (128). In what follows we show that, in one dimension, there is a sharper upper bound on this number, leading to the sub-linear growth announced in Section IV.3. For this, we suppose without loss of generality that $\mathcal{O}$ is supported on site [math] and $h_{x}$ on sites $x$ and $x+1$. Then it is not hard to see that a $\ell_{x_{n}}\cdots\ell_{x_{2}}\ell_{x_{1}}\mathcal{O}\neq 0$ only if for all $k=1,\dots,n$,

[TABLE]

We define $\mathcal{P}_{n}$ to be the set of $(x_{1},\dots,x_{n})$’s that satisfy (130) and denote its size by $P_{n}:=|\mathcal{P}_{n}|$. Then, similarly to (129), we have

[TABLE]

Hence bounding $\mu_{2n}$ reduces to bounding $P_{n}$, which is a completely combinatorial problem.

To produce this combinatorial bound, we partition the set $\mathcal{P}_{n}$ as follows

[TABLE]

Intuitively, if the support of the operator grows to size $\ell+1$ after $n$ applications of Liouvillian, then $(x_{1},\dots,x_{n})\in\mathcal{P}_{n,\ell}$. By “size”, we mean the distance between the endpoints, disregarding the “holes” between them. In the 1d case, the operator size can only grow in two places: the left and right sides. Therefore, for any $(x_{1},\dots,x_{n})\in\mathcal{P}_{n,\ell}$, $x_{k}=L_{k}$ or $x_{k}=R_{k}$ must hold for $\ell$ values of $k$ among $1,\dots,n$: for each of such $k$’s, one has only two choices for $x_{k}$. For the remaining $n-\ell$, there are (at most) $\ell$ choices (by (130), minus 2 boundary choices). Summarizing, we have

[TABLE]

where the binomial coefficient counts the choices of the $\ell$ values. Combining this with (132), we have

[TABLE]

In the limit $n\gg 1$, the maximum is attained at $\ell=n/W(n)$ where $W$ is the product-log function defined by $z=W(ze^{z})$. For large $n$, $W(n)=\ln n-\ln\ln n+o(1)$, so

[TABLE]

where we used $n/W(n)=e^{W(n)}$ and Stirling’s formula. Therefore

[TABLE]

which grows more slowly than the moment asymptotics corresponding to a linear growth with rate $\alpha$ (75), $B_{n}\ll\left(\frac{4n\alpha}{e\pi}\right)^{2n}$, for any $\alpha>0$. So the Lanczos coefficients corresponding to (136) must be sub-linear.

What, then, is the fastest possible growth of the $b_{n}$’s in 1D? Although we cannot bound the individual Lanczos coefficients in a useful way from the bound on the moments, we can use the bound on their cumulative product $\ln\prod_{k=1}^{n}b_{k}^{2}\leq\ln\mu_{2n}$ (73) and differentiate with respect to $n$. As a result, we find

[TABLE]

The bound (73) (together with (136)) is satisfied asymptotically by the above choice of $b_{n}$ if and only if $A\leq 4\mathcal{E}/e$. Therefore, $b_{n}=ae^{W(n)}$ captures the correct asymptotic behavior of the upper-bound in the moments, and qualifies as the maximal growth rate of Lanczos coefficients in 1d.