Chaos in classical D0-brane mechanics

We study chaos in the classical limit of the matrix quantum mechanical system describing D0-brane dynamics. We determine a precise value of the largest Lyapunov exponent, and, with less precision, calculate the entire spectrum of Lyapunov exponents. We verify that these approach a smooth limit as N → ∞. We show that a classical analog of scrambling occurs with fast scrambling scaling, t∗ ∼ log S. These results confirm the k-locality property of matrix mechanics discussed by Sekino and Susskind.


Introduction and Summary of Results
This paper is devoted to a study of classical chaos in the classical limit of the matrix quantum mechanical system describing D0-brane dynamics.In particular we compute Lyapunov exponents in this system.
The motivation for this work flows from recent progress on the overlap between quantum chaos and quantum gravity.These developments have their origin in Quantum Information theory, and specifically in work done making good approximations to random unitary operators [1][2][3][4][5][6].Such approximations can be implemented very quickly, in a time proportional to log n, where n is the number of qubits (the analog of the entropy S in this system).
Hayden and Preskill [7] connected this timescale to one characteristic of black hole horizons [8,9] t * ∼ R log(M/m p ) ∼ R log S, where R is the Schwarzschild radius, M is the black hole mass, m p is the Planck mass and S is the Bekenstein-Hawking entropy of the black hole.This logarithm is a consequence of the exponential blueshift of modes at late Schwarzschild time near the horizon, following from its Rindler structure.They presented an example of a model typical of those discussed in Quantum Information: a Hamiltonian coupling pairs of qubits nonlocally with a random pattern, with the 2-qubit gates being chosen at random.It is easy to see that such a Hamiltonian will cause all qubits to become entangled with each other in a time of order log n, and reasonable to conjecture that chaos has set in by this time [1][2][3][4][5][6][7].This conjecture is supported by analysis of quantum circuits and a Lieb-Robinson bound [10].A crucial aspect of such Hamiltonians is "k-locality," where interactions can be nonlocal but only a finite number k of qubits are coupled together in each term of the Hamiltonian, independent of the total number of qubits in the system.
Sekino and Susskind made the connection between these ideas and gauge/gravity duality [11].They argued that matrix quantum systems behave similarly to k-local qubit systems: the matrix indices function like the qubit label, and the sum over indices in the matrix interactions couples a finite number of index pairs together nonlocally, but satisfying the k-local property.In some ways the simplest such system is maximally supersymmetric matrix quantum mechanics [12], which has M-theory as its infrared gravity dual in Matrix Theory [13] and type IIA string theory at somewhat higher energies [14].The horizons of the black hole duals in such systems are Rindler in nature, and so matrix quantum systems have the characteristic logarithmic time which they interpreted as a "scrambling time" t * ∼ β log S (here β ∼ R is the inverse Hawking temperature of the black hole).Sekino and Susskind went on to make the "fast scrambling conjecture," that in all reasonable physical systems chaos cannot set in faster than the logarithmic rate it does in black holes, in a time t * ∼ β log S.
The next stage in the analysis of this kind of quantum chaos was undertaken in [15][16][17][18][19][20]24] using holographic (and other) techniques.A sharp diagnostic of chaos is the growth of a commutator [21,22] of simple operators with time, C(t) = − [V, W (t)] 2 , where the brackets denote thermal expectation value.In a chaotic quantum system W (t) (in these contexts sometimes referred to as a "precursor" [23]) becomes more complicated with time due to the lack of cancellation between the first and last factors in W (t) = e iHt W e −iHt induced by the small perturbation W .On expanding out the commutator one finds that the quantity most sensitive to chaos is an out-of-time-order correlator, D(t) = V W (t)V W (t) .
As Larkin and Ovchinnikov [21] pointed out long ago, in few body quantum systems described schematically by a coordinate q and momentum p the commutator C(t) = −[p, q(t)] 2 goes over in the semiclassical limit to C(t) → 2 {p, q(t)} 2 where {•, •} is the Poisson bracket.This can be expressed as 2 ∂q(t) ∂q(0) 2 = 2 e 2λ L t , where λ L is the Lyapunov exponent.This motivates using the commutator as a diagnostic of chaos.
The quantities C(t), D(t) (and closely related thermofield double two-sided correlators) have been computed holographically in [15][16][17][18][19].The essential bulk phenomenon is a high energy collision between the quanta created by V and W (t) near the horizon.The perturbative strength of gravitational scattering in such a collision is of order G N s (in AdS units) where G N is Newton's constant and s is the center of mass energy squared.The center of mass energy is (up to order one constants) s = 1 β 2 exp 2πt β because of the Rindler nature of the horizon and the role of boundary time as Schwarzschild time.In the Einstein gravity limit the first term surviving in the commutator is second order in This becomes of order one at This is the precise large N holographic scrambling time for systems with a bulk Einstein gravity dual.Kitaev [19], building on [21], connected the exponential time behavior in (1) to Lyapunov behavior in chaotic systems.Here the Lyapunov exponent is given by λ L = 2π β = 2πT . 4This exponential behavior and the small 1 N 2 prefactor are the ingredients determining the fast scrambling conjecture timescale.
The authors of [20] were able to establish the Einstein gravity value λ L = 2π β = 2πT as a sharp upper bound on thermal quantum systems with a large number of degrees of freedom and a large hierarchy between scrambling and dissipation times.The argument uses only general principles of quantum mechanics and plausible physical assumptions about the decay of time-ordered correlators.This bound does not enable one to compute the value of Lyapunov exponents in a given quantum system.A suggestion was made in [17] about how to compute λ L at weak coupling, motivated by the BFKL ladder summation for high energy scattering.Stanford [25] has recently succeeded in implementing this calculation in matrix φ 4 quantum field theory.
Kitaev [24] has shown how to compute λ L in a strongly coupled large N fermionic quantum mechanics system related to the Sachdev-Ye model [26,27].He proceeds by summing ladder diagrams that in this case give the exact large N solution to the model.In the limit of strong coupling the exponent saturates the bound -a remarkable result.
Direct numerical work on this aspect of quantum gauge systems seems challenging.Here we follow a different approach, exploring the classical dynamics of such a system.In particular we explore the classical dynamics of the maximally supersymmetric matrix quantum mechanics in 0+1 dimensions, in the large N limit.The Lagrangian is Here X i (i = 1, . . ., 9) are N × N traceless Hermitian matrices and is the covariant derivative, where A t is the SU (N ) gauge field.We take the large N limit with the 't Hooft coupling λ = g 2 N .The remaining terms in (3) involve fermions, which do not contribute in the classical limit.
At large N and low temperature, the theory under discussion is holographically dual to a black hole in classical gravity.We will focus on the large N , high temperature classical limit, where the dual black hole description is no longer valid.The dimensionless parameter λ eff = λ/T 3 characterizing the large N dynamics goes to zero in this limit.(Previous numerical studies confirmed that there is no phase transition which separates the low and high temperature regions in this theory [37].We therefore expect some qualitative features of the black hole, such as fast scrambling, to survive at high temperature.) The high temperature limit of a quantum mechanical system is well approximated by its classical dynamics.This statement is only true for quantum mechanics, not quantum field theory -high-temperature field theory does not have a good classical limit because of the UV catastrophe.Indeed, in high-temperature quantum field theory the occupation numbers of typical field modes are of order one, while classical equations of motion approximate quantum fields with large occupation numbers. 5revious numerical studies [33,34,36] showed that for generic initial conditions the classical system thermalizes into what can be thought of as a bound thermal state of N D0-branes.
In this work we compute the Lyapunov exponents of this system by solving the equations of motion numerically.For the leading exponent we give a precise result, while for the spectrum of subleading exponents we get a semi-quantitative estimate.The classical system has a phase space with dimension that is of order N 2 and has the same number of Lyapunov exponents.At large N we find that they converge to a continuous spectrum with a finite maximum value.That the chaotic dynamics has a smooth large N limit provides support for the k-locality of matrix interactions, as discussed by Sekino and Susskind [11].
In particular we find that that the largest Lyapunov exponent λ L approaches a finite value in the large N limit, λ L → 0.292 λ 1/4 eff T .Note that this is parametrically smaller than the bound λ L ≤ 2πT established in [20] in the classical limit λ eff → 0. This determines the fast scrambling time, t * ∼ 1 λ L log N 2 , confirming that this model is a fast scrambler.In classical systems the Lyapunov exponents are related to the Kolmogorov-Sinai (KS) entropy, which measures the rate of growth of coarse-grained entropy when the system is far away from equilibrium.Pesin proved that the KS entropy is equal to the sum of positive Lyapunov exponents, and this result allows us to compute the KS entropy in the matrix theory.Our result that the Lyapunov spectrum converges to a smooth density at large N implies that the KS entropy is proportional to N 2 .
The paper is organized as follows.In Section 2 we present the matrix model and describe its classical limit.In Section 3 we review the classical theory of Lyapunov exponents, and explain how it applies to the classical matrix model.The main difficulty here is in dealing with the gauge symmetry of the model.In Section 4 we present numerical results for the Lyapunov exponent in the large N limit, using various methods to compute the exponent.Then, in Section 5 we present the computation of the Lyapunov spectrum in this system.Section 6 includes a discussion of the results, and several appendices present some technical details of the computation.

D0-Branes at High Temperature
The model we consider is the low-energy effective theory that lives on a stack of N D0branes [38].It can be obtained by dimensionally reducing super Yang-Mills in 9+1 dimensions to zero space dimensions.This is a supersymmetric quantum mechanics with a U (N ) gauge symmetry and an SO(9) global R-symmetry.Its degrees of freedom include nine N × N Hermitian matrices X i ab , i = 1, . . ., 9, a, b = 1, . . ., N , as well as 16 fermions ψ ab in the spinor representation of SO (9), and a gauge field A ab .The action is The covariant derivative is • ], and summation over repeated SO(9) indices is implied.In this work we take the matrices X i to be traceless because the trace mode is decoupled.When the matrices X i are diagonal, their N eigenvalues correspond to the positions of the D0-branes in 9-dimensional flat space.Off-diagonal elements correspond to open string degrees of freedom that are stretched between different branes.
Let us take the large N limit, keeping the 't Hooft coupling λ = g 2 N fixed.The coupling λ is dimensionful, and at finite temperature T we can define the dimensionless coupling λ eff = λ/T 3 which controls the size of loop corrections.We will take the limit of small λ eff , which is the weak coupling / high-temperature limit where classical dynamics provides a good approximation.There, the fermions do not contribute to the dynamics of X i , so we can discard them [39].We choose the gauge to be A t = 0. Integrating out the gauge field leads to the Gauss law constraint, which should be preserved due to gauge invariance.Fixing A t = 0 does not completely fix the gauge; the residual gauge freedom corresponds to global (i.e.time-independent) SU (N ) transformations.
We will work in an ensemble with fixed energy E, and where the conserved angular momentum is set to zero. 6Averages in this ensemble will agree with thermal averages in the thermodynamic limit N → ∞; the corresponding temperature T is given as follows.
The equipartition theorem for this system relates temperature, energy and number of degrees of freedom as The total energy is E = K + U where K is the kinetic energy and U is the potential energy.n dof is the number of physical degrees of freedom.Naively the total number of degrees of freedom is d(N 2 − 1), where d = 9 is the number of matrices, but in accounting for the Gauss law constraint (5) and the residual gauge symmetry we have to subtract (N 2 − 1).Furthermore, the conservation of the angular momentum Tr(X i Ẋj − X j Ẋi ) should be taken into account, reducing the number of degrees of freedom by d(d − 1)/2. Therefore, In the weak coupling limit we can use the classical approximation to describe the real-time evolution of the system, at least for typical states at a given temperature. 7Thermodynamic properties can then be computed using ergodicity, which we assume.(Numerical results are consistent with this assumption.)The scalar equation of motion in our gauge is Equations ( 5) and ( 8) fully describe the system in the classical approximation.Notice that the equations do not depend on the coupling.Therefore, due to the form of the action (4), classical observables may depend on the temperature and the coupling only through the combination λT = λ eff T 4 ; the power of this combination is then determined by dimensional analysis.From now on we set T = 1 without loss of generality.

Discretization
In order to study the time evolution numerically we discretize the equation of motion (8) while preserving the constraint (5) exactly.For this purpose we write the equations of motion as The discretized evolution with time step δt is taken to be of order δt 2 [36].It is given by It is easy to check that this prescription preserves the Gauss law constraint, namely that if the constraint i [X i , V i ] = 0 holds at time t, then under the evolution (11) it also holds at t + δt. 8 All that is left is to ensure that the initial conditions obey the constraint and have zero angular momentum.We do this by initially setting V i = 0 while taking X i to have random (Gaussian) matrix elements.
In order to control the discretization error after evolving for time t, we use two different time steps: δt = 10 −4 and δt = 5 • 10 −4 , and compare the results.We compared several quantities such as the norm of the perturbation |δX|, whose definition will be given later in this paper, as well as Tr(X 2 i ) and Tr([X i , X j ] 2 ).We found agreement for t 60.A similar comparison with the same discretization has been performed previously; see Fig. 2 of [36].

Lyapunov Exponents
In this section we briefly review the theory of Lyapunov exponents in classical systems, and its application to the matrix model.We stress the complexities that arise due to the gauge symmetry of our model.Consider a Hamiltonian system with a phase space M of dimension n.Hamilton's equations define the mapping of a point x 0 in phase space to a point x(t) after time t.By linearizing Hamilton's equations we can define a linear operator U (t; x 0 ) (the transfer matrix), that maps a tangent vector δx 0 (i.e. an infinitesimal perturbation) at x 0 to a final vector δx(t) at x(t).
The signature of a chaotic system is the exponential growth of perturbations.In order to discuss this growth we introduce a Riemannian metric g on phase space.In a chaotic system, a typical perturbation grows as |δx(t)| ∼ |δx 0 |e λ L t , where |δx| = g(δx, δx).We define the Lyapunov exponent that is associated with the initial perturbation by Note that there is no natural choice for g on phase space, but if phase space is compact then the Lyapunov exponents are independent of g; see Appendix A. If phase space is noncompact then the exponents will not be well-defined in general.
In an ergodic system, the Lyapunov exponents λ L can take up to dim(M) = n distinct values [40].The largest exponent is the one that is referred to as 'the' Lyapunov exponent, because it dominates the growth of typical (non-fine-tuned) perturbations.The spectrum of Lyapunov exponents is determined by the size of g (U (t; x 0 )δx, U (t; x 0 )δx) = g δx, U † (t; x 0 )U (t; x 0 )δx , namely by the eigenvalues of U † (t; x 0 )U (t; x 0 ).Equivalently, the spectrum can be determined by performing a singular-value decomposition (SVD) on the transfer matrix; here we choose an orthonormal basis for the tangent space (with respect to g), and write the transfer matrix in this basis as where W, V are unitary and Σ = diag(σ 1 , . . ., σ n ) is positive-definite, with The Lyapunov exponents λ 1 , . . ., λ n are then given in terms of the decomposition by For ergodic systems, λ i = λ i (x 0 ) is independent of the starting point x 0 .Phase space carries with it a symplectic structure (a closed, non-degenerate 2-form ω), and the transfer matrix is a symplectic transformation.Therefore, the Lyapunov exponents are paired: For every exponent λ i there is a corresponding exponent −λ i [41].We will be interested in taking the limit in which the dimension of phase space n goes to infinity (this will correspond to a 't Hooft limit of our matrix model).As we will see, in the matrix model the set of discrete exponents λ i approach a distribution ρ(λ) in this limit.The distribution is supported on [−λ L , λ L ] where λ L is finite.

Finite Time Approximation
In a numerical calculation of the exponent based on (12), time must be kept finite.We define the time-dependent exponents λ i (t; x 0 ) by They converge to the Lyapunov exponents λ i as t → ∞.Due to the symplectic structure, the exponents are paired: λ i (t; x 0 ) and −λ i (t; x 0 ) appear together.
Let δx 0 be a generic perturbation of unit norm.Given the decomposition (13), let {v i (t)} be the column vectors of V (t) such that where w i (t) are the columns of W (t) (from now on the dependence on x 0 will be implicit).
Expand the initial perturbation as δx 0 = i c i (t)v i (t).The evolved perturbation then has squared norm In the last step we used the fact that for a typical vector δx 0 we expect that |c i (t)| 2 ≈ 1/n.The Lyapunov exponent (defined in (12)) is then approximated at finite times by In Hamiltonian systems, it was argued that individual exponents typically approach their asymptotic values as λ i (t) ∼ λ i + a i t after averaging over initial conditions [42]. 9 In the matrix model, it will turn out that the individual exponent evolution is well-approximated by λ i (t) ∼ λ i + a i t + b i log t t .We will also find that the effective exponent λ L (t) approaches its asymptotic value λ L much faster than do the individual exponents.

Matrix Model Application
Let us now consider the Lyapunov exponents in the context of the D0-brane system.The phase space M of this system (after gauge fixing to A t = 0) is a vector space with coordinates (X i , V i ) and a symplectic structure that is given by ω = dX i ab ∧ dV i ba .As explained above, in order to have well-defined Lyapunov exponents the space should be compact.Let us restrict ourselves to a subspace with fixed energy.This subspace is still noncompact due to the existence of flat directions: A configuration of the form has energy that is independent of the brane position y i .However, as we show in Appendix E, simple estimates suggest that even for finite N the fixed-energy phase space has finite volume, and this is confirmed by the equilibration of the system even for N = 2. 10Therefore, the Lyapunov exponents are effectively well-defined for this classical system.
The next problem we face is that of gauge redundancy.Having fixed the gauge to A t = 0, all physical configurations must satisfy the Gauss law constraint (5).Let us restrict our space to the constraint surface11 Restricting to M 0 is not sufficient because M 0 is not a phase space (in general it does not admit a symplectic structure), and also because of residual gauge symmetries.To see this, let us define a Riemannian metric g on the phase space M by g(δx, δx ) = g(δX, δV ; δX , δV ) ≡ Tr(δXδX ) + Tr(δV δV ) .
Here, δx = (δX, δV ) denotes a vector in phase space.This metric is invariant under the residual gauge transformations (with respect to the gauge A t = 0) However, the metric ( 21) leads to a non-zero geodesic distance between gauge-equivalent configurations, namely between two configurations that are related by the transformation (22).Therefore, using the phase space M (or the constrained space M 0 ) with the metric (21) to define the Lyapunov exponents will lead to 'spurious' exponents that correspond to pure gauge modes rather than to physical perturbations.
The solution to this problem is to define a physical phase space from which the pure gauge modes have been modded out.This procedure is known as the symplectic reduction of M, and it is explained in detail in Appendix B. The upshot it that the physical Lyapunov exponents are obtained from a modified transfer matrix given by where P (x) is a projector that projects out vectors that do not obey the Gauss law constraint, as well as vectors that correspond to pure gauge transformations.The gaugeinvariant exponents are obtained as before by a singular value decomposition of U phys .
The presence of residual gauge transformations does not affect the leading exponent, essentially because perturbations corresponding to gauge transformations do not grow with time.In the following section we will compute the leading exponent, so we will be able to ignore this issue.In Sec. 5 we will compute the full spectrum of exponents, and there the prescription (23) will be used.

Leading Exponent Computation
In this section we compute the leading Lyapunov exponent of the classical matrix model by following diverging trajectories in phase space.Our main result is that the exponent converges at large N .One important corollary is that the classical matrix model is a fast scrambler, namely that the classical analog of scrambling time (defined below) scales as log N 2 .Finally, we compute the exponent using an alternative method by considering gauge-invariant correlation functions, and find good agreement.
The computation of the Lyapunov exponent consists of three steps.
1. 'Thermalize' the system by evolving it for a long enough time.
3. Evolve both the original and the perturbed configurations, measuring the exponential rate at which they diverge.
Let us discuss each step in detail.
We begin by choosing an initial state where the X variables are random and traceless, and where Ẋ = 0.This initial state satisfies the Gauss law constraint, and also has vanishing momentum and angular momentum.We then evolve the system for a sufficiently long time, so that it reaches a 'typical state' that is uncorrelated with the (atypical) initial conditions.This is the ergodic equivalent of thermalization.How long do we need to evolve for in order to thermalize the system?Fig. 1 shows the resulting Lyapunov exponents as a function of thermalization time t 0 .(We will explain how the exponents are evaluated shortly.)We see convergence for t 0 2000, and in what follows we set t 0 = 4000.Note that this is much longer than the thermalization time typically needed for other observables, and for observables previously studied in the literature; see e.g.[33,36].The origin of this slow relaxation phenomenon is mysterious and is an interesting topic for future research.
Given a thermalized configuration (X, V ), we perturb it slightly while preserving the Gauss law constraint (5) by using the method described in Appendix C. Having obtained the reference configuration (X, V ) and the perturbed configuration (X , V ), we evolve both together and compute the distance between them.The distance function we use is The distance grows exponentially, where λ L is the Lyapunov exponent. 12he evolution of |δX(t)| is shown in Fig. 2. Exponential growth sets in quickly, and continues until the size of the perturbation becomes of the order of the system size, at around t 60.We shall call this the 'scrambling time' t * of the perturbation.In principle, the Lyapunov exponent can be extracted directly from the exponential growth.As discussed in Sec.3.1, the accuracy of this calculation is limited by the finite time of the perturbation growth.For this reason we now consider Sprott's algorithm [43], which is an alternative method for computing the exponent.The algorithm is explained in Appendix D. It allows us to study the growth at arbitrarily long time scale (we used t = 10 4 ), and to extract the largest Lyapunov exponent.Fig. 3 shows the convergence of the exponent computed using this algorithm.Notice that convergence to within one percent occurs at t 100.This suggests that the Sprott result should be close to the direct fit result, and this is indeed what we find.In Sec. 5 we determine the subleading exponents and give more detail on how this agreement is achieved.
The measured Lyapunov exponent for various values of N is shown in Fig. 4. We find  that the large N behavior is given by13 The dependence on the temperature and the coupling follows from dimensional analysis, as explained in Sec. 2. The fact that the leading correction goes as N −2 is consistent with the 't Hooft counting.

Fast Scrambling
Quantum systems have a natural notion of 'scrambling time' t * , which is the time it takes a local perturbation to completely de-localize, or become 'scrambled'.In classical systems we can only discuss the scrambling time of a given perturbation (rather than as a property of the system itself).This is because we can make the growth time of a perturbation arbitrarily large by taking the initial perturbation to be small (in quantum systems we are limited by uncertainty).Earlier we defined the scrambling time to be the time at which |δX(t)| stops growing.We can then consider the N scaling of the scrambling time by scaling N while keeping the size of the initial perturbation fixed (quantum mechanically, the minimal size of a perturbation is O(N 0 )).
Let us now show that our classical system is a 'fast scrambler', namely one in which the scrambling time t * scales as log N 2 .The typical value of |δX(t)| when it stops growing can be estimated by picking two random configurations X and X from the ensemble and calculating the difference between them, |X − X | = Tr((X − X ) 2 ) ∼ √ N .We therefore expect the scrambling time t * to be given by We have already seen that λ L is independent of N to leading order.It is left to show that the perturbation indeed grows to be of order √ N .Fig. 5 shows the late-time evolution of |δX| for various N values.One can verify numerically that at late times |δX| ∼ √ N as expected to within less than one percent.This establishes fast scrambling in the classical matrix model.

Lyapunov Exponent from Poisson Brackets
The calculations described so far were classical in nature, relying on the time evolution of nearby points in phase space.On the other hand, computations of the Lyapunov exponent in quantum systems rely on commutators and out-of-time-order correlation functions [21,22].In this section we bridge this gap by extracting the exponent from the classical limit of commutators -Poisson brackets.The results agree with the ones obtained using the previous method.
To motivate the method from a classical perspective, consider a particle in D-dimensional space with spatial coordinates x I and momenta π I , where I = 1, . . ., D. One can use the classical correlator 14   {x to give an equivalent definition of the Lyapunov exponent λ L [21].Here we take I = J to ensure that the 1-point function vanishes. 15We expect that correlators of the form {V (t), W (0)} 2 p.b. (where V, W are operators that are local in time) exhibit the same exponential growth as (29).
In the matrix model we would like to focus on gauge-invariant correlators that have a good large N limit.We will first consider the correlator O ij (t, 0) 2 (with no summation over i, j), where 14 Classical correlators are defined by time-averaging (assuming ergodicity): Here Π i is the conjugate momentum to X i .We set i = j so that the one-point functions O ij (t, 0) vanish by SO(9) symmetry.The growth of the correlator is driven by the derivatives in (30), which are analogous to the derivative in (29).We therefore expect the correlator to grow as where λ L is the Lyapunov exponent of the matrix model.
Computing the correlator consists of the following steps.First, thermalize the system as before by evolving a random initial configuration for time t 0 = 4000 to obtain a reference configuration (X, V ).Next, define the perturbed configuration (X , V ) = (X + δX, V ) where δX i is a polynomial in V i with small, random coefficients.Given the reference configuration (X, V ) and the perturbed configuration (X + δX, V ), evolve both in time and compute O ij (t) 2 (the derivatives in (30) are approximated by replacing ∂X(t) → X (t) − X(t)).Finally, average the results over different choices of i = j (which are related by SO(9) symmetry), as well as over different initial values and choices of perturbation.
The resulting correlator ( 31) is shown in Fig. 6.An initial transient is followed by exponential growth, which saturates when the distance between the reference and perturbed configurations becomes of the same order as the system size.The fact that the growth stops at t 60 is an artifact of our approximation of the derivative in (30) using finite distances; the exact correlator keeps growing indefinitely.Fig. 7 shows the Lyapunov exponents we get by fitting the growing part of the curves. 16The large N behavior is given by 17 This is consistent with the previous result (26) obtained using Sprott's algorithm.
As mentioned above, we expect the Lyapunov exponent to not depend on the operators 16 For each sample of O ij (t, 0) 2 , we average the values at each t over i = j and then fit an exponent.The fitting window is between 10 −3 and 10 −11 times the saturation (late time) value of the correlator in the averaged sample.We then average the exponents that are obtained in this way from a few tens of samples with given N value.The error bars in Fig. 7 denote the statistical errors from this averaging of exponents. 17As in (26), the uncertainties quoted here do not take into account the error bars.
The 1-point function of O ijk vanishes for any choice of i, j, k.The result is shown in Fig. 6, and the Lyapunov exponent we obtain from this correlator is consistent with the previous results.

Lyapunov Spectrum Computation
In this section we go beyond the largest exponent and study the full spectrum of Lyapunov exponents [40], as defined in Sec. 3. The evolution of a perturbation δX i , δV i is given by the linearization of the equations of motion (10).Explicitly, where After discretization, the perturbation evolves according to where U (δt; x(t)) is the transfer matrix for a single time step.Our discretization (11) preserves the Gauss law constraint, and therefore the discretized transfer matrix should preserve the linearized constraint: 3. P gauge (x) project out pure gauge modes.The pure gauge modes at x = (X, V ) are spanned by the vectors By using an orthonormal basis of this space { w a }, the projector can be defined as We now define P (x) ≡ P gauge (x) • P Gauss (x) • P U (1) (x).It is easy to verify that P (x) is an orthogonal projector. 18The physical transfer matrix is then defined by (c.f. ( 57)) This physical transfer matrix has n = 2(d − 1)(N 2 − 1) nonzero singular values, and the n physical Lyapunov exponents can be computed from these by using (14).Fig. 8 shows the spectrum of the time-dependent exponents (15) for N = 6.Numerics19 limit us to this modest value of N .But the rapid convergence to the infinite N limit displayed above indicates that these results will be meaningful.Notice that the largest exponent is larger than the Lyapunov exponent λ (N =6) L 0.28, and that it decreases with time.In Fig. 9, the spectrum at a single time step t = δt is shown.The largest exponent is larger by an order of magnitude compared to t → ∞.What causes this suppression?Consider the singular vector v(t) of U phys (t; x 0 ) that corresponds to the maximal singular value.If v(t) stayed roughly constant, then the perturbation δx would quickly align itself with v(t), and the Lyapunov exponent would correspond to the maximal short-time exponent.Instead, our numerical results suggest that v(t) evolves quickly in time, such that the perturbation cannot become aligned with it.This suppresses the exponent over time, leading to a smaller λ L .At t 10, the spectrum is well described by the ansatz where λmax and γ both depend on time.Fig. 10 shows the finite-time positive Lyapunov spectrum for N = 6 and a fit to the ansatz (43).(Note that this λmax is a fitting parameter and is not exactly the same value as the largest Lyapunov exponent measured in the simulation.)We can see that λmax decreases with t (see also Fig. 8), while γ is consistently close to 0.5.More generally, we found that γ = 0.5 ± 0.1 in all checks we made.
There are two exponents at finite t which should both converge to the Lyapunov exponent as t → ∞: The largest exponent determined from the transfer matrix, which we call λ max (t), and the 'effective' exponent calculated in Sec. 4, which is defined by As shown in Sec.3.1, for generic perturbations we can approximate this exponent by Fig. 11 compares these exponents.It is surprising that λ L (t) quickly approaches its asymptotic value and then remains essentially constant, while λ max (t) converges much more slowly.We do not have an explanation for this behavior.It is consistent with the clean exponential growth of perturbations that we observed in Sec. 4. We tried several fitting ansatz to match the evolution of λ max (t) such as a + b t , a + b t + c t 2 , and a + b t c .It turned out a + b t + c log t t fits the data very well at a wide time window, and attains the correct late-time value λ max (t = ∞) 0.28 determined in Sec. 4. 20 By recalling the fact that γ stays close to 0.5, we can expect the late-time behavior to be where λ L 0.28 is the largest Lyapunov exponent determined in Sec. 4 and γ 0.5.The relation to random matrix theory has not escaped our attention, although we do not see how to make it precise.It would be interesting to see whether it holds for larger values of N as well.Fit of λ max :  0.87 t + 0.47 log(t)

Discussion
The above results show that the Lyapunov spectrum approaches a smooth N -independent limit as N → ∞ in this classical system.This is consistent with the observation of Sekino and Susskind [11] about the k-locality of the matrix quantum mechanics Hamiltonian in matrix index space.
In addition, these results bring into focus the existence of a whole spectrum of Lyapunov exponents in these large N systems.This raises the question about the nature of the corresponding spectrum in quantum large N systems and their meaning in the bulk gravity dual.There is one indication of this in existing work: In [17] it was pointed out that stringy corrections to Einstein gravity modify the rate of exponential growth of scrambling, tending to decrease it.In particular the result found there is Here k T is the transverse momentum of the perturbation in the field theory space and c 1 and c 2 are known positive constants.This gives a spectrum of Lyapunov exponents Sprott's algorithm.To compute the leading k exponents one chooses k initial orthonormal vectors, evolves each one according to Sprott's procedure, and also re-orthogonalizes the vectors after each time step.
indexed by k T .That chaos develops at different rates for different modes gives a physical interpretation of the diffusion of chaos from a localized perturbation found in [17].The analog of k T and its role in the spectrum in the setup described in this paper where there is not a field theoretic transverse space are interesting questions for future work.
Another important question is the quantum field theory and holographic origin of the order N 2 different Lyapunov exponents found in the classical system.Weakly coupled field theory analyses [17,25] should shed some light on the QFT question.We hope to return to this in subsequent work.
As discussed above, the existence of a spectrum of Lyapunov exponents immediately brings to mind Pesin's theorem and the Kolmogorov-Sinai entropy one gets by adding up the positive ones.This raises the important question of the meaning of the KS entropy in large N quantum systems and in their holographic duals.Suggestions for definitions of quantum KS entropy have been made in the literature [44][45][46][47][48][49].These should certainly be explored further given recent developments.We would expect any entropy, including the KS entropy, to be proportional to N 2 , sharpening the question mentioned above.
A simple model connecting KS entropy and entanglement entropy has been constructed in [50], motivated by [51].Building on these ideas the authors of [50] have recently constructed a toy field theory model relating entanglement entropy and KS entropy and conjectured a bound on the rate of growth of entanglement that follows from the bound on Lyapunov exponents [20].
The question about the holographic role of KS entropy and its relation to rates of increase of various entropies has been raised by a number of authors [52][53][54].Entanglement entropy growth is a natural candidate that has been discussed in detail [53].One hint of such a connection is an observation of Stanford [55].In Einstein gravity the butterfly velocity v B describing the transverse spread of chaos [15,18] is the same as the saturation velocity 22 in the entanglement tsunami picture of [56] that describes the rate of growth of the spatial region where entanglement entropy reaches its saturated value.This connection occurs because the Ryu-Takayanagi surface computing the entanglement entropy in this regime dips very close to the horizon, the region where the exponential blueshifts responsible for holographic chaos occur.
Let us show that λ x 0 (δx 0 ) = λx 0 (δx 0 ).Define r + > 0 by Here, T x is the tangent space at x.The maxima are well-defined because the norm is continuous, and both M and the unit ball at each point are compact.For any x ∈ M and any v ∈ T x , we then have The other inequality can be obtained using the definition This completes the proof.

B Lyapunov Exponents and Gauge Symmetry
In this section we construct a 'physical' phase space M phys for the matrix model by following the procedure of symplectic reduction (or Marsden-Weinstein reduction).The physical phase space is free of gauge redundancy.We then equip the space with a Riemannian metric, which allows us to define gauge-invariant Lyapunov exponents.
To construct the physical phase space we begin with the total space M, parameterized by (X, V ) and equipped with the symplectic form ω = dX i ab ∧ dV i ba .The dimension of M is 2d(N 2 − 1), where d = 9.As explained in Sec.3.2, gauge redundancy affects this space in two ways.First, Gauss's law restricts physical configurations to lie on the constrained surface Second, points on M 0 that are related by a residual gauge transformation (22) are physically identical.
We define the space M phys by identifying points on M 0 that are related by a gauge transformation.The physical phase space is the coset space M phys ≡ M 0 /∼, where for any x, x ∈ M 0 we say that x ∼ x if these points are related by a gauge transformation of the form (22). Points on M phys will be denoted by [x] where x ∈ M 0 .M phys will generally have a complicated global structure that includes singularities.However, a typical point on a given fixed-energy subspace has a smooth neighborhood, and we will only be interested in the local structure at such points.The dimension of M phys at such points is 2(d − 1)(N 2 − 1).
The tangent space at a point [x] ∈ M phys is obtained from the tangent space at x = (X, V ) ∈ M 0 by modding out infinitesimal gauge transformations.The subspace of gauge transformations at (X, V ) is spanned by the vectors (δX H , δV H ) where and H is any traceless, Hermitian N × N matrix.Vectors on the physical tangent space, which we denote by [δx], obey [δx] = [δx + (δX H , δV H )] for any Hermitian H.
In order to regard M phys as a phase space, we must equip it with a symplectic structure.Wherever M phys is smooth, we define the natural symplectic form Here [δx], [δx ] are vectors at [x] ∈ M phys , and δx, δx are chosen representative vectors at x.It is easy to verify that this definition is independent of the choice of representatives because of the Gauss law constraint.
To define Lyapunov exponents we must also equip M phys with a Riemannian metric.Let us use the metric g on M (c.f. ( 21)) to define a metric g phys on the physical phase space.First, restrict g to M 0 and consider the tangent space at a point x = (X, V ) ∈ M 0 .Let P gauge denote the orthogonal projection operator (with respect to g) that projects out the pure gauge vectors.We now define the inner product of two vectors [δx] = [(δX, δV )], [δx ] = [(δX , δV )] on M phys by g phys ([δx], [δx ]) ≡ g(P gauge (δX, δV ); P gauge (δX , δV )) .
On the right-hand side we have chosen representative vectors (δX, δV ), (δX , δV ) at x.The metric is well-defined, in that it is independent of the choice of vector representatives and of x.In particular, notice that the problem that prompted the introduction of the physical metric is now solved: Two points on M 0 that are related by a gauge transformation are equivalent on the reduced M phys , and have vanishing distance under the physical metric.

B.1 Gauge-Invariant Exponents
Lyapunov exponents can now be defined for fixed energy subspaces of M phys using the physical metric, and they will be independent of our choice of metric as shown in Appendix A. The first step is to define a transfer matrix U phys that only propagates physical modes.It can be done by a projection where P (x) ≡ P gauge (x)P Gauss (x)P U (1) (x) is an orthogonal projector defined in Sec. 5. Given a generic initial vector on M, P Gauss (x 0 ) restricts the perturbation to lie on M 0 , and P gauge (x 0 ) removes the pure gauge modes.This chooses a representative vector on M phys .The vector then propagates with the usual transfer matrix U (t; x 0 ).After propagation we project again.
To compute the Lyapunov exponents, perform the singular value decomposition of U phys .There are 2d(N 2 −1) singular values, of which 2(N 2 −1) vanish due to the projections.The gauge-invariant Lyapunov exponents are computed from the remaining (positive) singular values by using ( 14).
As we now show, the physical transfer matrix U phys is symplectic with respect to ω phys .As a result, the physical Lyapunov exponents are paired.To show that U phys is symplectic, we need to show it obeys the equation Here we introduced the notation U phys (t; x) ≡ U phys (x → x(t)) for clarity.ω phys (x) and ω phys (x ) are matrix representations of ω phys using the same bases we use to represent U phys .They are given by ω phys (x) = P (x) • ω(x) • P (x), where ω(x) represents the symplectic form on the total phase space (we may choose ω(x) to be constant, but this is not necessary).
Notice that the matrix ω phys (x) generally depends on x.Equation (58) can be written as Now we claim that the P (x ) factors on the left-hand side re redundant.To see this, first note that P (x )U (x → x )P (x) = P gauge (x )U (x → x )P (x), due to the fact that time evolution preserves the Gauss law constraint.Further, the remaining factor of P gauge can be dropped because, after reducing to the Gauss-constrained subspace, pure gauge perturbations vanish automatically in the symplectic form.We therefore are left with the equation which follows immediately from the fact that U is symplectic with respect to ω.This concludes the proof that U phys is symplectic on the physical space, and therefore the physical Lyapunov exponents are paired.
C Perturbation Compatible with the Gauss Law Constraint Given a thermalized configuration X(t), V (t), we would like to perturb it slightly while preserving the Gauss law constraint (5).We will do this by deforming the potential energy with additional interaction terms that preserve the constraint, evolving the system for a short time to obtain a perturbed configuration X (t), V (t), and then restoring the original Lagrangian.We add the following term to the potential, The force is modified from F i (t) to where {• , •} is the anti-commutator.The Gauss law constraint is still preserved, because the modified force still satisfies the relation i [X i (t), F i (t)] = 0.In practice we choose k 0 = 2, the coefficients c k are chosen randomly from N (0, 10 −8 ), and we evolve the deformed system for time t 1 = 1 before turning off the deformation. 23

D Sprott's Algorithm
In this section we describe Sprott's algorithm [43], which we use in Sec. 4 to compute the leading Lyapunov exponent of the classical matrix model.The basic idea behind the algorithm is to rescale the perturbation at each time step such that it stays small, the linear approximation continues to hold, and the growth never saturates.The evolution can then continue until the measured exponent converges with some chosen precision.
4. Define the new configuration x n by The difference has been rescaled such that x n and x n are again at distance d 0 .
5. Repeat steps 2-4.The leading Lyapunov exponent is given by Note that the rescaling in step 4 implies that the new configuration x n does not satisfy the Gauss law constraint.However, the violation is subleading in the size of the perturbation, and we verified numerically that the violation remains negligible over the entire time evolution.

E Finite Volume of Classical Flat Directions
Consider the U (N ) theory with d matrices (d = 9 in the case considered in this paper).First let us consider the simplest situation, where one of the D0-branes is separated by a distance L from the other (N − 1) branes that are close to each other.By using the residual gauge symmetry and SO(d) symmetry, we can take X d N N L and have all other matrix components be much smaller than L. Then the potential energy coming from N -th row and column is approximately 1 aN | 2 (here we are neglecting contributions from elements that do not scale with L).This contribution must be smaller than the total energy E (which is kept fixed), and therefore possible values of X i aN are suppressed as L becomes large, as d−1 i=1 N −1 a=1 |X i aN | 2 g 2 E/L 2 .Hence the phase space volume for L > L 0 is suppressed by at least The factor L −2(d−1)(N −1) comes from the integral with respect to (d − 1)(N − 1) complex variables X i aN (a = 1, . . ., N − 1) and L d−1 comes from SO(d) rotational symmetry.As L 0 goes to infinity, the volume vanishes except for when d = 2, N = 2.In other words, this flat direction occupies a finite volume in phase space unless d = 2, N = 2.
When more eigenvalues are separated from the rest, more off-diagonal elements have to become small and hence such configurations have even smaller measures.

Figure 1 :
Figure 1: Lyapunov exponent as a function of thermalization time t 0 .

Figure 2 :
Figure 2: Time evolution of |δX(t)| for N = 16.Here t = 0 is the time of the initial perturbation.

Figure 3 : 4 eff T = 1 .Figure 4 :
Figure 3: The exponent estimated using the Sprott algorithm as a function of time, for N = 20 and at λ 1/4 eff T = 1.The band represents the statistical fluctuations of different samples.

Figure 8 :
Figure 8: Distribution of Lyapunov exponents with N = 6 at different times, both normalized with unit area.

Figure 10 :
Figure 10: (a) Positive Lyapunov spectrum for N = 6 and a fit to the ansatz (43) at the largest t we have studied, and (b) the fitting parameters γ versus t.Here we normalize the area of only the positive part of the spectrum to unity, and multiply the right-hand side of (43) by 2 accordingly.

Figure 11 :
Figure 11: The largest and average exponents as a function of time, with N = 6.Data is an average over 80 samples.