Detecting few-body quantum chaos: out-of-time ordered correlators at saturation

We study numerically and analytically the time dependence and saturation of out-of-time ordered correlators (OTOC) in chaotic few-body quantum-mechanical systems: quantum Henon-Heiles system (weakly chaotic), BMN matrix quantum mechanics (strongly chaotic) and Gaussian random matrix ensembles. The growth pattern of quantum-mechanical OTOC is complex and nonuniversal, with no clear exponential regime at relevant timescales in any of the examples studied (which is not in contradiction to the exponential growth found in the literature for many-body systems, i.e. fields). On the other hand, the plateau (saturated) value of OTOC reached at long times decreases with temperature in a simple and universal way: $\exp(\mathrm{const.}/T^2)$ for strong chaos (including random matrices) and $\exp(\mathrm{const.}/T)$ for weak chaos. For small matrices and sufficiently complex operators, there is also another, high-temperature regime where the saturated OTOC grows with temperature. Therefore, the plateau OTOC value is a meaningful indicator of few-body quantum chaos. We also discuss some general consequences of our findings for the AdS/CFT duality.


Introduction
Recent years have seen a renewed interest in classical and quantum chaos in the context of high-energy physics, black holes and AdS/CFT, thanks to the relation of chaos to quantum information theory and the information problems of black holes. Sharp and reasonably rigorous results such as the celebrated MSS chaos bound [1] and its subsequent refinements [2,3] establish a connection between chaos and the fundamental properties of gravity and black holes [4,5]. Maximal chaos, with the Lyapunov exponent λ = 2πT at temperature T , is reached for strongly coupled field theories in the large N limit, which have a classical gravity dual with a black hole. In [6] and other works it is explicitly shown how the Lyapunov exponent changes with finite N effects.
However, it has been pointed out many times, also in the pioneering MSS paper [1], that the multiple notions of quantum chaos in the literature mean different things. The out-of-time ordered correlation function (OTOC), given by the expectation value of the commutator of some operators A and B at times 0 and t: is a natural quantity in quantum field theories, i.e. many-body systems, and defines the quantum Lyapunov exponent λ as the exponent of the time growth of OTOC. However, the classical limit of this exponent does not necessarily have much to do with the classical Lyapunov exponent λ class , obtained by solving the variational equations [7][8][9]. The reason is the noncommutation of the three limits to be taken: the classical limit → 0, the long-time limit t → ∞, and the small initial variation limit δx(0) → 0. The crucial insight of [8,9] is that the mechanism of scrambling may be the chaotic dynamics, in which case λ class is related though still not identical to the OTOC exponent (quantum Lyapunov exponent) λ, or it may originate in local instability (hyperbolicity), in which case even regular systems may have a nonzero λ exponent and likewise chaotic systems may have λ which is completely unrelated to the classical counterpart.
This mismatch between the classical and quantum Lyapunov exponent is just the tip of the iceberg. The problem is twofold: not only what is the relation between the quantum (OTOC) exponent and classical chaos, but also what is the relation between the quantum Lyapunov exponent λ and other indicators of quantum chaos such as, first and foremost, level statistics. The bread and butter of quantum chaos is the famous Dyson threefold way leading to the Wigner surmise, the level repulsion statistics determined solely by the time reversal properties of the Hamiltonian [11], which follows from the random matrix approximation of chaotic Hamiltonian operators [12]. It is no secret for several years already that the black hole quasinormal mode spectra follow the random matrix statistics [13], and the OTOC of a Gaussian unitary ensemble (GUE) has been computed analytically in terms of Bessel functions in [14,15]; the outcome is close to the expected behavior of large-N field theories only at long timescales, longer than the scrambling time; at shorter timescales there are important differences. The authors of [14] have reached a deep conclusion in this respect: random matrices have no notion of locality as the correlation of any pair of eigenvalues is described by the same universal function. This is why the OTOC of a GUE system deviates from that of a local field theory at early times, when the perturbation in field theory has not had time to spread yet (i.e. when it is still localized). Therefore, the level repulsion does not imply the usual picture of the chaotic (exponential) OTOC behavior. However, we do not know yet how this correlates to the behavior of few-body or more precisely few-degrees-of-freedom quantum systems as opposed to the large-N field theories with a gravity dual. In few-body systems the notion of locality (and a classical gravity dual) does not exist anyway and the main problem found by [14] is irrelevant; at the same time, such systems are often very well described by random matrix statistics, i.e. Wigner-Dyson statistics [11]. In this paper we aim to understand the behavior of OTOC in such systems. Running a bit forward, we can say that the growth of OTOC is rather unremarkable: we find no universal trend, and little connection to level statistics. This confirms the results found for specific examples in [9,16]. The relation of OTOC, level statistics and the classical Lyapunov exponent was studied for few-body systems (quantum mechanics) in [9,[16][17][18][19] and the picture is inconclusive. One can have a nonzero growth exponent in integrable systems, 1 whereas fast scrambling with the exponent close to 2πT or at least growing linearly in T has not been found even in some clearly chaotic systems [9,10]. Arguments for many-body systems such as spin chains even suggest that quantum-chaotic systems with Gaussian spectral statistics generically never show fast scrambling [20], but no claims of such generality have been tested or formulated for few-body quantum chaos.
Various indicators of chaos relevant also for small systems, and their relation to OTOC and scrambling were studied by [21][22][23][24][25][26][27][28] among others. In particular, in [27] some important insights can be found: even in small systems devoid of the locality notion, OTOC can be interpreted as a measure of delocalization of a state in phase space, and the oscillatory component of the OTOC dynamics has to do with the power spectrum of the system. This last insight provokes a more general question: can we learn something from the quasistationary regime of OTOC, where no systematic growth is present but only oscillations? In this paper we provide a partial answer from a detailed study of this saturated (asymptotic, plateau) OTOC regime: the magnitude of the OTOC average at the plateau has a simple temperature dependence, and apparently can differentiate between weak chaos (dominantly Poissonian level distribution with some admixture of the Wigner-Dyson statistics) and strong chaos (clear Wigner-Dyson level repulsion). We will demonstrate this on three representative systems: the quantum Henon-Heiles Hamiltonian, whose classical limit has mixed (regular/chaotic) phase space and thus we expect on average weak chaos, a simplified BMN matrix model (at small N ) exhibiting strong chaos for most initial conditions, and Gaussian random matrices, the prototype of strong quantum chaos. The long-time limit of OTOC behaves in subtly different ways in each case.
Before we start, one caveat is in order (we will consider this issue in more detail later on): one might think that the saturated OTOC value is always trivially determined by the system size. We typically assume that the OTOC function C as defined in (1.1) behaves roughly as C(t) ∼ c/N 2 × exp(λt) with c of order unity, so when t ∼ t * ≡ log N 2 /λ the growth of C(t) stops and OTOC approximately reaches unity (when appropriately normalized). But the twist is precisely that c is system-specific and in general poorly known. The leading N 2 behavior indeed determines the OTOC values for N large, but when N and c are comparable within an order of magnitude the effects of fluctuations and finite N corrections are significant. This is at the root of our observations in this work.
The plan of the paper is as follows. In section 2 we recapitulate and generalize some results on computing OTOC in quantum mechanics, and show how OTOC sensitively depends on both the Hamiltonian and the operators A, B from the definition (1.1). In section 3 we apply the general formalism to random matrix ensembles and show that the OTOC growth is a complicated and nonuniversal function but that its asymptotic value behaves in a rather simple way. Section 4 discusses the behavior of OTOC for deterministic quantum-chaotic Hamiltonians. Section 5 sums up the conclusions. exponential growth even in absence of chaos.

OTOC in quantum-mechanical systems
Consider a four-point time-disordered correlation function for a quantum-mechanical system in 0 + 1 dimensions at temperature T = 1/β. Starting from the usual definition (1.1) as the squared module of the commutator of the two operators A and B, we can write it out as where the averaging is both thermal and quantum mechanical: . . . = tre −βH vac| . . . |vac . We can pick a basis of states and express the above defining expression in terms of matrix elements of the operators (this closely follows the derivation in [9,17]): where we have inserted the completeness relation 1 = m |m m|. For a single element c mn (t) one gets: where in the second line we have again inserted a completeness relation and in the third line we have used the fact that we work in the energy eigenbasis. The outcome is expressed in terms of the matrix elements a mn , b mn of the operators in the energy basis. In practice, it may or may not be possible to compute these analytically. Specifically, for A = x, B = p, we get the analogue of the classical Lyapunov exponent. From now on we call this the kinematic OTOC as it is directly related to the classical trajectory. Let us now see what general bounds can be put on (2.3) from the properties of quantum-mechanical Hamiltonians.

An upper bound on OTOC saturation
We begin with a very general and very formal result, which immediately makes it clear that in a generic quantum-mechanical system (integrable or nonintegrable) OTOC can be bounded from above by a quantity which solely depends on the energy spectrum of the  4) where N is the matrix size. In the second and third line we have used the inequality between the arithmetic and harmonic mean. Now we can bound the value of C(t): This means that C(t) is bounded at all times by an oscillatory function of time, whose frequencies are linear combinations of three eigenenergies Such a combination is generically always nonzero for a chaotic system except when the energies coincide, e.g. E m = E n = E k (according to the non-resonance condition). Therefore, since OTOC is typically a non-decreasing function of time, the behavior of C(t) for t large is roughly its maximum value and is likely close to the right-hand side in (2.5). This suggests that the OTOC dynamics after saturation likely consists of a very complex oscillatory pattern (with ∼ N 3 frequencies if the Hilbert space has dimension N ) superimposed on a plateau. The numerics will indeed confirm such behavior. Another estimate, which is time-independent and relevant for our main result -the magnitude of the saturation (plateau) OTOC value, is obtained from the triangle and mean inequalities: where we have used the obvious relations k C mnk = (A·B) nm and k D nmk = (B·A) nm = (A · B) * mn = (A · B) nm , assuming also the hermiticity of the operators. For some models (e.g. random matrices, Henon-Heiles), this sum can be estimated in a controlled way and provides an approximation for the plateau of OTOC. These estimates are obviously very simple and very weak (in the mathematical sense) but provide us with a framework into which we can insert specific A, B and H (the Hamiltonian with energies E n ) and perform back-of-the-envelope calculations which explain the numerical findings.

OTOC for random matrix ensembles
Random matrix theory [11,12] provides a highly detailed and rigorous (within its starting assumptions) stochastic effective description of the few-body quantum chaos, and allows an analytic calculation of OTOC along the lines of (2.3). Let us focus on Gaussian ortohogonal ensembles of size N ×N , appropriate when there is full time reversal invariance. It is known [11] that the joint distribution all the elements of all eigenvectors is obtained simply from the statistical independence of the eigenvectors from each other and of the elements in each eigenvector (and the orthogonality of the eigenvectors): where i = 1 . . . N is the component of the eigenvector and 1 ≤ n, m ≤ N count the eigenvectors themselves, i.e. the energy levels; so the n-th eigenvector |n is represented by the column vector ψ (n) with the elements (c n 1 , . . . c n N ). Special cases like the probability distribution for the p-tuple of elements of a single eigenvector are obtained from (3.1) by integrating out all the other elements [11]. We will also need the probability distribution of the energy levels {E} = E 1 , E 2 , . . . E N , the celebrated Wigner-Dyson distribution function [11]: where σ is the standard deviation, fixing the unit of energy, and b = 1, 2 or 4 for orthogonal, unitary and symplectic ensembles respectively. Most of our work is independent of the symmetry class, however our default class will be the Gaussian orthogonal ensemble (GOE) with b = 1 when not specified otherwise.

Estimate of the OTOC and its plateau
The idea is to use the results recapitulated in the previous section to find the ensemble expectation value of OTOC from the "master formulas" (2.2-2.3). Representing the eigenvectors and the operators as matrices in some (arbitrary) basis we can obviously write out and similarly for b nk . Inserting the above expression for the matrix elements into (2.3), multiplying c mn (t) by its complex conjugate taking into account the reality of the eigenvectors and relabelling the indices in the sums where convenient we find (denoting the average over the random matrix ensemble by C(t) ): where {c} determines the whole set of N 2 random elements c (n) j with j, n = 1 . . . N and likewise {E n } is the whole set of eigenenergies. All the sums run from 1 to N . The integral over the eigenvector elements {c} in (3.4) produces only an overall constant as these coefficients do not couple to the other quantities (in fact the integral d N 2 {c} is a textbook Jeans integral, but we do not need its value as it only produces an N -dependent, T -independent constant). The remaining integral, over the eigenenergies, is again a sum of products of Jeans-type integrals but with an additional linear term −βE in the exponent. Notice that the imaginary (sine) terms in (3.4) cancel out when the sum is performed; this is a consequence of the module squared in |c mn | 2 , i.e. of the reality of OTOC. Now we see that (3.4) becomes a sum where each term is a product of factors of the form where p i is some polynomial and s may be zero or unity, i.e. some terms have this factor and some do not. Every such term is a Jeans-type integral. The number of terms in P ({E}) equals the number of partitions of N (N − 1)b/2, and the sums over the coefficients {c} bring alltogether N 12 terms. When everything is said and done (for details see Appendix A), the final outcome, ignoring the multiplicative constant factors, reads: where 1 F 1 is the confluent hypergeometric function. The sum runs over all partitions of N (N −1)b/2, and the product has four terms as each factor |c mn | 2 has four matrix elements of A and B.

Kinematic OTOC
In order to move further we need to specify at least to some extent the operators A and B. We will consider (1) the kinematic OTOC, with A = x, B = p (2) generic sparse operators, with O(N ) nonzero elements in the matrices a mn and b mn , and (3) dense operators A and B, with O(N 2 ) nonzero elements, in particular the case when the operators A, B are themselves represented by Gaussian random matrices. Let us estimate OTOC for each case. For the kinematic OTOC, A ij = x i δ ij is diagonal and in the large-N limit B can be approximated as B ij ∼ δ ij /x i . The Kronecker deltas reduce the number of terms in the sums over {c} to N 4 , the number of partitions j α j = n can be approximated as p(n) ∼ exp(π 2n/3)/ √ n, and the general expression (3.4) becomes 2 where W 0,1,2 are polynomials in σβ of degree N (N − 1)b/2 ∼ N 2 b/2, Q 1 is an even polynomial (with only even powers) of the same degree, and Q 2 is an odd polynomial of the same degree. Each coefficient in the polynomials W 0,1,2 comes from ∼ N 2 terms (Appendix A), therefore the size of the coefficients scales approximately with N 2 . Eq. (3.6) is a very complicated oscillating function as many terms are involved. But if we are only interested in the average value of C(t) at long times, we may simply ignore the oscillations (which in the first approximation average out to some value of order unity) and write the estimate for the long-term, saturated or plateau OTOC value that we denote by C ∞ : We deliberately do not write lim t→∞ in the above definition as the limit in the strict sense does not exist because of the oscillatory functions, and in addition our derivation is obviously nothing but a crude estimate. A similarly rough estimate of the temperature dependence of C ∞ can be obtained in the following way. For sufficiently large σβ, roughly σβ/N 2 > 1, the polynomial W 0 is dominated by the highest-degree term and we have, from (3.7): where in the second step we have assumed β 1 so that the power-law prefactor β N 2 b/2 becomes negligible compared to the exponential. We deliberately emphasize that there are other terms in the expansion (. . .), including also a constant term (from the zerothorder term in W 0 ). This is important as it tells us that the scaling is in general of the form C ∞ ≈ const. + exp(σ 2 /4T 2 ), i.e. the temperature dependence is superimposed to a constant. This is also expected as the (appropriately normalized) saturated value C ∞ should always be of order unity, and the temperature dependence will only account for the relatively small differences between the plateau values of C(t), as we will see later in Figs. 1 and 2.
On the other hand, for sufficiently small σβ, the polynomial W 0 can be estimated as a geometric sum of monomials in −σβN 2 (remember the terms in W 0 have alternating signs): We have now reached an important point: the plateau OTOC falls off exponentially with 1/T 2 at low temperatures 3 and grows as a function of 1/T at high temperatures (we are not sure which function, as there are higher order terms in addition to the one written in (3.9), and there is no clearly dominant term like the exponential at large β), with the crossover temperature: 4 T c ∼ σN 2 . (3.10) If we consider a pair of arbitrary sparse operators A and B, the whole above reasoning remains in place, except that the products of matrix elements such as remain as arbitrary constants. Therefore we get the same qualitative behavior with two regimes and a crossover between them. The crossover temperature is very high for typical N 1 (otherwise the random matrix formalism makes little sense) and finite σ (again, σ → 0 makes little sense). In particular, in the N → ∞ limit the crossover temperature becomes infinite and the only regime is the exponential decay.

OTOC for dense and/or random operators
Now consider the case when the matrix elements in (3.5) are generically all nonzero (and for now nonrandom, i.e. we fix the operators and do not average over them). The large-t limit yields the expression where q 0,1 and w 0 are polynomials of degree N (N −1)b/2 ∼ N 2 b/2, and w 1 is the polynomial of the same total degree of two variables, σβ + ıt and σβ − ıt. The coefficients of w 0,1 are proportional to products of matrix elements but now a typical coefficient of the polynomials w 0,1 behaves as N 2 (|A||B|) 8 . Therefore, the scaling in the low-temperature regime remains the same as (3.8): But the high-temperature regime yields therefore the crossover now happens at and therefore may be lower than the very high value (3.10), depending on the norm of the operators A and B. Finally, if the operators A and B are both random Hermitian matrices (for concreteness, from the Gaussian unitary ensemble with the distribution function Π and the standard deviation ξ), we need to average also over the distribution functions for A and B and work with the double average C(t) : This estimate is very crude, based simply on the fact that the distribution functions Π have ∼ N 2 /2 pairs of the form (a i − a j ) 2 . The important point is that the scaling from (3.13) that behaves essentially as ∼ ξ 16 now becomes ∼ ξ N 2 , therefore the crossover temperature is significantly reduced compared to (3.14) and behaves as T c ∼ σN 2 λ N 2 . So for random operators the crossover may happen at temperatures that are not very high and thus can be clearly visible in the numerics and experiment.

Numerical checks
Now we demonstrate numerically that the crude estimates from the previous subsection indeed make sense and describe the characteristic behavior of OTOC. Our chief goal is to understand the behavior of C ∞ , however it is instructive to start from the time dependence of the kinematic OTOC (Fig. 1). We find the expected pattern of early growth followed by a plateau, however the growth is closer to a power law than to an exponential; this follows from the polynomial terms in (3.6), although the power law is not perfect either, as we see in the panels (A, C). This is in line with the prediction of [14], where the authors find for a slightly different ensemble of random matrices (J 1 is the Bessel function of the first kind). This function is also neither an exponential nor a power law but at early times it is best approximated by a power law at leading order (at long times it falls off exponentially but the saturation is reached already prior to that epoch). In Fig. 2 we focus on the plateau behavior. It has the form of a constant function with superimposed aperiodic oscillations, and the differences of the plateau values are the subject of our theoretical predictions.
These are relatively small and become important only when N is finite and not very large.
In Fig. 2 we plot again the time dependence of the kinematic OTOC but now we focus on long timescales, to confirm that the plateau is indeed stable, and to show the very complex oscillation pattern superimposed on the plateau. At small inverse temperatures there is some deviation from the linear scaling law but this we also expect. Looking now at the OTOC for a pair of random Hermitian operators in Fig. 4, we detect also the other regime at small enough temperatureslog C ∞ decays with the inverse temperature. This regime is likely present also in Fig. 3, but only at extremely high temperatures (which we have not computed in that figure).

OTOC for weakly and strongly chaotic Hamiltonians
For quasi-integrable few-degrees-of-freedom Hamiltonians one would naively expect that OTOC closely resembles the Lyapunov exponent, at least for high quantum numbers, approaching the classical regime. As we have already commented in the Introduction, it is known that this is not true in general [7,9,10,17,18] and that both chaotic systems with zero quantum Lyapunov exponent and regular systems with a nonzero exponent exist. We will now try to find some common denominator of OTOC dynamics in (deterministic) quantum-mechanical systems. It will quickly become clear from our general analysis of the master formula (2.3) that the function C(t) is as complicated as for random matrices (indeed, even more so). But we will again construct an upper bound for the saturated OTOC value and arrive at a rough scaling estimate.

Hamiltonian
a well-known paradigm for classical chaos with applications in galactic dynamics and condensed matter. For = 0 it obviously reduces to a 2D linear harmonic oscillator and becomes integrable. As we know, nonintegrability does not always imply chaos; indeed, this is a typical system with mixed phase space, with both chaotic and regular orbits. Chaotic orbits proliferate only when the perturbation is larger than some finite c ; they are almost absent at low energies, numerous at intermediate energies and again absent at very high energies when the potential energy is negligible compared to the kinetic energy [29,30]. For such a quasi-integrable system our idea is to apply elementary perturbation theory in the occupation number basis to estimate OTOC starting from (3.4). We will present the perturbation theory in a fully general way, for an arbitrary Hamiltonian H 0 + V , and some of the conclusions will also turn out to be quite general. Only at the end we will show the quantitative results for the Henon-Heiles system (4.1).
Let us write the perturbative expression for OTOC. Replacing the initial basis states |n with the first-order 5 perturbatively corrected states |n + δn and introducing likewise the perturbative corrections δa mn , δb mn for the matrix elements of A and B, the equation Now we insert this result into (2.2) and apply the Cauchy-Schwarz-Bunyakovski inequality: whereρ ≡ diag(exp(−βE n )) is the non-normalized density matrix. In the above derivation, we have also used the definition of trace and the definition of thermal expectation values A ≡ Tr (ρ · A). The norm of a matrix is defined as usual by |A| 2 ≡ TrA † A. This estimate manifestly replaces C(t) by its asymptotic (maximum) value, as we have replaced the terms containing the time-dependent phase factors by their time-independent norms. In order to move further, notice first that |ρ| 2 = n exp(−2βE n ) = Z(2β) for a canonical ensemble with the diagonal density matrix that we consider here. This means, from (4.3), that the temperature dependence is encapsulated in the ratio Z(2β)/Z(β). The prefactor will again differ between dense and sparse A, B and V . For sparse matrices, we can write |A| 2 ∼ N a 2 whereas for dense matrices we have |A| 2 ∼ N 2 a 2 , assuming that all matrix elements have some characteristic magnitude a. Obviously, if this is not true the outcome will be more complicated, but it seems this does not influence the temperature dependence. For concreteness we assume sparse A and B. For sparse V with nonzero elements of order v concentrated near the diagonal (this is true for the Henon-Heiles Hamiltonian and in general for perturbations expressed as low-degree polynomials in coordinates and momenta), we can estimate |δ| 2 ∼ N v 2 /ω 2 . Here we assume an approximately equidistant spectrum of H 0 with frequency (neighboring level spacing) ω. This yields: For a dense perturbation V , the only factor that changes is |δ|. Assuming again the utterly simple situation where all matrix elements of V are roughly equal v, we have δ mn ∼ v/E mn ∼ v/(ω(m − n)), which yields a series that can be summed analytically. However, we will not pursue this further as the temperature dependence is universal in all cases, given by the simple ratio of the partition functions: (4.5) The first simplification holds when H 0 is a 1D harmonic oscillator, and the second one when N → ∞. But the basic result (the ratio of partition functions) always holds. We are in fact more interested in the 2D harmonic oscillator, which is the integrable part of the Henon-Heiles Hamiltonian. For that case, we plot the sum (for finite N ) in Fig. 5.
Of course, the analysis of the function R(β) is trivial -we plot it in the figure merely to emphasize the qualitative agreement with the numerics.  Figure 5. The thermal dependence factor of log C ∞ for weak perturbative chaotic systems, given by Z(2β)/Z(β) with β = 1/T the inverse temperature, here given for a 2D linear harmonic oscillator with the frequencies ω x = ω y = 0.1, with N = 20 levels (blue) and with N = 150 levels (red). We also plot the sum over N = ∞ levels (black). In (A) we zoom in at high β/low temperatures, and in (B) we focus on smaller β/larger T . The N = 150 plot is aready quite close to the monotonic N = ∞ dependence but at high temperatures there is always a region decaying with β, before the approximately linear log R(1/T ) dependence sets in, just like in the numerical results. At very low temperatures the ratio saturates, as we see in the panel (A). This ensures that our estimate for the saturated OTOC has a finite limite at zero temperature. The overall scale is arbitrary as the R factor is always multiplied by various other factors.
As a final remark, what we have found is the correction of the OTOC plateau δC ∞ . There is still the unperturbed value of C ∞ for the integrable Hamiltonian H 0 . We know that this can be nonzero and even quite large because of local instability [7,8,10]. We are mainly interested in the opposite situation, when the scrambling chiefly comes from chaos so that OTOC does not grow when H = H 0 . In this case C (1) ≈ δC ∞ and the temperature dependence is primarily determined by (4.4). In the next subsection we will look both at the Henon-Heiles system where this holds, and a perturbed inverse chaotic oscillator where H 0 is unstable.

Weak chaos: examples and numerics
As our main example we can now study the Henon-Heiles system of Eq. (4.1). Starting from the nonperturbed Hamiltonian (2D harmonic oscillator), we will express the nonzero elements of c mn (t) exactly, i.e. we will not use the estimates (4.3) as the perturbation is quite simple and amenable to analytic treatment. Denoting a basis state by the quantum number n = (n x , n y ), we can write the amplitudes c nxnyn x n y as products of amplitudes of the two decoupled subsystems c nxnyn x n y = C nxn x C nyn y , with and all other elements are zero; for the y quantum numbers the coefficients are the same with (n x , ω x ) → (n y , ω y ). For nonzero , the off-diagonal matrix elements can be represented exactly as c nxn x nyn y (t) = δ |nx−n x |−2 δ |ny−n y |−1 m x (m x − 1) m y + 1, m x,y ≡ min(n x,y , n x,y ).
(4.7) The state vectors are now easily calculated in textbook perturbation theory. We have compared the analytic calculation to the numerics and find that they agree within a relative error ≤ 0.04; therefore, one may safely use (4.6-4.7) in order to speed up the computations and avoid numerical diagonalization of large matrices.
The magnitude of the plateau value of C(t), computed by long-time averaging similarly to the random matrix calculations in Section 3, are given in Fig. 6. At large T values, C ∞ decays with 1/T , at intermediate values it shows an exponential growth with 1/T just like R(1/T ) in Fig. 5, and as the temperature goes to zero it reaches a finite value. In Fig. 7 we consider a system with much reduced state space, with N = 25. We expect that for small N the existence of two regimes is more clearly visible, and that the crossover temperature is higher. This is indeed what happens, although the exact form of the function C ∞ (1/T ) is not very well described by the analytical result. As we have made many crude approximations, this is not surprising: our analytical result still explains the existence of two regimes and the crossover between them. 6 One unexpected finding is that the hightemperature regime is apparently universal for all perturbation strengths and scales as C ∞ ∝ exp(−4π/T ). This is probably specific for the Henon-Heiles system; we do not understand it at present.
It is instructive to look at the energy level statistics of the Henon-Heiles system for the same parameters that we have used for the OTOC calculation, in order to understand the relation of OTOC to chaos. In Fig. 8  1 T log C Figure 6. (A) The saturated kinematic OTOC value C ∞ for a range of inverse temperatures β = 10 −4 , 10 −3 , 5 × 10 −3 , 0.01, 0.02, 0.05 and a range of perturbation strengths = 1, 2, 5, 10, 15, 20 (black, blue, green, magenta, red, orange). Here we see the scaling C ∞ ∝ exp(c/T ), with c growing with . In the (B) panel we zoom in the high-temperature region, to show that for ≤ 5 there is also the other regime where C ∞ grows with T . Since the number of points in this interval is small it is not easy to judge the form of T -dependence. In (C) we focus on the opposite regime, at very low temperature, showing that C ∞ saturates as T → 0. This is again in accordance with the β → ∞ limit of Z(2β)/Z(β). The system is truncated to N = 144 levels.  Figure 7. The saturated kinematic OTOC value C ∞ (T ) for the truncated Henon-Heiles model with N = 25 levels. The perturbation strength is = 0.1, 0.3, 0.5, 0.7, 0.9, 1.5, 3.0 (black, blue, green, magenta, red, orange, cyan). Already in (A) we see that for ≤ 0.5 there is a finite crossover temperature T c so that C ∞ grows wth temperature for T > T c . Since T c goes down when the Hilbert space is reduced, we can observe the high-temperature regime very clearly and see that it collapses to a universal law C ∞ ∼ exp(−4π/T ). This is seen in the panel (B) where we zoom in at the interesting region.
Even for large , the regular (Poisson) component is dominant over the chaotic (Wigner-Dyson) component. In other words, the classically mixed phase space, with the increasing chaotic component, is almost completely regular in the quantum regime; quantum chaos is "weaker", as is often found in the literature [34]. For us, the fact that the system is outside the Wigner-Dyson regime means that indeed the behavior of C ∞ is a good litmus test of quantum dynamics, behaving (at low temperatures) as exp(1/T 2 ) or exp(1/T ) for strong or weak chaos respectively. 7 Indeed, we would not expect that a system which is well described by perturbation theory shows strong level repulsion.
Finally, it is instructive to look at the inverse Henon-Heiles system, with (ω 2 x , ω 2 y ) → (−ω 2 x , −ω 2 y ), so that H 0 is the inverse harmonic oscillator. As already found in the literature, scrambling is significant already at = 0, and this contribution dominates even at high , at all temperatures. In other words, neither the perturbation nor the temperature have a significant influence over C ∞ . This is fully in accordance with the result (4.3) and the morale is that OTOC directly describes scrambling, and chaos only indirectly, through the scrambling, if the scrambling originates mainly from chaos; if not, OTOC is largely insensitive to chaos. Therefore, the temperature dependence of the OTOC value, derived from the assumptions about the dynamics (perturbative chaos or strong, random-matrix chaos) cannot be seen when there is a strong non-chaotic component to scrambling. 1 T log C Figure 9. Temperature dependence of the asymptotic OTOC C ∞ for the inverse Henon-Heiles Hamiltonian, with ω 2 x = 4ω 2 y = −1, and perturbation strength = 0, 0.1, 0.5, 0.9, 1.5 (blue, green, magenta, red, orange). The curves are more or less flat and without a clear trend, as the scrambling is rooted in local instability, not chaos.

Strong chaos: numerics and the return to random matrices
As a final stroke, we will now examine a strongly chaotic system which is also relevant for black hole scrambling and related problems in high energy physics. This is the bosonic sector of the D0 brane matrix model known as the BMN (Berenstein-Maldacena-Nastase) model [31], obtained as a deformation of the BFSS (Banks-Fischler-Shenker-Susskind) model [32] by a mass term and a Chern-Simons term. This model has been studied in detail in the context of non-perturbative string and M theory. It is known to describe the dynamics of M theory on pp-waves and is also related to the type IIA string theory at high energies; for details one can look at the original papers or the review [33]. Following [34][35][36][37], we focus solely on the bosonic sector which is enough to have strongly chaotic dynamics with equations of motion that are not too complicated. The quantum-mechanical Hamiltonian of the BMN bosonic sector reads: 8 where Π i are the canonical momenta, X i the canonical variables, ε abc is the Levi-Civitta tensor, and ν 2 > 0 is the mass deformation which also determines the size of the Chern-Simons deformation (the last term in (4.8)). Following [36], we will study the "mini-BMN" model with X α = 0, so we effectively only have three degrees of freedom. The matrices X 1,2,3 and P 1,2,3 are N × N matrices. For this example we have to abandon the master formulas for OTOC (2.2-2.3) as it is very difficult to find the quantities c mn -for this we would have to perform exact diagonalization of the Hamiltonian (4.8). Instead, we follow [37] and write a truncated system of equations directly for the expectation values X a and P a and the two-point correlators X a X b , Π a X b and Π a Π b , where the expectation value is obtained through the trace over the density matrix: X a ≡ Tr(ρX a ).
The equations read (for their derivation see [37]): where the last line contains the leading quantum corrections: all possible terms with a single contraction of the classical equation of motion, and the summation over repeated indices is understood. The equations of motion for the two-point correlators are obtained again by writing the classical equations of motion for the bilinears Π a Π b , Π a X b and X a X b and taking all possible single contractions in each term. This yields: As explained in [37], this truncated system is obtained by assuming a Gaussian approximation for the wavefunctions. Therefore, we solve the truncated quantum dynamics of the mini-BMN model -essentially a toy model, but it will serve our purpose. Now that we have set the stage, we can express the kinematic OTOC as C(t) = X a Π b − Π a X b and study its dynamics. The outcome is given in Fig. 10. We are essentially back to the random matrix regime of Section 3 -there is a clear scaling C ∞ ∼ exp(1/T 2 ) (we do not see the other regime, but again it may well be there for sufficiently high temperatures), and the level distribution is a near-perfect fit to the Wigner-Dyson curve. Therefore, if a Hamiltonian is strongly chaotic, then both the level distribution and the OTOC plateau are well described by the random matrix theory.

Discussion and conclusions
In this paper we have formulated a somewhat unexpected indicator of quantum chaos, useful mainly in few-body (few-degrees-of-freedom) systems. While OTOC has become the quintessential object in the studies of quantum chaos and information transport, characterized mainly by its growth rate -the (quantum) Lyapunov exponent, in our examples its growth pattern tends to be quite nonuniversal and "noisy" (in the sense that it depends sensitively on the system at hand and the operators we look at). Our analytic treatment of OTOC dynamics is quite sketchy, however both analytical and numerical results strongly suggest there is no clear exponential growth. At first glance, one might think that this finding is completely at odds with the established wisdom, however this is not true. In the literature, exponential growth is mainly characteristic for systems with a classical gravity dual (and reaches its maximum when the dual contains a thermal black hole horizon).
There are abundant examples of quantum chaotic systems which do not have an exponentially growing OTOC (we especially like [20] but there are many other published examples). The exponential growth follows, in the AdS/CFT picture, from the shock wave dynamics in a classical gravity background, and need not exist when the background is not classical or when the gravity dual does not exist at all. This is precisely what happens here: the Henon-Heiles Hamiltonian is certainly nothing like a strongly coupled large N field theory, while the truncated mini-BMN model comes closer (it is actually related to discretized Yang-Mills) but we tackle it at finite N and thus away from the fast scrambling dual. For random matrices, our findings for C(t) are in line with the rigorous results of [14]. As pointed out in that work, the crucial difference between random matrices and strongly coupled field theories is that the former have no notion of locality neither in time nor in space. In our small systems, the spatial locality is irrelevant anyway but if the system is not sufficiently chaotic there will still be long-term temporal correlations in dynamics (this indeed gives rise to different scaling regimes for strong and weak chaos).
On the other hand, what we have found is that the long-time OTOC behavior, when it becomes essentially stationary, with a complex oscillation pattern, is surprisingly regularbehaving as exp(1/T 2 ) and exp(1/T ) respectively in strong and weak chaos. This indicator seems to have a stronger connection to quantum chaos in the sense of level statistics than the Lyapunov exponent; in all examples we have studied the exp(1/T 2 ) regime and the Wigner-Dyson level distribution go hand in hand. At very high temperatures we detect also a different regime, when the OTOC plateau grows with temperature. This regime seems less universal, and we do not understand it very well. One might think that the plateau value should not carry any useful information; it is often laconically stated that OTOC reaches saturation when the initial perturbation has spread all over the system and that this saturation value is unity when OTOC is appropriately normalized. This is roughly true, however "spreading all over the system" is not a rigorous notion -depending on the system and the operators A, B in OTOC, the perturbation may never spread completely due to symmetry constraints, specific initial conditions, quasi-integrals of motion etc. Such factors are particularly important in finite systems (quantum mechanics as opposed to quantum field theory) that we study. Looking at the figures, one sees that differences in the asymptotic OTOC value C ∞ tend to be small, and C ∞ tends to be about the same to an order of magnitude in all cases. We conjecture that such differences would dwindle to zero in the field limit.
A simple intuitive explanation for the falloff of asymptotic OTOC with temperature is the following: we expect that higher temperatures lead to faster information spreading and quicker equilibration. Therefore, it is logical that the plateau value will be lower, so that the system needs less time to reach it, i.e. it needs less time to equilibrate.
We note in passing that we have confirmed that scrambling can originate from at least two distinct mechanisms: local instability and chaos, so in the former case the relation of OTOC to chaos is largeley lost. This is a known fact in many examples already [7,8,10,17,27] and we emphasize it here merely as a reminder to the reader that the OTOC-chaos connection is really a relation of three elements: OTOC-scrambling-chaos, and if the second link is missing no attempt should be made to understand chaotic dynamics from OTOC.
We conclude with some speculations. The OTOC plateau value, as we found, is a rather universal function of temperature, and it is essentially a finite-size fluctuation of the correlation function, when the system is small enough that the relative size of fluctuations does not go to zero. We may then look for universality and the connections to chaotic dynamics in other similar quantities, e.g. the average fluctuation of the expectation value of some operator during thermalization. Such a quantity remains nonzero also in AdS/CFT at large N , and may relate our results to the more familiar fast scrambling, strongly correlated holographic systems.

Acknowledgments
This work has made use of the excellent Sci-Hub service. Work at the Institute of Physics is funded by the Ministry of Education, Science and Technological Development and by the Science Fund of the Republic of Serbia, under the Key2SM project (PROMIS program, Grant No. 6066160).

A Detailed structure and calculation of OTOC for Gaussian ortohogonal ensembles
In this Appendix we consider the calculation of OTOC for random matrix systems in some more detail, and describe the detailed structure of the correlation function C(t). Let us first denote, for the sake of brevity: Denote also the products of matrix elements of the operators A, B entering the expression (3.4) by χ 1 , χ 2 , χ 3 , χ 4 . Now the expression for C(t) can be written as: (A.1) As we have noticed in the main text, the integral over C yields just a numerical constant. Let us therefore evaluate the energy integral I 1 = d N {E}P({E})e −βEn χ 1 e i(E k −E k )t . We have: The absolute values of the differences can be written out in the obvious way: where all α i,j are some (positive) integer exponents and π(i) is the appropriate sign factor −1 or 1. Therefore, I 1 can be reorganized as: (A.4) Note that the part n<m |E n − E m |, is not essential for the general behavior, since the singular integral E which yields the closed-form expression for the temperature dependence of I 1 : I 1 ∼ δ k,k e β 2 /4 + 1 − δ k,k e β 2 /4 e −t 2 /2 ., (A.7) where δ k,k is the Kronecker delta, reminding us that the main contribution comes from the terms with E k = E k which generically means k = k . It is clear that a similar calculation holds for the other parts of C(t) . This produces the temperature scaling found in the main text for random matrices, of the form C(t) ∼ e 1/4T 2 . But the time dependence is more complicated. In order to see this, we look at the structure of the polynomial factors in I 1 in some more detail. We see immediately that C(t) will also have dependence on t 2n , β n . Start from where E i = u − it/2. Let us look at two cases: α i even and α i odd. For any α i the polynomial will have the form: Assume first that α i is even. This means that j and α i − j are of same parity. For even j the Gaussian integral evaluates to some constant, but we will also have the prefactor of (it/2) α i −j , for all even j ≤ α i . The odd powers (j odd) will disappear because of the symmetric domain of integration. For α i odd, j and α i − j will be of different parity so again, only even j give a nonzero integral. In conclusion, the integral (A.6) with polynomial prefactors included will have the form: where Q(t 2n ) is a real polynomial depending on even powers of t, and 2n ≤ α i . Alternatively, for α i odd, we get: where R(t 2n+1 ) is a real polynomial depending on odd powers of t, and 2n + 1 ≤ α i . Analogous logic holds for the β dependence. Now we look back at I 1 : When we write out the products of energies, we have the following types of monomials in the resulting polynomial: A.1 The large matrix limit In the limit N −→ ∞ we can say more on the structure of OTOC. We can first schematically rewrite (A.15) together with any exponentially suppressed corrections as Here we have first absorbed all time and β dependence of (A.15) into the functions Q and W respectively, and then we have included the exponentially suppressed correction coming from the k = k terms in the integrals I 1 and K n . By L 1 , L 2 we denote the constant (timeand temperature-independent) factors. In general one can write L 1 as We can easily estimate the second sum. Namely, for N −→ ∞. Exactly the same logic goes for L 2 .
In the large matrix limit it is possible to show explicitly what we know has to happen: OTOC reaches a plateau. Looking at (A.16), the condition to reach the plateau for times longer than some scale t 0 is e −t 2 2 Q(t 2n )(L 1 + L 2 e −t 2 /2 ) = const. t > t 0 .
(A. 22) It is more convenient to look at the forms given in (A.15). First let us look at the condition Q(t 2n ) = const. × e t 2 /2 . The exponential term can be represented as a series; equating it with Q(t 2n ) we get j α j t 2j = const. × j t 2j 2 j j! , (A.23) thus, we need α j ∼ 1 2 j j! , which we know is the case from (A.9). For the second term the situation is similar: j β j t 2j = const. × j t 2j j! , (A. 24) so we need to have β j ∼ 1 j! ; this is true by cos(βt/2) = Q(t 2n )W (β 2n ) and sin(βt/2) = Q(t 2n+1 )W (β 2n+1 ), since the terms in the Taylor expansions of the left-hand sides behave as ∼ 1 j! . We can also look at the opposite limit in which t −→ 0. Let us rearrange (A.16): Now, simply using the definition of Q and expanding into a series we get: After some algebra we get: We see now that OTOC behaves in a very simple way for early times; this expansion is also consistent with the result (3.16) of [14].