Euclidean operator growth and quantum chaos

We consider growth of local operators under Euclidean time evolution in lattice systems with local interactions. We derive rigorous bounds on the operator norm growth and then proceed to establish an analog of the Lieb-Robinson bound for the spatial growth. In contrast to the Minkowski case when ballistic spreading of operators is universal, in the Euclidean case spatial growth is system-dependent and indicates if the system is integrable or chaotic. In the integrable case, the Euclidean spatial growth is at most polynomial. In the chaotic case, it is the fastest possible: exponential in 1D, while in higher dimensions and on Bethe lattices local operators can reach spatial infinity in finite Euclidean time. We use bounds on the Euclidean growth to establish constraints on individual matrix elements and operator power spectrum. We show that one-dimensional systems are special with the power spectrum always being superexponentially suppressed at large frequencies. Finally, we relate the bound on the Euclidean growth to the bound on the growth of Lanczos coefficients. To that end, we develop a path integral formalism for the weighted Dyck paths and evaluate it using saddle point approximation. Using a conjectural connection between the growth of the Lanczos coefficients and the Lyapunov exponent controlling the growth of OTOCs, we propose an improved bound on chaos valid at all temperatures.


I. INTRODUCTION AND RESULTS
Operator spreading, or growth, in local systems is a question of primary interest, which encodes transport properties, emergence of chaos and other aspects of many-body quantum dynamics [1][2][3][4][5][6][7]. A classic result of Lieb and Robinson [8] (see also [9,10] for recent progress) establishes that under conventional time evolution the fastest possible spatial spreading of local operators is ballistic. There is no norm growth in this case since the time evolution is unitary. Ballistic spreading of operators, and signals, has been established for many models [3,[11][12][13][14][15] and seems to be a universal feature of local systems in any dimensions. At the same time, evolution of local operators in Euclidean time which we study in this paper, is much more nuanced.
Since the Euclidean evolution is not unitary, the norm of A(−iβ) quickly grows with β. Moreover, as we explain below, the operator growth is not universal and reflects if the system in question is integrable or chaotic. We start in section II by deriving a bound on |A(t)| valid uniformly for |t| = β by expanding (1) in Taylor series and bounding corresponding nested commutators. For local H there is a combinatorial problem of counting contributing nested commutators, which we solve exactly for short range systems defined on Bethe lattices, which includes local systems in 1D. In higher dimensions we conjecture an asymptotically tight bound. Hence, we expect our bounds on operator norm to be optimal in the class of Hamiltonians we consider -lattice Hamiltonians with local interactions. We find that maximal rate of growth is very different in 1D, where it is at most doubleexponential, and in higher dimensions or Bethe lattices, where the norm can become infinite in finite Euclidean time. We extend the analysis to include spatial growth in section III, where we find that in 1D operators spread at most exponentially, while in higher dimensions, including Bethe lattices, they can reach spatial infinity in finite Euclidean time. When the 1D system is finite, the minimal time necessary for an operator to reach the boundary is logarithmic, which may explain logarithmic convergence of the numerical Euclidean time algorithm proposed in [16]. We further speculate in section V that the timescale originating from the Euclidean Lieb-Robinson bound might be related to the Thouless energy of the corresponding quantum many-body system [17].
In section IV the results on norm growth are used to constrain individual matrix elements. We find that matrix elements in energy eigenbasis E i |A|E j must decay at least exponentially with ω = |E i − E j |, while in 1D the decay must be faster than exponential, as provided by (44) and (47). We also establish a number of bounds on the auto-correlation function at finite temperature C T (t), and its Fourier transform -the power spectrum Φ T (ω), ρ ∝ e −H/T , Tr(ρ) = 1.
The bounds have integral form, see (56,57) and (59). At the physical level of rigor, they suggest that Φ T (ω) decreases exponentially with ω in D ≥ 2, while in 1D the arXiv:1911.09672v1 [cond-mat.stat-mech] 21 Nov 2019 decay at large frequencies is superexponential. This emphasizes that one-dimensional systems are indeed very special, and many numerical results established for one dimensional systems may not necessarily apply to real higher-dimensional systems. The bound on |A(t)| established in section II depends only on the absolute value |t|. Obviously, it is overly conservative for real t when the time evolution is unitary. We argue, however, in section V that it does correctly capture the Euclidean growth t = −iβ of chaotic systems. We also consider system size dependence of |A(t)| and find it to be consistent with the Eigenstate Thermalization Hypothesis (ETH). For the integrable systems we find the growth of |A(t)| to be much slower than maximal possible, and in particular spatial growth of A(−iβ) in this case is not exponential but polynomial.
The bound on |A(−iβ)| can be translated into a bound on the growth of Lanczos coefficients b n , appearing as a part of the recursion method to numerically compute C T (t). This is provided we assume that asymptotically b n is a smooth function of n. To perform this calculation, we introduce a formalism of summing over weighted Dyck paths in section VI, and evaluate the corresponding path integral via saddle point approximation.
The obtained bound on Lanczos coefficients growth(75) is valid at all temperatures.
Translating it into a bound on Lyapunov exponent of OTOC, we find a new bound on chaos whereβ is such that C T (t) is analytic inside the strip | (t)| ≤β(T ). For local systems we findβ(T ) ≥ 2β * with β * given by (27) for all T . We illustrate this bound for SYK model in section VI, see Fig. 2.
We conclude with a discussion in section VII.

II. BOUND ON OPERATOR NORM GROWTH IN EUCLIDEAN TIME
Our goal in this section is to bound the infinity norm of a local operator evolved in Euclidean time Here f (β) is a bound which depends on the inverse temperature β, the strength of local coupling J and geometrical properties of the underline lattice model. We conjecture that our bound (15) (for 1D systems) and (26) (for higher dimensions) is optimal for the class of models characterized by the same strength of the local coupling constant J and lattice geometry encoded in the Klarner's constant λ and animal histories constant ε which we introduce later in this work.
For simplicity, we first consider nearest neighbor interaction Hamiltonian in 1D where each h I acts on sites I and I +1 and for all |h I | ≤ J for some J [18]. Any nearest neighbor interaction spin chain would be an example. The operator A will be an one-site operator. An example with L + 1 = 6 sites is shown in Fig. 1.
Euclidean time-evolved A(−iβ) can be expanded in Taylor series Using decomposition (5) operator A(−iβ) can be represented as a sum of nested commutators of the form [19] A Here the sum is over all sets of indexes {I 1 , . . . , I k } which satisfy the following "adjacency" condition: first index I 1 must be adjacent to the site of A, I 2 must be adjacent to the endpoints of I 1 (which include the site of A), I 3 is adjacent to the endpoints of the union of I 1 , I 2 , etc. In other words, any subset of bonds I 1 , I 2 , . . . , I for ≤ k defines a connected cluster. Otherwise, the commutator in (7) vanishes. A connected cluster of bonds of any particular shape is called a bond lattice animal. In 1D, all lattice animals consisting of j bonds are easy to classify: they are strings of consecutive bonds from some I to I + j. In higher dimensions, the number of different bond lattice animals of j bonds quickly grows with j.
Each set {I 1 , . . . , I k } from (7) defines a lattice animal, but not the other way around. This is because indexes can repeat and appear in different orders, subject to the constraints outlined above. If we think of the set {I 1 , . . . , I k } as a "word" written in terms of "letters" I , then corresponding lattice animal defines the alphabet.
There is a more nuanced characteristics of index sets from (7), the order in which new indexes appear. Namely, we take a set {I 1 , . . . , I k } and going from left to write remove indexes which have already appeared. In this way we obtain a new (shorter) set which satisfies the same adjacency condition described above. A particular order is called "history." For example, two sets {2, 3, 2, 4, 3} and {3, 3, 4, 2, 4} define the same lattice animal consisting of bonds I = 2, 3, 4 but different histories, {2, 3, 4} and {3, 4, 2} correspondingly, see Fig. 1.
Going back to the sum (7), to bound the infinity norm of A(−iβ) we can bound each nested commutator by (2J) k |A|. Then and the non-trivial task is to calculate the number of sets {I 1 , . . . , I k } for any given k, which satisfy the adjacency condition. Evaluating sum (9) can be split into two major steps. The first is to calculate the total number φ(j) of animal histories associated with all possible lattice animals consisting of j bonds. Second step is to calculate the sum over sets {I 1 , . . . , I k } associated with any given history {J 1 , . . . , J j }. This last problem can be solved exactly in full generality. Let's assume we are given a history -a set {J} = {J 1 , . . . , J j } which satisfies the adjacency condition. We want to know the number of different sets {I} = {I 1 , . . . , I k } for k ≥ j satisfying the adjacency condition such that {J} is the history of {I}. We denote this number by S(k, j). An important observation here is that any given set {I} defines a partition of {1, 2, . . . , k} into j groups labeled by elements from {J} by assigning each number 1 ≤ i ≤ k to a group specified by I i . And vice verse, each partition of {1, 2, . . . , k} into j groups defines a proper set {I} satisfying the adjacency condition. To see that we need to assign each group a unique label from {J}. We do it iteratively. The element 1 belongs to a group, which will be assigned the label J 1 . Then we consider element 2. If it belongs to the same group labeled by J 1 we move on to element 3, otherwise we assign the group it belongs label J 2 . Then we consider elements 3, 4 and so on. In this way all j groups will by labeled by the unique elements from {J} such that the adjacency condition is satisfied.
In other words, we have established a one-to-one correspondence between the space of proper sets {I} for the given history {J} with the space of partitions of k elements into j groups. The number S(k, j) of such partitions is the Stirling numbers of the second kind which admits the following representation [20] To evaluate (9) we still need to know the number of lattice animal histories for a given j. In case of 1D systems, those can be calculated exactly, while in higher dimensions we propose an asymptotically tight bound. Hence, we consider these cases separately.

A. 1D systems
In one dimension, all lattice animals consisting of j bonds are simply the strings of j consecutive bonds. There are N (j) = j + 1 such animals which include the site of the operator A. A convenient way to enumerate them is to count the number of bonds j 1 and j 2 , j 1 + j 2 = j, to the left and to the right of A, respectively. For the given j 1 , j 2 there is, obviously, only one animal, N (j 1 , j 2 ) = 1.
For any given j 1 , j 2 we denote by h(j 1 , j 2 ) the number of histories associated with this animal, i.e. the number of different sets {J} = {J 1 , . . . , J j } such that each J i belongs to the animal, all J i in the set are unique and {J} satisfies the adjacency condition. Each history {J} can be completely parametrized by the order in which the cluster "grew" in left and right directions, for example histories {2, 3, 4} and {3, 4, 2} from Fig. 1 can be parametrized as "left,right,right" and "right,right,left" correspondingly. In other words histories with given j 1 , j 2 are in one to one correspondence with strings of j elements, each element being either "left" or "right," and there are in total j 1 and j 2 elements of each kind. Obviously, the total number of such strings is Combining all ingredients together, we find from (9) By definition k ≥ j. Crucially, expression (10) vanishes for 1 ≤ k < j. Therefore the sum over k can be extended to go from k = 1 and can be easily evaluated, where j := j 1 + j 2 . The sum over s can be evaluated explicitly, Here f (j 1 , j 2 , β) is the bound on the norm of the part of A(−iβ) supported on the cluster of size j 1 + j 2 . We can sum over different j 1 , j 2 corresponding to the total cluster size, Summing over all different cluster sizes j gives the full bound Re-expanding (15) intro Taylor series in β, where B k are Bell polynomials, we find a bound on the norm of individual nested commutators, These results are trivial to generalize to the case of an infinite 1D lattice with a boundary.

B. Higher dimensional systems
The calculation of previous section can be extended to an arbitrary lattice system. If we fix a particular animal lattice history {J 1 , . . . , J j } and sum over all possible sets {I 1 , . . . , I k }, k ≥ j, associated with this history with the weight (2J|β|) k /k!, by repeating the steps of the previous section we find [21] k≥j {I1,..., If we now introduce the total number of lattice animal histories φ(j), associated with all possible lattice animals consisting of j bonds, then f (β) defined in (9) is trivially given by The number of different lattice animal histories is difficult to evaluate exactly. But it is known that the number of different lattice animals consisting of j bonds (which include a particular site) grows rapidly in higher dimensions. While the exact formula is not known, the asymptotic growth is known to be exponential with growth controlled by the so-called Klarner's constant λ, By introducing a sufficiently large but j-independent constant C we can uniformly bound the number of lattice animals consisting of j bonds by [22] The number of histories for any given animal is the number of different sets {J 1 , . . . , J j } where all indexes are distinct, subject to the adjacency condition. Let us denote by h(j) the average number of histories for all animals consisting of j bonds. Then, it is trivially bounded by h(j) ≤ j!. It can be shown that for sufficiently large j [19] for some a > 1. We, therefore, conjecture that for local systems h(j) is uniformly bounded by for some ε > 1 and a j-independent constant C ≥ 1. This bound is trivially satisfied for ε = 1. The nontrivial part here is the expectation that (23) correctly captures the leading (exponential) asymptotic behavior of h(j) with some ε > 1, i.e. (23) is the optimal bound which can not be further improved. We therefore introduce here the constant ε which we call animal histories constant and conjecture that it is strictly larger than 1.
In the end of this section we also derive a lower bound on ε/λ. By combining (21) together with (23) we find the bound Here f (β) is defined to be larger than the sum in (9). The coefficient characterizes lattice geometry. Unlike in 1D, where (14) has an additional factorial suppression factor, f (j, β) in higher dimensions grows exponentially for sufficiently large β. Summing over j yields In contrast to 1D, while (15) is finite for all β, (26) is finite only for While (26) is only a bound on f (β) defined in (9), location of the singularity in both cases is the same because it is only sensitive to the asymptotic behavior of N (j) and h(j). Expanding (26) in Taylor series where P k are the polynomials defined via yields a bound on individual nested commutators The divergence of bound (26) at |β| = β * is not an artifact of an overly conservative counting, as confirmed by a 2D model introduced in [19], for which |A(−iβ)| is known to diverge. We will argue in section V that the growth outlined by the bounds (15,26) reflects actual growth of |A(−iβ)| in non-integrable systems and singularity of (26) at finite β is a sign of chaos. We also note that in case of 1D systems the bound (15) ensures that the operator norm of A(t) remains bounded for any complex t. This is consistent with analyticity of correlation functions in 1D [23]. On the contrary, in higher dimensions, physical observables may not be analytic. We discuss the relation between the singularity of |A(−iβ)| and non-analyticity of physical observables due to a phase transition in section V and show that they have different origin.
To conclude this section, we point out the advantage of counting lattice animal histories in (19) and derive a bound on ε/λ. There is a straightforward way to estimate the number of sets I 1 , . . . , I k in (9) from above by counting the number of ways a new bond can be added to the set at each step. Provided the lattice has coordination number z, starting from the site of A, there are z ways to choose I 1 , at most z(2z) ways to choose I 2 , z(2z)(3z) ways to choose I 3 and so on. As a result we would have an estimate for f (β), This result was previously obtained in [24,25]. This gives the following estimate for the location of the pole Since the bound above ignores the adjacency conditions, it can not be better than true f (β), and in particular the location of the singularity (32) can not exceed β * defined in (27). This yields a constraint ε/λ ≥ e 1/z − 1.
A result analogous to (26) has been previously established in [26]. There q 0 was lattice animal constant, which is Klarner's constant introduced in previous section, q 0 = λ. Crucially, we improve this result to account for proper lattice animal histories by introducing animal histories constant ε in (33). As we have shown above without ε the critical value of β where f (β) diverges may actually be worse than the approximation (32).

C. Bethe lattices
We saw above that the behavior of f (β) differs drastically in one and higher dimensions. To better understand this difference we consider an "intermediate" scenario of a short range Hamiltonian define on a Bethe lattice of coordination number z [27]. Namely, we assume that each h I from (5) "lives" on a bond an acts on the Hilbert spaces associated with the two vertexes adjacent to that bond. For any finite k in the Taylor series expansion (7) only finite number of bounds are involved and the corresponding lattice animals (clusters) live on the Cayley tree. Thus, similarly to 1D, there are no loops, but the total number of lattice animals consisting of j bonds grows exponentially, (34) as in higher dimensions.
The total number of lattice animal histories φ(j) can be calculated exactly in this case (see the appendix A), leading to the bound In other words, the total number of histories grows as a factorial, exactly as in higher dimensions. It is also interesting to compare (35) with (24), which yields lattices animal histories constant ε for Bethe lattices, It can be easily verified that exact ε and λ satisfy (33) and ε > 1 for any z ≥ 2, supporting our conjecture that z is always strictly larger than 1.
As a final remark, we notice that taking z → 2 in (36) yields f (β) = e 2q , in full agreement with (15). In this limit β * becomes infinite while the approximation (32) remains finite.

III. SPATIAL GROWTH IN EUCLIDEAN TIME
While deriving the bound on the norm of local operators evolved in Euclidean time, (15) and (26), we obtained a stronger result -a bound f (j, β) (or f (j 1 , j 2 , β) in 1D) on spatial growth of A(−iβ). It can be immediately translated into the Euclidean analog of the Lieb-Robinson bound [8] on the norm of the commutator of two spatially separated local operators. If B is an operator with finite support located distance away from A (measured in the Manhattan norm in case of a cubic lattice), then in D ≥ 2 where we assumed that |β| < β * . For larger |β| there is no bound as the sum does not converge. This result means that the local operator can spread to the whole system, no matter how large or even infinite that is, in finite Euclidean time β = β * . We will argue in section V that this is the true physical behavior in the chaotic case and therefore the bound can not be improved to get rid of the divergence at |β| = β * in full generality.
In 1D the situation is very different. Assuming local operator B is located bonds away from A we find (38) (If the system is infinite only in one direction and A is sitting at the boundary, one factor of e q should be removed.) Qualitatively the RHS of (38) behaves as for q +1, and asymptotes to 2|A||B|e 2q for q +1. This means a local operator spreads exponentially fast, to distances ∼ e 2Jβ , in Euclidean time β.
Exponential spreading of operators in 1D seems to be in agreement with the convergence of the Euclidean variational algorithm of [16] in logarithmic time. The connection between Euclidean Lieb-Robinson bound and the convergence time is intuitive, but difficult to establish rigorously, in particular, because the latter is sensitive to the choice of initial wave-function. For the integrable models, for which the spreading of operators is at most polynomial (see section V), convergence time might be even shorter because of a well-tuned initial wave-function. For the chaotic systems we expect no fine-tuning of the initial state and hence a direct relation between the convergence time and Euclidean Lieb-Robinson bound.
Another possibly intriguing connection is with the studies of Thouless times in chaotic Floquet systems without conserved quantities [17]. There, it was noticed that in 1D Thouless time is logarithmic in system size (see also [28]), and finite in D ≥ 2 (see, however, [29]). That is exactly the same behavior as in the case of Euclidean operator spreading. One potential interpretation would be that Thouless time can be associated with the slowest Euclidean mode propagating in the system. Under Euclidean time evolution with a time-dependent random Hamiltonian our extension of Lieb-Robinson bound holds. We also surmise that in this case spatial growth of all quantities, including the slowest, is qualitatively and outlined by the bound with some effective J, q 0 . When the system in question has a local conserved quantity, the slowest transport mode is diffusive, leading to L 2 scaling of Thouless time [30]. Thus, to compete this picture it would be necessary to establish that under Euclidean time evolution time necessary for a diffusive mode to travel across the system is the same as in the Minkowski case, i.e β ∼ L 2 , where L is the system size.
Finally, we notice that the Euclidean analog of the Lieb-Robinson bound in 1D (38) looks similar to the conventional Minkowski bound [10] with 2Jt substituted by q(β).

A. Individual matrix elements
Constraints on the infinity-norm of A(iβ) † = A(−iβ) provide an upper bound on the magnitude of matrix elements A ij = E i |A|E j in the energy eigenbasis. Starting from This inequality holds for any β and we therefore can optimize it over β. Using explicit form of the bound (15) in 1D we find optimal value of β to be (without loss of generality we assumed ω = E i − E j ≥ 0) This yields where These results shows that in 1D for large energy difference ω = |E i −E i | J off-diagonal matrix elements A ij decay faster than exponential. For ω ≤ 4J the bound trivializes to |A ij | ≤ |A|.
In higher dimensions the bound on A ij from (42) can not be better than exponential. This is because f (β) is a monotonically increasing function of β which diverges for some |β| = β * . In particular for any β and ω ≥ 0. To find leading exponent we optimize (42) over β to find, and Taking ω → ∞ limit, we find that the asymptotic exponential behavior is given by (45), where C is some ω-independent constant. Constraints on individual matrix elements (45) and (45) only depend on energy difference ω. In the case when the system satisfies ETH, off-diagonal matrix elements for i = j are known to be exponentially suppressed by the entropy factor, |A ij | 2 ∼ e −S . Therefore for the chaotic systems the bound will be trivially satisfied unless ω is extensive.

B. Constraints on power spectrum
Bounds on individual matrix elements found above can be extended to the autocorrelation function of a Hermitian local A, and its power spectrum Here ρ is an arbitrary density matrix which commutes with the Hamiltonian, ρ = i p i |E i E i |, Trρ = 1.
Although bounds on moments M k derived below are universal for all ρ, in what follows we will be most interested in the case when ρ is the Gibbs ensemble ρ = e −H/T /Z, in which case autocorrelation function and power spectrum will be denotes by C T and Φ T correspondingly. As a function of complex argument C T satisfies, First we notice that for any complex t, which guarantees analyticity of C(t) for 1D systems on the entire complex plane.
Using the bound on individual nested commutators (17) and (30) one can bound the growth of Taylor coefficients of C, To obtain an optimal bound, nested commutators should be split equally between two A's using cyclicity of trace Here R k = B k (2) for infinite 1D system, R k = B k for semi-infinite 1D system with a boundary, and R k = C P k (q −1 0 ) for D ≥ 2. Using the asymptotic behavior of Bell polynomials [31] B n (x) ∼ n 1 + o(1) e log(n/x) n , n x, and the Stirling approximation formula, the bound on moments for k 1 can be rewritten as (for the infinite 1D system) It is easy to see that the Taylor series of C T (t) converges in the whole complex plane, as was pointed out above. In higher dimensions, to find asymptotic behavior of P k (x) for large k, we use the following representation Substituting the sum over j by an integral and taking saddle point approximation gives Focusing on the case when ρ = e −H/T /Z, (57) guarantees that Taylor series of C T (t) converges absolutely inside the disc |t| ≤ 2β * . By representing C T as a sum over individual matrix elements it is easy to see that if the sum for C T (−iβ) is absolutely convergent, then it is absolutely convergent for any C T (t), (t) = −iβ. Therefore C T (t) is analytic inside the strip 2β * > (t) > −2β * . Because of reflection symmetry (50) function C T (t) must be analytic inside a wider strip 2β * > (t) > −2β * −1/T [32]. Hence symmetrically ordered autocorrelation function is analytic inside the strip 2β * + 1/(2T ) > (t) > −2β * − 1/(2T ), which is wider than the strip of analyticity of C T (t), and indicates a more rapid exponential decay of the power spectrum Φ W T of (58) in comparison with Φ T (ω).
The logic above is general and does not require any specific details of M k . Using reflection symmetry (50) we have shown in full generality that if C T (t) develops a singularity at t = ±i2β * , then C W T (t) is analytic at least inside the strip | (t)| ≤ 2β * + 1/(2T ).
There is another integral bound on power spectrum, valid for any density matrix ρ which commutes with H. By integrating (49) we find the following inequality Here ω is non-negative and in the second equality we used (41) with an arbitrary positive β. Now we can use |A(−iβ)| ≤ |A|f (β) and optimize over β, yielding Function κ is given by (45) and (47) in D = 1 and D ≥ 2 correspondingly. When ρ is maximally mixed state, i.e. temperature T is infinite, the bound can be strengthen to We would like to emphasize that all bounds discussed above, i.e. bounds on M k and (59), are integral in form. We do not know a rigorous way to directly constrain asymptotic behavior of Φ(ω). At the same time at physical level of rigor, if we assume that Φ(ω) is a smoothly behaving function at large ω, analyticity of C(t) inside the strip | (t)| < 2β * immediately implies that power spectrum in D ≥ 2 is exponentially suppressed by In 1D we similarly find super-exponential suppression The bound on moments for large k (56,57) and the integral bound (59) for large ω follow from here via saddle point approximation. Superexponential suppression of Φ(ω) emphasizes peculiarity of one-dimensional systems. In particular, it implies that high frequency conductivity [33] and energy absorption [24] for such systems will be superexponentially suppressed. This is a very special behavior, which should be kept in mind in light of the numerical studies, which are often limited to one dimensions, and therefore may not capture correct physical behavior.
An exponential bound on the integral of Φ(ω) was first established in [24], where the authors also noted superexponential suppression in 1D, albeit without proposing an explicit analytic form.

V. FINITE SIZE SCALING AND CHAOS
The bounds obtained in section II correctly account for the number of non-trivial nested commutators [h I k , [. . . , [h I1 , A]]] but do not take into account peculiarities of individual local Hamiltonians h I . We therefore expect our bound to be strongest possible among the uniform bounds for the entire family of local shortranged Hamiltonians defined on a particular lattice. We further assumed that each nested commutator is equal to its maximal possible value (2J) k |A|. This is certainly too conservative, but for the chaotic systems, i.e. in absence of some additional symmetries, we expect a finite fraction of nested commutators to grow as a power of k. We therefore expect that for large β our bounds (15) and (26) to qualitatively correctly describe actual Euclidean growth in local chaotic systems with some effective values of J, as it happens in [19]. In particular in one dimensions we expect |A(−iβ)| to grow double-exponentially, and in higher dimensions we expect |A(−iβ)| to diverge at some finite β * .
We similarly expect the bound on spatial growth outlined in section III to correctly capture the spread of local operators when the system is chaotic. An indirect evidence to support that comes from the numerical results of [16], i.e. logarithmic convergence time of a numerical Euclidean time algorithm, in agreement with (39).
Below we further outline how |A(−iβ)| reflects chaos of the underlying system when the system size is finite. It follows from (19) that for large β animal histories with the largest number of bonds will dominate, where j is the total number of bonds in the system, i.e. j is proportional to the volume. Let us compare this behavior with the growth of the Frobenius norm, At large β leading behavior is where ∆E is the maximal value of ∆E = E i − E j such that corresponding matrix element A ij is not zero. (In other words ∆E is the support of Φ(ω).) For the chaotic systems satisfying Eigenstate Thermalization Hypothesis we expect most matrix elements to be non-zero, even for extensive ∆E, matching extensive behavior of 2Jj in (63). Assuming qualitative behavior of (26) is correct for non-integrable systems, going back to thermodynamic limit in D ≥ 2, we expect singularity of |A(−iβ)| and C(−2iβ) at some finite β. This singularity has a clear interpretation as a manifestation of an ergodic behavior of A. We first interpret (A|B) := Tr(A † B)/Tr(1) as a scalar product in the space of operators and denote corresponding Frobenius norm of A by |A| F ≡ (A|A) 1/2 . Then if A were typical, i.e. Haar random, in the space of all operators, Euclidean time evolution can be split into two parts, A(−i(β + β )) = e β H A(−iβ)e −β H such that At time β = 0 we start with a local operator, which is not typical. In principle A(−iβ/2) only explores a particular trajectory in the space of all operators, and therefore can not be fully typical at any β. Yet, if we assume that at time β the trajectory of A becomes ergodic and A(−iβ/2) can be considered typical enough, we obtain Taking into account that free energy ln(Z) is extensive, we immediately see that (67) diverges for any β > 0. Hence, the singularity of C(−iβ) and thus also of |A(−iβ/2)| marks the moment when A(−iβ/2) becomes typical.
It is interesting to note that since C(t) is analytic for local one-dimensional systems, for such systems, even non-integrable, A never becomes typical and hence these systems can not be regarded as fully chaotic.
We separately remark that the conventional time evolution C(t) = (A(t/2)|A(−t/2)) does not have an interpretation as the Frobenius norm-squared of A(t/2), therefore (66) does not apply and even if A(t/2) becomes sufficiently typical at late t, the analog of (67) may not hold.
If the system is finite, at large β free energy simply becomes ln Z(β)/Z(0) ∼ −βE m , where E m is extensive (minimal or maximal ) energy of the system. Hence (67) will be proportional to e β ∆E , where ∆E is extensive, in full agreement with (65). This gives the following qualitative behavior of C(−iβ) when the chaotic system is sufficiently large but finite. For small β, ln C(−iβ) will behave as ∝ e q in 1D and ∝ ln(q 0 −q) in higher dimensions. This growth will stop at β ∼ log(L) in 1D or β ∼ β * in D ≥ 2, at which point in both cases ln C(−iβ) will be extensive. At later times ln C(−iβ) will grow as β∆E with some extensive ∆E. In the non-integrable case the transition between two regimes, "thermodynamic" when C(−iβ) has not yet been affected by the finite system size, and "asymptotic," is very quick, at most doublelogarithmic in 1D.
Behavior of chaotic systems described above should be contrasted with integrable models. In this case most matrix elements A ij are zero and for a wide class of systems, including classical spin models and systems with projector Hamiltonians, support of Φ(ω) remains bounded in the thermodynamic limit. (In terms of the Lanczos coefficients, introduced in the next section, this is the case of λ = 0.) For such systems the bounds (15) and (26) will be overly conservative. For sufficiently large systems and large β we expect (65) with a system size independent ∆E. This asymptotic behavior will emerge in finite system-independent Euclidean time. Infinity norm |A(−iβ)| will behave similarly. We further can use ( where last step assumes β∆E . This bound has the same structure as the conventional Lieb-Robinson bound in Minkowski space (40). Thus, in the case of non-interacting models or projector Hamiltonians (λ = 0 in the language of next section) we find ballistic spreading of operators for any complex t.
In the case of a general integrable model, the support of Φ(ω) is extensive and the behavior is more intricate. In many explicit examples in the thermodynamic limit Φ(ω) decays as a Gaussian, and C(−iβ) ∝ e (Jβ) 2 with some appropriate local coupling J [34][35][36][37][38]. (This is the case of λ = 1 in terms of the next section. See appendix B where we derive the Gaussian behavior starting from λ = 1.) Using the same logic as above this leads to the Euclidean Lieb-Robinson bound of the form which indicates a polynomial propagation of the signal ∝ β 2 . For a finite system of linear size L we may expect Gaussian behavior C(−iβ) ∝ e (Jβ) 2 up to the times β ∝ L 1/2 , after which the asymptotic behavior (65) should emerge. Although the model is integrable, ∆E is extensive, which implies the transition between "thermodynamic" and "asymptotic" behavior is long and will take up to β ∼ L. This indicates the qualitative difference between integrable and non-integrable (chaotic) models. When the system is finite in both cases the asymptotic behavior is given C(−iβ) ∝ e β∆E with an extensive ∆E (except for the λ = 0 case), but asymptotic behavior will emerge quickly, in finite (for D ≥ 2) or logarithmic (for D = 1) times in the non-integrable case, while in the integrable case asymptotic behavior will emerge much slower, after polynomial times in L.
A qualitatively similar picture will also apply if integrability is broken weakly, by a parametrically small coupling. For an operator initially characterized by λ = 0, the correlation function will first exhibit (65) with some sub-extensive ∆E, which will gradually grow to extensive values. It would be interesting to study this transition in detail, to see if the required times may be parametrically longer than β ∼ L.
We stress that non-analyticity of C(t) at imaginary times is due to A(−iβ) becoming typical and is not related to non-analyticity of free energy ln Z(β) due to a phase transitions at some temperature β. Indeed, C(t) for the SYK model is known to have a pole at imaginary time [39], while there is no phase transition and free energy is analytic. On the contrary, for the 3d Ising ln Z(β) is non-analytic due to a phase transition, but C(t) is entire, simply because A(t) explores only a very small part of the corresponding Hilbert space.
In conclusion, we note that the singularity of |A(−iβ)| and C(−iβ) at finite β in the thermodynamic limit has an IR origin. A straightforward attempt to extend the analysis of this section to field theoretic systems, which can be obtained from lattice systems via an appropriate limit, fails because both |A(−iβ)| and C(−iβ) are UV-divergent, and this obscures the IR divergence due to chaos. Formulating the criterion of chaos for QFTs using Euclidean operator growth thus remains an open question.

VI. CONSTRAINTS ON LANCZOS COEFFICIENTS
The bound on power spectrum established in section IV B can be used to constrain the growth of Lanczos coefficients. To remind the reader, Lanczos coefficients b n are non-negative real numbers associated with an orthonormal basis in the Krylov space A n generated by the action of H on a given operator A 0 = A. Starting from a scalar product and choosing A normalized such that |A| 2 = (A, A) = 1, Lanczos coefficients are fixed iteratively from the condition that operators A n defined via An autocorrelation function C W = (A(t), A), defined via scalar product (69), can be parametrized in a number of ways, via its power spectrum Φ W (ω), Taylor coefficients (moments) M k , or Lanczos coefficients b n . Schematically an asymptotic growth of b n for large n 1 is related to the behavior of M k , k 1, high-frequency tail of Φ W (ω), ω → ∞, or growth of C W (t) at the Euclidean time t = −iβ. But the detailed relation is not always trivial. Assuming exponential behavior of power spectrum at large frequencies it is trivial to obtain the growth of M k and C W (β) by calculating corresponding integrals over ω using saddle point approximation. Although much less trivial, but starting from the power spectrum (70), it is also possible to establish an asymptotic behavior of Lanczos coefficients [40] b 2 n ∝ n λ .
The converse relations between asymptotic behavior of b n , M k , Φ W (ω) and C W (−iβ) are much more subtle and may not hold. Thus, we show in the appendix C that smooth asymptotic behavior of M k does not imply smooth asymptotic of b n . It was proposed long ago that λ defined in (71) falls into several universality classes, characterizing dynamical systems [36]. In particular it was observed that λ = 0 for non-interacting and λ = 1 for interacting integrable models. (It should be noted that since λ characterizes a particular operator, the same system may exhibit several different values of λ.) Recently it was argued in [7] that λ = 2 is a universal behavior in chaotic systems in D ≥ 2.
To thoroughly investigate possible implications of this conjecture, it is desirable to derive the constraints on the behavior of C W , Φ W , and M k starting directly from the assumption that b n is a smooth function of n for large n. In full generality Lanczos coefficients b n are related to the moments M k via Here the sum is over Dyck paths parameterized by the sets satisfying h 0 = h k = 1/2, and h i+1 = h i ± 1, h i > 0. Assuming (71) out goal is to deduce an asymptote of M 2k using (72). We develop the approach of summing over weighted Dyck paths in the appendix B. Here we just mention main results. If b 2 n is asymptotically a smooth function of n, path integral over Dyck paths can be evaluated via saddle point approximation by identifying a trajectory in the space of indexes, which gives the leading contribution. Thus if b n is smooth, M k is also smooth. Furthermore, if λ = 2, b 2 n ∼ α 2 n 2 , and the leading order behavior is Thus, starting from the asymptotic behavior b 2 n ∝ α 2 n 2 we necessarily find that C W has a singularity at β = π/(2α), in full agreement with the conjecture of previous section that singularity in Euclidean time is the characteristic property of chaos.
Provided C W (t) is analytic inside a strip (t) ≤β W for someβ W would immediately imply a bound When ρ ∝ e −H/T , provided autocorrelation function C T (48) is analytic inside | (t)| ≤β(T ), function C W T defined in (58) will b analytic at least inside | (t)| ≤β W = β(T ) + 1/(2T ) (see the discussion in section IV B) and therefore We have also established in section II B thatβ(T ) ≥ 2β * for all T . The coefficient α has been recently conjectured to bound maximal Lyapunov exponent governing exponential growth of the out of time ordered correlation function (OTOC) [7,41], λ OTOC ≤ 2α. This leads to the improved bound on chaos which is stronger than the original bound λ OTOC ≤ 2πT of [42]. In the limit of quantum field theory, (β * ) −1 will be of the order of UV-cutoff, reducing (76) to the original bound. Yet the new bound is non-trivial for discrete models exhibiting chaos.
To illustrate the improved bound, we plot (76) in Fig. 2 for the SYK model in the large q-limit [43] against the exact value of λ OTOC , evaluated in [7,39]. We take 2β * = 1 to ensure that the autocorrelation function C T is analytic inside (t) < 2β * = 1 for all T . Temperature T is parametrized via 1 ≥ v ≥ 0, πvT = cos(πv/2) such that the exact Lyapunov exponent is λ OTOC = 2 cos(πv/2). (77) We have emphasized above that for 1D systems with short range interactions C T (t) has to be analytic in the entire complex plane. This imposes a bound on the growth of Lanczos coefficients. Assuming b n is a smooth function of n [7] proposed that the asymptotic growth in 1D non-integral systems will acquire a logarithmic correction b n+1 ≈ α n log(n/n 0 ) .
Using the integral over weighted Dyck paths in the appendix B, we find this to be consistent with the behavior of M k outlined in (56) provided Sum over Dyck paths in the case of λ = 1 associated with integrable systems is discussed in the appendix B. Since for the local models C W (t) is analytic inside a sufficiently small vicinity of t = 0, asymptotic behavior with λ > 2 in such systems is excluded.

VII. CONCLUSIONS
We have derived a number of rigorous bounds on the infinity norm of a local operator evolved in Euclidean time, and extended them to autocorrelation function (2). The novel ingredient of our approach is the counting of lattice animal histories and formula (19), using which we solved exactly combinational problem of counting nested commutators for Bethe lattices (and establish acorrect asymptotic for lattices in D ≥ 2). Some of the bounds derived in this paper were known before. We improved numerical coefficients, including the location of the singularity β * in D ≥ 2. Our results are strongest possible among the bounds uniformly valid for all local Hamiltonians characterized by the same |h I | ≤ J defined on a lattice of a particular geometry.
We have also established Euclidean version of Lieb-Robinson bound on the spatial operator growth. In 1D operators spread at most exponentially, while in D ≥ 2 operators can reach spatial infinity in finite Euclidean time. When the system is integrable, in all D operators spread polynomially.
As a main point of this paper, we advocated that Euclidean operator growth reflects chaos in the underlying quantum system. If the system is chaotic, the norm growth and spatial growth are maximal possible and the operator norm diverges at some finite Euclidean time. We interpreted this divergence as a consequence of ergodicity in the Krylov space.
There are several distinct characteristic properties of chaos for many-body quantum systems. One is the Eigenstate Thermalization Hypothesis [44], which is concerned with individual matrix elements. Another popular probe is out of time ordered correlation function, which extends the notion of exponential Lyapunov growth to quantum case. Its use as a characteristic of many-body quantum chaos was pioneered in [45][46][47] and brought to the spotlight by applications to quantum gravity [42]. Despite recent efforts [41,[48][49][50] there is no clear understanding of how to relate these two characteristics of chaos to each other. We hope that the Euclidean growth, which on the one hand is related to ETH via the behavior of C(−iβ) at large β, see (65), and on the other hand is related to OTOC via the bound (76), may provide such a bridge.
While this is not necessary for the bound on operator norm growth, for completeness we derive the number of lattice animals consisting of j bonds, N (j). We first consider all lattice animals which originate at the same vertex and extend into one particular direction ("branch") on the Bethe lattice. If the number of such animals is n(j), then it must satisfy the recursive relation n(j) = j1+...jz−1=j−1 n(j 1 ) . . . n(j z−1 ). (A6) It reflects the fact that we can "move" the initial point by one bond inside the branch and decompose j into j = z i=1 j i , j z = 1. In (A6) we also use that n(1) = 1. This gives in full generality .
The full number of lattice animals N (j) is related to n(j) via with the total number being . In the context of recursion method Lanczos coefficients b n define tri-diagonal Liouvillian matrix such that correlation function where scalar product of operators is defined in (69), and 0| . . . |0 denotes the upper left corner matrix element. By definition, moments M k are Taylor series coefficients of C W , From here and the tri-diagonal form of L it follows that while all odd moments vanish (this is also obvious from the symmetry C W (t) = C W (−t)). The sum above is over the sets h i such that h 1 = h 2k = 1/2, h i > 0, and h i+1 = h i ± 1. When k is large, sum over Dyck paths becomes a path integral, parametrized by a smooth function f (t), 0 ≤ t ≤ 1 [52], Function f (t) satisfies Derivative f (t) defines an average slope of a "microscopic" Dyck path around index i ≈ 2kt. The path is a sequence of "up" and "down" jumps with the probabilities p and 1 − p, which vary smoothly, such that 2p(t) − 1 = f (t). The number of different "microscopic" paths N [f (t)] associated with f (t) is given by the Shannon entropy of p(t), In other words N [f (t)] is the measure in the path integral over f (t). To verify this result we calculate the total number of Dyck paths, which is known to be given by the Catalan number, by evaluating corresponding path integral via saddle point approximation. By interpreting S 0 [f (t)] as a classical action, classical EOM is The only solution satisfying boundary conditions is f (t) = 0, which gives saddle point value This reproduces correct exponential behavior of Catalan numbers, C k ≈ 4 k /(k 3/2 π 1/2 ). Assuming b n+1 = b(n) is a smooth function of index, at least for large n, sum over weighted Dyck paths (B4) can be represented as an integral In case of asymptotic behavior b 2 (n) = α 2 n λ the EOM is . (B12) For general λ this equation can be solved in terms of an inverse of the Hypergeometric function. We are most interested in two cases, λ = 2 and λ = 1. In the latter case, b 2 (n) = α 2 n, the saddle point trajectory is f (t) = t(1 − t) leading to the asymptotic behavior of moments This gives an exponential growth of C W at larger β, In the "chaotic" case λ = 2 the solution satisfying boundary condition is and the saddle point value is Provided C W is analytic inside the strip (t) <β, the asymptotic growth constant has to be bounded by α ≤ π 2β . (B17) Finally we consider the scenario when the growth of Lanczos coefficients acquires logarithmic correction, b(n) = α n log(n/n 0 ) .
In other words at leading order the effect of logarithmic correction is in adjusting the scaling parameter λ. When λ ≈ 2, the solution of the EOM (B12) can be found in the power series expansion in 2−λ, with the leading term being simply f = sin(πt) π + O 1 log(2k/n 0 ) .
Comparing this with the asymptotic behavior of moments in 1D (56), we identify α = πJ/2, while matching n 0 would exceed the available precision of (56). behavior along the imaginary axis t = −iβ, its behavior along real axis is periodic and hence unphysical. Thus, it remains to be seen if for lattice models with local interactions b n is always asymptotically smoothly depend of n, or the behavior can be more complicated.