Quantum Lyapunov Spectrum

We introduce a simple quantum generalization of the spectrum of classical Lyapunov exponents. We apply it to the SYK and XXZ models, and study the Lyapunov growth and entropy production. Our numerical results suggest that a black hole is not just the fastest scrambler, but also the fastest entropy generator. We also study the statistical features of the quantum Lyapunov spectrum and find universal random matrix behavior, which resembles the recently-found universality in classical chaos. The random matrix behavior is lost when the system is deformed away from chaos, towards integrability or a many-body localized phase. We propose that quantum systems holographically dual to gravity satisfy this universality in a strong form. We further argue that the quantum Lyapunov spectrum contains important additional information beyond the largest Lyapunov exponent and hence provides us with a better characterization of chaos in quantum systems.


Introduction
Many-body quantum chaos is of fundamental interest in a variety of fields of physics, including condensed matter, quantum information, and quantum gravity. Considerable recent progress has come from the realization that it is possible, in some cases, to define a kind of quantum butterfly effect and a corresponding quantum Lyapunov exponent via so-called out-of-time-order correlators (OTOCs) [1,2]. It was also discovered that the exponent so defined obeys a universal bound, the Maldacena-Shenker-Stanford (MSS) bound, λ N ≤ 2πT [3], and that the bound is saturated by strongly coupled quantum systems holographically dual to Einstein gravity [4,5].
Given these rapid developments, it is natural and desirable to attempt to extend the notion of a classical Lyapunov spectrum to general quantum systems away from the classical limit. In this paper, we give one definition of a quantum Lyapunov spectrum which makes sense for any many-body quantum system. There are several motivations for this study. One question is do black holes become more or less chaotic as they grow in size? To answer this basic question, we need to define the strength of chaos precisely. As explained in Sec. 2, the KS entropy is a better indicator than the largest exponent alone. It is likely that, as a black hole grows, the KS entropy increases, while the largest Lyapunov exponent decreases [54]. Hence looking at the entire Lyapunov spectrum gives a different picture of chaos in black holes than the largest Lyapunov exponent alone. Another motivation is universality in the classical Lyapunov spectrum [55]: in some chaotic systems, the classical spectrum converges to a random matrix theory (RMT) spectrum, with the time scale for the onset of universality seemingly related to the strength of chaos. In particular, a matrix model of black holes [56,57,58] shows the universality from t = 0, indicating a possible signature of gravity in the Lyapunov spectrum.
In this paper, we introduce a generalization of the classical Lyapunov spectrum to quantum many-body systems, including systems far from any classical limit. Our definition has a natural physical interpretation and, when the classical limit can be taken, it reduces to the usual classical Lyapunov spectrum. Furthermore, the largest Lyapunov exponent agrees with the usual quantum Lyapunov exponent in the literature as obtained from OTOCs at sufficiently late times.
To elucidate the physics of our definition, we study two systems, the non-local Sachdev-Ye-Kitaev (SYK) model [5,59,60] and the local XXZ model with random magnetic field (see e.g. [61]). In the case of SYK, disordered couplings are part of the basic definition. In the case of XXZ, the model is considered in the isotropic (Heisenberg) limit with an additional disordered magnetic field. We analyze both by performing systematic exact diagonalization studies of the quantum Lyapunov spectrum including disorder averaging.
Our main results can be summarized as follows.
• We observe a period of approximately exponential growth, meaning an approximately time-independent Lyapunov spectrum, for the nonlocal SYK model. For the local XXZ model we observe a brief period of approximately exponential growth followed by a period of power law growth.
• We numerically demonstrate that a naive generalization of the classical KS entropy -just the sum of positive quantum Lyapunov exponents -is close to the production rate of the entanglement entropy, with a proper normalization needed for the connection to classical coarse-grained entropy.
• We also suggest that this 'quantum KS entropy' is maximized when the quantum system has dual gravity description on the black hole background. In other words, we suggest that black hole is not just the fastest scrambler but also the fastest entropy generator.
• For both SYK and XXZ models in the appropriate regime we observe random matrix theory (RMT) behavior for the quantum Lyapunov spectrum. This universality fails in expected ways, including for integrable and localized systems.
The outline of the remainder of the paper is as follows. Section 2 recalls the basic definition of the classical Lyapunov spectrum. Section 3 introduces the SYK and XXZ models. Section 4 defines the notion of a quantum Lyapunov exponent. Section 5 studies the growth characteristics of the quantum spectrum with time. Section 6 studies the distribution of the quantum Lyapunov spectrum, establishing a link to random matrix statistics. Section 7 contains concluding remarks and outlook.

Lyapunov Spectrum in Classical Chaos
First let us see how the classical Lyapunov exponents are defined. For simplicity we consider only Hamiltonian systems. Suppose the phase space is described by coordinate variables x i (i = 1, 2, · · · , N ) and conjugate momenta p i (i = 1, 2, · · · , N ). We use z i (i = 1, 2, · · · , 2N ) to denote x i and p i together. The sensitivity to the initial condition at t = 0 is captured by The finite-time Lyapunov exponents λ i (t) are defined by where s i (t) are singular values of M ij (t). These exponents converge to constant values at sufficiently late time. Note that they are calculated for each initial condition, which is analogous to the microstate in the quantum theory. Often the average over initial conditions is taken. Note also that we can calculate the exponents from the eigenvalues of which are (s i (t)) 2 .
Often the limit t → ∞ is considered. In this paper, we will be interested in the finite-t behavior.

Kolmogorov-Sinai Entropy
Suppose our knowledge about the initial condition is limited and we only know that it is in a small region in the phase space, say the blue disk in Fig. 1. As time passes by, this region is stretched along some directions and compressed along the others. If we introduce a grid like in Fig. 1 and count the number of cells needed for covering the region, more and more cells are needed at later time; the number scales as This exponential growth characterizes the loss of our knowledge about a given initial condition. The coarse-grained entropy, which is the log of the uncertainty, increases as h KS t. The Kolmogorov-Sinai (KS) entropy h KS is the growth rate of the coarse-grained entropy [62] (although for exceptions see Ref. [63]). 1 The quantum counterpart of the coarse-grained entropy is the entanglement entropy. Indeed, the growth of the entanglement entropy of a subregion of a system is an essential part of quantum thermalization. It would be desirable if the KS entropy h KS could be identified with the growth rate of the entanglement entropy as well, at least near the classical limit.

Application to Black Hole
In the sense discussed above, the spreading of information in phase space is better characterized by the KS entropy than the largest Lyapunov exponent alone. This is particularly so for a black hole. In order to understand why, let us consider the matrix model of black hole made of D0-branes and strings [56,57,58], which is the original setup used by Sekino and Susskind to argue for fast scrambling [7]. (Essentially the same argument applies to black hole described by other gauge theories as well.) This model contains nine N × N Hermitian matrices as bosonic degrees of freedom. Diagonal and off-diagonal elements are interpreted as the locations of D0-branes in nine-dimensional space and open strings connecting D0-branes. As shown in Fig. 2, a big black hole consisting of all N D0-branes and Figure 1: Suppose we know that the initial condition is contained in the blue disk in the upper left corner. As this region evolves with time, although the volume is conserved due to the Liouville theorem, the shape changes nontrivially, in particular, it is stretched exponentially in several directions. The rate is governed by the positive Lyapunov exponents. many open string excitations is described by fully excited N × N matrices, while block diagonal configurations describe several black holes and the diagonal matrices describe a gas of D0-branes without strings.
Suppose two black holes consisting of N 2 × N 2 blocks (the middle of Fig. 2) come close and merge to a single black hole (the left of Fig. 2). During this process, open strings stretched between two black holes (off-diagonal blocks) become shorter, and hence lighter, and eventually get excited. In this way the number of dynamical degrees of freedom are doubled, and hence temperature, which is roughly equivalent to the energy per degree of freedom, goes down [54]. For example, in the highly stringy region where the matrix model admits a classical description, the energy is E = 6N 2 T fin = 2 × 6 N 2 2 T init , and hence T fin = 1 2 T init . This is the gauge theory description of a key property of a black hole discovered in Hawking's seminal paper [64]: a black hole in flat spacetime cools down as it grows in size. As the black hole cools down, the largest Lyapunov exponent becomes smaller; in the highly stringy region, it scales as T 1/4 [65]. Hence the largest Lyapunov exponent becomes smaller. However the number of positive Lyapunov exponents increases because of the dynamical increase of the degrees of freedom. Hence there are two competing contributions -the decrease of the temperature pushes down the KS entropy, while the increase of the degrees of freedom pushes it up. A careful evaluation of these effects shows that the KS entropy actually increases as black hole grows [54]. Hence a bigger black hole is a faster entropy generator.

Universality in the Lyapunov Spectrum
The Lyapunov spectrum of the classical limit of the matrix model of black holes has been studied in Ref. [65] and Ref. [55]. In particular, the level correlations of the Lyapunov spectrum have been studied. Ref. [55] studied other classically chaotic systems as well, and found a universal random-matrix description: when the number of degrees of freedom is sufficiently large, the level correlation of the finite-time Lyapunov exponents converges to the one determined by a Gaussian random matrix theory. Furthermore, the matrix model of a black hole shows this universality already at t = 0.
If this universality can be generalized to quantum chaos, it should serve as a new characterization of quantum chaos. Note that this universality is different from (though it might be related to) the usual Wigner-Dyson universality of the energy spectrum. First of all we are looking at different things -Lyapunov exponents and energy levels -and furthermore the onsets of the universal distribution are observed at different time scales: while the latter sets in at rather late time (see [66] for a detailed discussion), the former can set in at earlier time, as it is there already at t = 0 for the matrix model of black holes. As we will see later in this paper, the SYK model shows this Lyapunov universality at t = 0 as well.

Technical Remarks
A few technical remarks are in order here. Firstly, there is an ambiguity in the definition of the Lyapunov exponents at finite time associated with a choice of variables, although there exists a unique t → ∞ limit. If we choose another basis z = z (z), M ij (t) = δz i (t) where J ij = δz i δz j is the Jacobian matrix. Still, this definition captures the Lyapunov growth. Furthermore, at sufficiently late time, the Lyapunov exponents converge to the same values. The choice of the basis may affect the details at early time, and hence we will choose a natural one which makes the physical interpretation more transparent.
Secondly, for the Hamiltonian system, M is symplectic, and the Lyapunov exponents form pairs of +λ and −λ. It provides us with an easy way to estimate the numerical error, namely this pair structure is gone when the error accumulates too much. This property does not necessarily hold for the quantum Lyapunov exponents introduced in Sec. 4.
Thirdly, the finite-time Lyapunov exponents depend on the initial condition z i (0). In the thermodynamic limit, as long as a sufficiently generic initial condition is chosen, the same global structure (overall distribution) and microscopic spectral properties are obtained. Practically, because our simulation is always with finitely many degrees of freedom, we consider the average over many samples. See Refs. [55,65] for details.

Models
We will study two physically distinct models. The first will be Sachdev-Ye-Kitaev (SYK) model of four-local Majorana fermions coupled randomly and non-locally. The second model is the XXZ spin chain composed of nearest neighbor spin interactions and random magnetic field along the z-axis.

Sachdev-Ye-Kitaev (SYK)
The SYK model [5,59,60] (see Ref. [67] for a recent review) with N Majorana fermions is given byĤ where the anti-commutation relations are given by and J ijkl is random Gaussian with mean zero and standard deviation J. We set J to be 1, so all times are measured in units of 1/J. The dimension of the Hilbert space is 2 N/2 . We will also consider slightly more general version where K ij is also a random Gaussian variable with mean zero and standard deviation K. At K = 0, this model is 'maximally chaotic' at low temperatures, in the sense that the MSS bound [5,59] is saturated. When K > 0, it is not chaotic at sufficiently low temperature, while it remains chaotic at high temperature [68,69].

XXZ Spin Chain
We will also consider one-dimensional spin chain with periodic boundary condition σ N site +1 = σ 1 . w i is the random magnetic field along zdirection, chosen to be uniform random number between [−W, +W ]. At W 2.75, most of the energy eigenstates are in the many-body localized (MBL) phase [61,70]. (For the physics of the MBL phase, see e.g. [71,72,73,74]; for OTOC calculations see [75,76,77,78,79].) As W is lowered, ergodic phase expands from the center of the spectrum, and gradually the system becomes ergodic except for a small region at low and high energy regions. The boundary between the ergodic and MBL phases can be obscure when the system size is small.

A Definition: Quantum Lyapunov Exponents
whereΠ j is the canonical conjugate ofẑ j . For a given state |φ ,M ij (t)|φ can grow exponentially, and (where * is a conjugate as an operator acting on the Hilbert space) is a natural counterpart of L ij (t) in the classical theory. 3 From this we can define the Lyapunov exponents λ  i (t), because we do not think there is a risk of confusion. We will take |φ to be energy eigenstates.

Lyapunov Exponents in SYK
As a natural counterpart of Eq. (10) in a fermionic system, we can usê We will again take |φ in Eq. (11) to be energy eigenstates. When K is zero, for N not a multiple of eight, the eigenstates of the Hamiltonian are doubly degenerate due to a symmetry. If N is not a multiple of four, the degeneracy occurs within each parity sector. In order to avoid this uncertainty, we will not consider K = 0, instead we consider very small but finite value of K. For each |φ , we obtain N Lyapunov exponents, λ 1 ≤ λ 2 ≤ · · · ≤ λ N .

Lyapunov Exponents in XXZ
Some caution is in order when we define the Lyapunov exponents for the XXZ spin chain. Because the total z-spin S (total) z = 1 2 i σ z,i commutes with the Hamiltonian, it is better to defineM so that it commutes with S (total) z as well, in order to avoid a mixture of different spin sectors which may complicate the analysis. Furthermore σ x , σ y and σ z are redundant, in the sense that σ x σ y = √ −1σ z . One possible option is to use the fermion representation obtained by the Jordan-Wigner transformation,ψ which satisfy the standard anticommutation relation Still, it is probably better ifM is compatible with the locality which is manifest in terms of the σ i variables. Then there are several other options such as and where σ ± = σx+iσy 2 . Note that the latter is neither Hermitian or skew-Hermitian. The former vanishes at t = 0, which makes it difficult to define the Lyapunov growth precisely at early time, while the latter gives [σ +,i (0), σ −,j (0)] = σ z δ ij . Here we use the latter: Then L (φ) The physical interpretation is clear: σ +,i (t) and σ −,j (0) creates/annihilates the z-spin at point i, j and time t, 0, respectively.

Kolmogorov-Sinai Entropy
For classical systems, the Kolmogorov-Sinai entropy h KS is defined as the sum of all positive Lyapunov exponents. We use the same definition here, by using the Lyapunov exponents we have defined in this paper: Away from the classical limit, the properties of this quantity are not immediately clear. It will be studied in Sec. 5

Is the 'perturbation' actually small?
In the study of quantum chaos based on OTOCs, and also in our approach, one usually assumes that the perturbation -multiplication of local operatorsV andŴ -does not change the state too much. This should be the case when the system size is sufficiently large provided W and V are few-body operators. However it is not necessarily the case in actual numerical calculations. Below we see that, with realistic system sizes within our reach, this assumption is not valid for the SYK model and hence some care is needed when we extract physics from the numerical data.

SYK
When an energy eigenstate |E is 'perturbed' by the multiplications ofψ's, the 'perturbation' is small when the energy is still well localized around E after the 'perturbation'. To see it quantitatively, we take |E to be the ground state |E 0 , and plot 16 and 20 in Fig. 3. Here the energy eigenvalues are ordered as We can see rather large deviations of d 1 (j) and d 2 (j) from 1 even at large values of j, which means thatψ k |E 0 andψ kψk+1 |E 0 involve large contributions from the excited states. When K is close to zero, they are almost uniform superposition of all eigenstates. Therefore, with the resources we used for this paper, it is difficult to study physics at different energy scales separately, especially when K is close to zero.
According to Ref. [68,69], the model under consideration is integrable at sufficiently low temperature, when K > 0 and N = ∞. As we have seen already, it is hard to study the properties of the integrable phase and chaotic phase separately. However by varying  the value of K we can change the numbers of integrable and chaotic states; as K becomes larger, more energy eigenstates belong to the integrable sector. Hence we can learn about the difference between two phases by looking at the way the property of the mixture changes.

XXZ
In order to estimate the size of the perturbation, we calculated Here |E i,s is the energy eigenstate in the total spin s sector, ordered as E 0,s ≤ E 1,s ≤ · · · . The results are shown in Fig. 4. We can see that the perturbations are actually small, unlike the case of SYK.

Relations to Other Approaches
Usually the Lyapunov exponent is defined in terms of OTOC, whereV andŴ are arbitrary local operators and · β stands for the thermal average with temperature T = β −1 . The idea is that the Lyapunov growth with the largest exponent can be captured byM ij with any combination ofẑ i andΠ j , and we do not need to take a specific basis. By assuming the eigenstate thermalization hypothesis, the thermal average should be indistinguishable from the expectation value taken in a typical energy eigenstate, and then the components of L Averages over 100 samples (N site = 6, 8) and 5 samples (N site = 10). The horizontal axis is j/L Sz=1/2 , where L Sz=1/2 is the dimension of S z = 1/2 Hilbert space.
grow as e 2λ i t , and the largest eigenvalue should grow as e 2λ N t . At late time, the largest exponent dominates the growth. Hence the largest Lyapunov exponent obtained by using our definition should be the same as the Lyapunov exponent defined from OTOC at sufficiently late time.
One might consider i L ; if a given theory were gauge theory, it would be a gauge invariant quantity. Actually the spectrum λ i can naturally capture the sub-leading contributions in the Lyapunov growth, because In Ref. [27], an operator which is essentially the same as i L (φ) ii (t) has been considered, and the interpretation as the growth of a size of local operator has been explained.
Note that, even when one is interested only in the largest Lyapunov exponent, a use of the spectrum can have technical gain for numerical calculations. If one uses the usual OTOC to define the Lyapunov exponent as e 2λ (OTOC) t ≡ 1 the contribution from the smaller exponents can have non-negligible contribution at early time. Numerically, it is not easy to study sufficiently late time where the largest exponent dominates. By calculating the spectrum, it is possible to extract the largest exponent even at early time. We will demonstrate this in Sec. 5 ( Fig. 9 and Fig. 13).

Lyapunovian
A related notion called the 'Lyapunovian' was recently proposed in Ref. [53]. Suppose we have a system with onex and oneΠ. By using energy eigenstates |E m , we can make the 'Lyapunovian matrix' Then one can study the eigenvalue statistics of this matrix. A natural counterpart for the many-body case is This is complementary to our approach, in which we tried to see the fixed-energy physics for each reference state. Note also that the size of this Lyapunovian matrix in the many-body case is exponentially large, unlike our approach where the size of the matrix is of order the number of degrees of freedom. We can also consider a hybrid, Then there is no ambiguity associated with a choice of a reference state |φ . It would be interesting to study the properties of the spectrum obtained from this matrix.

Lyapunov spectrum from projection
Another recently proposed approach due to Ref. [52] considers projecting the manybody Schrodinger equation onto a subspace of states in the full Hilbert space, specifically a set of low-entanglement matrix product states. The resulting projected dynamics can be viewed as a classical nonlinear dynamical system with a symplectic structure. As such, it can exhibit classical chaos and has a notion of Lyapunov spectrum. This auxiliary classical problem gives another way of associating a Lyapunov spectrum with an arbitrary quantum system.
Based on the results of Ref. [26], we expect that this projected low-entanglement dynamics can accurately capture the long-distance early growth of OTOCs. Hence, one might expect that part of their spectrum agrees with our definition, although we have not checked this. However, there has also been some work advocating caution with such an approach [86]. More generally, it is not clear to us how their full spectrum relates to the spectrum we defined. It would be interesting to determine if their spectrum also exhibits random matrix statistics, as suggested by our results.

Lyapunov Growth
In this section, we present numerical results for the Lyapunov growth in SYK and XXZ models. We then compare the Kolmogorov-Sinain entropy growth to the entanglement entropy growth inspired by classical analogy of KS entropy and coarse grained entropy.

Lyapunov Growth in SYK
Let us start with the 'usual' OTOC (20), which is In the left panel of Fig. 5, we plot λ (OTOC) t at β = 1/T = 0. We can see the exponential growth followed by the saturation at λ (OTOC) t ∼ log N as in Fig. 6. 4 Two red vertical lines show the times at which this growth reached 20% and 80% for N = 20. Between them, the slope is approximately constant for each N . We can see it more clearly in the middle panel, where λ (OTOC) is shown. At 6 ≤ N ≤ 24, the exponent λ (OTOC) changes substantially with N , and it is not easy to take the large-N limit. In the right panel, we have shown the energy dependence of λ (OTOC) , by taking |φ to be the energy eigenstates. Strangely, the exponent in the ground state is larger than the one in the excited states. As we will explain later, this is because the finite-N effect is large. Next let us see the Lyapunov spectrum obtained from our definition. In the first two panels of Fig. 7, we have plotted λ N t for several values of K and N = 14, 16 at β = 0; we can see the exponential growth followed by the saturation. The third and fourth panels in Fig. 7 show the N dependence of λ N t and λ N for K = 0.01. Compared to λ (OTOC) , the value at each N is larger (see Fig. 9), because λ (OTOC) contains the 'contamination' from the smaller exponents λ 1 , · · · , λ N −1 . As we will see shortly, the finite-N corrections to the smaller exponents are larger than the one to the largest exponent. For this reason, λ (OTOC) depends more severely on N . Note also that the N -dependence of λ N is not smooth, but rather it shows sensitive dependence on N mod 8, which suggests that λ N captures the finer detail of the theory at finite N .
In Fig. 8, all the exponents λ i are shown for N = 16. When K is small, all the exponents are positive. For larger K (K 10), beyond t 1, a gap emerges between the larger half and smaller half of the exponents, and the density distribution of the lower eight exponents become increasingly sharper and get closer to zero. Two red vertical lines represent the 20% and 80% saturation of λ N t. Between them the exponents are almost constant.
In Fig. 9, the largest exponent λ N is compared with λ (OTOC) . As explained around Eq. (20), λ (OTOC) is contaminated by the smaller exponents; indeed, we can see clear dif- 4 The fit value is close to e 2λ (OTOC) t = N 2 . This is value can be explained as follows. At late time, all φ|{ψ i (t), ψ j (0)} 2 |φ in (24) give the same contribution.
Each of them contains four terms, φ|ψ j (0)ψ i (t)ψ i (t)ψ j (0)|φ , φ|ψ i (t)ψ j (0)ψ j (0)ψ i (t)|φ , φ|ψ i (t)ψ j (0)ψ i (t)ψ j (0)|φ and φ|ψ j (0)ψ i (t)ψ j (0)ψ i (t)|φ . The first two terms are 1 4 , while the latter two terms are suppressed at large N . Hence e 2λ (OTOC) t → 1 N N i,j=1 2 × 1 4 = N 2 up to a small correction.  ference in the plot. It is interesting to note that the difference between λ N and λ (OTOC) becomes smaller as N increases. In Sec. 5.1.3, we will study this point further. Fig. 10 is made in order to see the effect of the choice of the reference state. The left panel is the energy vs the largest exponent λ N at N ≤ 24. There are two peculiar points here (note that we observed the same for λ (OTOC) in Fig. 5): The growth rate does not seem to depend heavily on the choice of the energy eigenstate. Even when we take |φ to be the ground state, we can still see the exponential growth. Furthermore, the ground state gives faster growth. Seemingly it has a tension with the large-N result at low temperature, λ = 2πT . Presumably this is because the values of N studied here are so small that the energy excitation caused by operatorM is not negligible, as we have seen in Sec. 4.4. At least, as we can see from the right panel, the exponents calculated by using the ground state |E 0 become smaller than the ones obtained from the state at the center of the energy spectrum |E L/2 at N 26. Still, our numerics is not good enough to show whether λ N defined from the ground state vanishes in the large-N limit, as expected from both the usual intuitive picture and the MSS bound.

Kolmogorov-Sinai and Entanglement Entropy
For classical systems, the Kolmogorov-Sinai entropy h KS is defined as the sum of all positive Lyapunov exponents. We use an analogous definition here for the quantum Lyapunov exponents we have defined in this paper; see Eq. (17).

Kolmogorov-Sinai Entropy (KS)
In the middle panel of Fig. 10, the energy dependence of h KS /N is shown. At small N , the curve is concave due to the finite-N effect. As N becomes large, it will become convex   as we can see from the right panel of Fig. 10.

KS vs EE
As discussed in Sec. 2.1, in the classical limit the growth rate of the coarse grained entropy should agree with the Kolmogorov-Sinai entropy h KS . Therefore, we expect that, in quantum theories, N S EE /|A|, 6 which corresponds to the coarse grained entropy of the system, and h KS t should agree up to an additive constant. As we have emphasized, for the values of N studied in this paper, h KS captures contributions from almost all energy eigenstates.
In the left panel of Fig. 11, we have plotted N S EE /|A| and h KS t obtained from the SYK model at β = 0. At early time, they show similar growths; indeed, as we can see the right panel, at 1 t 2 they agree very well just by a constant shift. 7 This results are not 5 Here we are defining fermions in terms of spin variables via the Jordan-Wigner transform. The spin Hilbert space then has a sensible tensor product structure and this decomposition is what we are referring to. 6 Note that we are proposing to use the entanglement entropy as a quantum analog of classical course grained entropy. 7 This shift can be understood as the ambiguity of the size of the cell in Fig. 1.
conclusive, however, they are suggestive that the KS entropy can actually be understood as the entropy production rate in quantum systems.  for all j = 1, 2, · · · , N 2 has been used.
[Right] Essentially the same plot, but S EE is shifted by a constant.

Fastest entropy generator?
As we have seen in Sec. 5.1.1, the difference between λ N and λ OTOC becomes smaller as N increases. This means the largest exponent λ N and the smallest exponent λ 1 get close. We can actually numerically confirm that λ N − λ 1 scales as 1/N at t 2. This strongly suggest that the Lyapunov spectrum peaks like the delta function in the large-N limit. This is consistent with the previous analysis on the OTOC at N = ∞: if the spectrum has nontrivial distribution, it can give a power law correction t ν e λ OTOC t with ν > 0. But such correction has not been found [59].
It is interesting to compare this behavior with the usual (weakly coupled) string theory dual. The Lyapunov spectrum of the D0-brane matrix model in the classical limit (highly stringy region) converges to the semi-circle distribution with an O(N 0 ) width [65], unlike the weak coupling region of SYK. Hence the large-N limit of the D0-brane theory does not by itself make the distribution peaked, unlike in the SYK model. However the analyses on the dual gravity side including the finite coupling correction at large N (the α correction in gravity side) to the Lyapunov exponents [43] suggest that the Lyapunov exponents peak at the largest possible value (the MSS bound) at strong coupling (for example, a power-law correction has not been seen). The same seems to be true for other quantum field theories which can be analyzed with dual gravity calculations [43].
Assuming this is true, both in the SYK model and quantum field theory with the usual string theory dual, all exponents saturate the MSS bound at strong coupling and large N . Therefore, both of them appear to take the largest possible Kolmogorov-Sinai entropy, or  the largest possible entropy production rate. 8 For the canonical ensemble, the largest possible value of each exponent is 2πT . Apparently, when all exponents saturate this bound, the sum is maximal. Hence, it would be natural to conjecture that black hole has the largest possible KS entropy. Note that we have the single black hole configuration (the leftmost figure in Fig. 2) in mind.
For the microcanonical ensemble, it would be natural to conjecture that the entropy generation rate increases as the black hole grows, as demonstrated for a simple case in Sec. 2.1.1. We can also show the same pattern for more generic initial conditions, and we can also show that the KS entropy decreases as black hole evaporates [54]. Therefore we conjecture that the KS entropy is maximal when all the degrees of freedom are absorbed in one black hole and thermalized.

Lyapunov Growth in XXZ
Because the total z-spin S (total) z commutes with the Hamiltonian (9), we focus on the Lyapunov growth in the zero-spin sector, S (total) z = 0. As we have seen in Sec. 4.4, the multiplication of σ can actually be regarded as a small perturbation, and hence it makes sense to study the temperature dependence, unlike the case of SYK.
In Fig. 12 we have plotted λ N site t and λ (OTOC) t as functions of t for N site = 12 and λ N site t for various temperatures. For each N nsite , λ (OTOC) t converges to the same value,   Fig. 14, the deviation from the late-time plateau in the MBL phase is plotted in the log-log scale. We can see a straight line, which means the late-time behavior is A − Bt −p . This is consistent with the theoretical expectation in Ref. [45]. With the range of N site available at this moment, it is hard to take the large volume limit, N site → ∞.
A possible explanation of this pattern is as follows. In the classical theory, the perturbation at t = 0 can be sent arbitrarily small, and the exponential growth can continue forever. However when the perturbation is finite, the exponential growth stops at a finite time; otherwise the causality is broken! (In the nonrelativistic theory the speed is not lim-9 For each energy eigenstate |E , terms of the form E|σ +,j (0)σ −,i (t)σ +,i (t)σ −,j (0)|E = ||σ +,i (t)σ −,j (0)|E || 2 and E|σ −,i (t)σ +,j (0)σ −,j (0)σ +,i (t)|E = ||σ −,j (0)σ +,i (t)|E || 2 give dominant contributions at late time. Because we are taking |E to be in the total spin zero sector, when σ −,j (0) is multiplied on |E , half of the terms in the z-spin basis -terms with down spin at j-th sireis annihilated. Hence σ −,j (0)|E is roughly norm 1/ √ 2, and consists of terms with N site /2 + 1 down spins and N site /2 − 1 up spins. Then when we further multiply σ +,i (t), (N site /2 + 1)/N site terms survive. Hence ||σ +,i (t)σ −,j (0)|E || 2 (N site /2 + 1)/2N site . For the same reason, ||σ −,j (0)σ +,i (t)|E || 2 (N site /2 + 1)/N site . Hence e 2λ (OTOC) t  ited, but once the initial condition is set, then the energy conservation sets the upper limit of the speed at later time.) In the quantum theory, the perturbation is necessarily finite and hence the exponential growth has to stop at some point. In nonlocal systems like the matrix model and SYK model, the exponential growth stops when the 'local' perturbation (say the multiplication of ψ 1 ) affects all other degrees of freedom substantially. This is typically O(log N site ) time. For local systems like the spin chain, the exponential growth has to stop at O(1) time, because σ j (t) and σ i (0) commute when t is smaller than |i − j| divided by the butterfly velocity; when the exponential growth stops, only O(1) number of degrees of freedom talk to each other. (The system size does not matter, otherwise the causality or the Lieb-Robinson bound is broken.)M ij (t) is a banded matrix with width w ∼ t, and if we use a very rough approximation that all the nonzero entries are of order one, then the singular values scale as √ w. It leads to a late time behavior λt ∼ 0.5 log t, which is in the right ballpark compared to 0.35 log t and 0.32 log t in the left panels of Fig. 14. As a related example, let us consider planar black p-brane (p > 0), which is described by U(N ) super Yang-Mills on R 1,p . How is a localized perturbation scrambled in this theory? Firstly the fast scrambling with λ ∼ 2πT mixes the information among the gauge degrees of freedom; then the information gradually spreads along R p . In terms of gravity, the horizon has a topology of S 8−p ×R p , and the fast scrambling takes place along S 8−p while the growth along R p is slower and dominant at late time. The 1d spin chain is analogous to p = 1 and very small N .
In the explanation above, only the local physics is important for the early-time exponential growth. Hence the same pattern is expected both in the ergodic and MBL phases. The time scale of the saturation of the power growth can be different; in the ergodic phase the saturation time scale should increase with the system size, while in the MBL phase it is independent of the system size, instead only the volume of the region affected by the perturbation matters. It is consistent with the numerical results: in Fig. 14 the saturation time scale increases with the system size, while in Fig. 14 it seems to be insensitive to the system size.
For local quantum systems, the absence of the exponential Lyapunov growth should be generic. As we have seen above, it is not easy to distinguish the ergodic and MBL phases just from the power growth. However, as we will see in Sec. 6.2, the statistical features of the Lyapunov exponents are clearly different in these two phases.

Random Matrix Statistics of Lyapunov Spectrum
In this section we study the statistical properties of the quantum Lyapunov spectrum, motivated by the universality in the classical Lyapunov exponents explained in Sec. 2.2.

Lyapunov Spectrum vs RMT in SYK
According to Ref. [68,69], the q = 2 deformed SYK model (Eq. (8)) is integrable at sufficiently low temperature, when K > 0. Therefore, by carefully choosing the energy eigenstates, we can study the statistical features of the Lyapunov spectra in the integrable and chaotic phases, in principle. However, as we have seen already, perturbations by multiplications ofψ is not really 'small' at the system size we can study numerically, and hence it is hard to study the properties of the integrable phase and chaotic phase separately.
Below, in addition to the nearest-neighbor level correlation, we study the nearestneighbor gap ratio, In the upper row of Fig. 16, we have shown the nearest neighbor level separation estimated after the standard unfolding procedure. Namely, we have estimated the distribution of the Lyapunov exponents by using energy eigenstates within certain range (between 5% and 10% in these specific plots), numerical fit it by polynomial of degree 10 and used it for the unfolding. We can see good agreement with GUE at small K, but there is a small but visible deviation.
A possible flaw of this method is that, when N site is small, peaks arising due to the level repulsion can be visible and the unfolding can eliminate them as well, so that the universal random matrix behavior is erased. To circumvent such possibility, we tried another unfolding prescription as well: normalize the gap g i = λ i+1 − λ i so the average is 1. Definẽ g i ≡ g i / g i and look at the distribution ofg i . We call it 'fixed-i unfolding'. The result is shown in the lower row of Fig. 16. Compared to the standard unfolding (the upper row of Fig. 16), the agreement with GUE is improved substantially for K = 0.01. On the other hand, for K = 10, there is no improvement.
In Fig. 17, we show the time dependence of r at β = 0. (As discussed around Fig. 8 for the case of N = 16, a gap develops between λ N/2 and λ N/2+1 for larger K as t is increased. In order to check that this gap does not affect the result, we have also calculated r using only the larger N/2 exponents. We can see that the early rime behaviors are almost identical.) At K = 10 and 100, large deviation from the GUE value is observed before the Lyapunov growth (with almost constant Lyapunov exponent) sets in. At K = 1, a large deviation is observed before the Lyapunov growth ends. At K = 0.01 and 0.1, the r-parameter stays close to the GUE value even at t = 100. In Fig. 18, we show essentially the same plot, but using |E 0 , |E L/2 and |E L−1 . We don't see a significant change as expected; for the value of N we study, the perturbation is too large, so that we can only see a mixture of almost all states. At sufficiently large N we expect different behaviors depending on the energy.
As K becomes larger, more energy eigenstates belong to the integrable sector. That the Poisson statistics sets in with larger K suggests that the spectrum in the integrable sector follows Poisson. Note that the GUE appears as t → 0 even at large K. We do not have understanding about this property. [Lower] The fixed-i unfolding has been performed. When K is small, good agreement with GUE is observed, until very late time. When K is large, small deviation from GUE can be seen at early time, and the deviation grows quickly at later time.  The r-parameter r with fixed-i unfolding, obtained using all exponents for the N site = 14 SYK model, for the ground state (left), the state at the center of the energy spectrum (center) and the highest energy state, plotted as a function of t for K = 0.01, 0.1, 1, 10.

Lyapunov Spectrum vs RMT in XXZ
Because the total z-spin S (total) z commutes with the Hamiltonian (9), we focus on the Lyapunov spectra obtained by using eigenstates with S (total) z = 0. As explained in Sec. 3.2, this theory has ergodic and MBL phases. We will study W = 0.5 and 1.0, which are mainly in the ergodic phase except for the edges of the energy spectrum, and W = 4.0, which is dominantly in the MBL phase.
In Fig. 19 we have plotted P (s) obtained by using eigenstates with the energy within 45% -55% from the lower edge of the spectrum, 10 with the fixed-i unfolding introduced in Sec. 6.1. We can see a good agreement with GUE for W = 0.5 and 1.0 (the ergodic phase) at sufficiently late time, while the Poisson distribution is favored for W = 4.0 (the MBL phase). In Fig. 20, we fixed W = 0.5 and varied the energy band. We can see the time evolution strongly depends on the choice of the energy. We can also see that the GUE is not obtained near the ground state, which is close to the MBL phase. Note that we have shown two results, one obtained by using all exponents and the other obtained by only the largest three exponents. The reason is as follows. When we plot the sample averaged values of λ i against the eigenstate index, for the XXZ model with smaller values of W , at short times large gaps between the smaller, nearly twofold degenerate exponents are observed for lower energy eigenstates. For larger N/2 exponents, the averaged values are evenly distributed for t 5. Therefore, in order to make sure the universal behavior can be seen regardless of the choice of the exponents, we have shown two results. Below, we will also study the gap ratio r. This is more sensitive to the change of the gap size, and hence, we will use the largest three exponents for safety.
In order to see the time and energy dependence quantitatively, we have plotted r with the fixed-i unfolding in Fig. 21, for W = 0.5, for various energy bands. We can see better agreement with GUE (both the value and time window) at the center of the energy spectrum. Recalling that the middle of the energy spectrum is in the ergodic phase except that the small region near the edges remain MBL, this is consistent with the interpretation that GUE is obtained for the ergodic phase but not for the MBL phase.
The N site -dependence of r is shown in Fig. 22. The GUE can be seen with good precision at N site ≥ 8. We studied the values of r for N = 6, 8, 10 and 12 until very late time (t 10 8 ). For these values of N , r becomes almost constant, which is different from the GUE value, at t 10. Therefore we expect the importance of the large-N limit before t → ∞ for the emergence of the universality.

Conclusion and Outlook
In this paper we have proposed a generalization of the Lyapunov spectrum to quantum theories, and studied its properties by using the SYK model and the XXZ model with    [upper] for the two gaps between the three largest exponents and [lower] for all gaps have been used for the Lyapunov exponents obtained using 0% -5% (left), 45% -55% (middle; the same as the left panel in Fig. 19), and 85% -90% (right) of the spectrum for each sample.  Figure 21: Plot of r for the two gaps between the three largest exponents for XXZ, N site = 10, W = 0.5. Better agreement with GUE (both the value and time window) can be seen at the center of the energy spectrum, which is consistent with the fact that the center of the energy spectrum is in the ergodic phase while the edges remain MBL.  random magnetic field as examples. By definition, the Lyapunov spectrum contains more information than just the largest exponent.
The KS entropy -which we defined by the sum of the positive exponents -is likely to be a better characterization of the strength of the chaos, because it can describe the entropy production rate. We conjectured that the black hole maximized the KS entropy.
We also found the numerical evidence for the universality of the Lyapunov spectrum. (Previously, this universality has been observed in classical chaos as well [55].) It is interesting if we could understand the meaning of the onset of the universal RMT behavior. It should have something to do with holography, because special theories which are dual to quantum black holes -the SYK model, and classical D0-brane matrix model, as demonstrated in Ref. [55] -show the universality already at t = 0. We propose that the quantum systems holographically dual to Einstein gravity satisfy this 'strong' universality.
The conjectures above are based mainly on the numerical observations for limited number of theories. It is important to study more examples, and also, to develop the understanding on the gravity side. Another important issue is how the universality class is determined. For the examples studied in this paper we observed only GUE ensemble, regardless of the system size.
It is also important to apply the method presented in this paper to various physical systems, especially in the contexts of condensed matter and quantum gravity. It should be possible to get new insight into scrambling and thermalization by observing the Lyapunov growth, and with various examples we might be able to understand the meaning of the characteristic time scales associated with the Lyapunov growth and the onset of the universal spectral behavior. Toward the study of full string theory, probably the weakly coupled region of D0-brane matrix model [87,88] is a good place to start.
It would also be interesting to develop a measurement protocol for the Lyapunov spectrum along the lines on Ref. [13]. A brute force way to approach the problem is to consider performing a whole set of many-body interference experiments to measure the various matrix elements needed to construct the spectrum-defining matrix. A detailed study of the feasibility of the this approach, and the search for more economical approaches, is left to future work.
String Theory 2018" were useful to complete this work. We acknowledge JSPS KAKENHI Grants 17K14285 (M. H.) and 17K17822 (M. T.). This work was also partially supported by the Simons Foundation via the It From Qubit Collaboration, and by the Department of Energy, award number DE-SC0017905.