Universality in Chaos: Lyapunov Spectrum and Random Matrix Theory

I Introduction and Summary

In this paper we suggest that the statistical property of the Lyapunov spectrum in classical chaotic systems with a large number of degrees of freedom is described universally by Random Matrix Theory (RMT). More precisely, we consider the spectrum of the finite-time Lyapunov exponents, which is defined from the growth of small perturbations during a finite time interval $t$. Unlike the majority of the previous references in which $t\to\infty$ is taken first, we will take the limit of large number of degrees of freedom at each finite $t$ footnote_limit . This is a natural limit which leads to various universal results such as the universal bound on the Lyapunov exponent Maldacena:2015waa .

Our initial motivation was in a different kind of universality in quantum many-body chaos, which has been a hot topic in string theory and quantum information communities in recent years (see e.g. Sekino:2008he ; Maldacena:2015waa ). It has been argued that the largest Lyapunov exponent $\lambda_{\rm max}$ has to satisfy a certain bound, and the black hole in general relativity saturates the bound Maldacena:2015waa . In this context G. Gur-Ari, S. Shenker and one of the authors (M. H.) have studied Gur-Ari:2015rcq the Lyapunov exponents of a classical matrix model (the D0-brane matrix model) deWit:1988wri ; Witten:1995im ; Banks:1996vh ; Itzhaki:1998dd which is related to a quantum black hole with stringy corrections via the gauge/gravity duality Maldacena:1997re ; Itzhaki:1998dd . They found that the global distribution of the Lyapunov exponents follows the semi-circle law near the edge, which is a characteristic feature of the energy spectrum of RMT. This suggested the existence of certain universal behaviors in the Lyapunov spectrum of such systems.

Motivated by this observation, we studied the statistical property of the Lyapunov spectrum in the matrix model footnote_reference . As we will show, its statistical property is described by RMT for all $t$. When we introduce the mass deformation, the RMT description is lost for small $t$. However, it does emerge for large $t$. The spectrum of the product of random matrices, which has been studied as an analytically tractable model of chaos, admits the same RMT description. This is true in other models as well; some examples will be reported in HST_to_appear . Based on these results, we conjecture that the Lyapunov exponents of a large class of many-body chaos, both deterministic and nondeterministic, are described by RMT at late time.

II Lyapunov exponent and Lyapunov spectrum

Let us consider the phase space consisting of $K$ variables, $\phi_{i}$ ($i=1,2,\cdots,K$). By solving the equations of motion, the classical trajectory $\phi_{i}(t)$ is obtained depending on the initial condition at $t=0$. When a small perturbation is added at $t=0$, $\phi_{i}\to\phi_{i}+\delta\phi_{i}$, the time evolution of the perturbation can be evaluated by solving the equations of motions with the perturbed initial condition. When $\delta\phi_{i}$ is infinitesimally small, the evolution is described by the transfer matrix $T_{ij}(t,t^{\prime})$ ($t>t^{\prime}$) as $\delta\phi_{i}(t)=\sum_{j}T_{ij}(t,t^{\prime})\delta\phi_{j}(t^{\prime})$. Let $a_{1}(t,t^{\prime})\geq a_{2}(t,t^{\prime})\geq\cdots\geq a_{K}(t,t^{\prime})>0$ be the singular values of $T_{ij}(t,t^{\prime})$. The time-dependent Lyapunov exponent $\lambda_{i}(t,t^{\prime})$ is defined by $\lambda_{i}(t,t^{\prime})=\frac{\log a_{i}(t,t^{\prime})}{t-t^{\prime}}$.

When the trajectory is bounded, the exponents have unique limits $\lim_{t-t^{\prime}\to\infty}\lambda_{i}(t,t^{\prime})$. Usually they are called the Lyapunov exponents. An existence of a positive exponent characterizes the sensitivity to the initial condition, which is a necessary condition for the chaos.

In this paper we consider the finite-time exponents, and study their statistical properties at large $K$. Note that we take the large-$K$ limit for each fixed time interval $t-t^{\prime}$, and use many samples which are generated from different initial conditions. Two limits, $K\to\infty$ and $t-t^{\prime}\to\infty$, may or may not commute, depending on the systems footnote_limit . In chaotic systems, generic initial states evolve to ‘typical’ states after some time, and the statistics is dominated by them. We will pick up only typical states. It can be achieved by taking $t$ to be sufficiently late time. For the simplicity of the notation, we will redefine the time and set $t^{\prime}=0$, and call $\lambda_{i}(t,0)$ as $\lambda_{i}(t)$.

In order to compare the statistical property of the Lyapunov spectrum with RMT, we use the standard unfolding method Brody:1981cx . Note that $\{\lambda_{i}(t)\}$ and $\{a_{i}(t)\}$ lead to the same unfolded distribution. Hence the universality of the Lyapunov exponents discussed in this paper is equivalent to the universality in the singular values of the transfer matrix describing the linear response.

III D0-brane matrix model

In Gur-Ari:2015rcq , the classical limit of the matrix model of D0-branes has been considered footnote_BFSS_reference . The Lagrangian is given by

where $X_{I}$ $(I=1,\dots,d)$ are $N\times N$ traceless Hermitian matrices; $D_{t}X_{I}=\partial_{t}X_{I}-[A_{t},X_{I}]$, where $A_{t}$ is the $SU(N)$ gauge field. The number of the traceless Hermitian matrices is $d=9$. This system has a scaling symmetry which relates solutions with different energies. We will employ a natural energy scale $E=6(N^{2}-1)-27$ footnote_normalization , which corresponds to the unit temperature, $k_{\rm B}T=1$. We use the same simulation code as in Gur-Ari:2015rcq .

In the $A_{t}=0$ gauge, the equation of motion is

supplemented with the Gauss’s law constraint

By following the procedures explained in Gur-Ari:2015rcq , we can study the Lyapunov exponents. In Gur-Ari:2015rcq , it has been observed that the spectrum of $\lambda$ is well approximated by

We have studied the Lyapunov spectrum for $0\leq t\leq 10$ with $N=4,6,8$. The number of the Lyapunov exponents, which appear in pairs of positive and negative ones with the same absolute value, is $K=16(N^{2}-1)$ footnote_DOF . We ordered the positive exponents as $\lambda_{1}\geq\lambda_{2}\geq\cdots$, and studied the distribution of the level spacing $s_{i}\equiv\lambda_{i}-\lambda_{i+1}$. From these exponents, the distribution $P(s)$ of the unfolded level separation can be obtained. (For the detail of the analysis, including the error estimate, see the supplementary materials.) It agrees well with the nearest-neighbor level statistics of the GOE ensemble, which we denote by $P_{\rm GOE}(s)$ Dietz:1990 , as shown in Fig. 1, for all values of $t$. Already at $t=0$, the spectrum agrees very well with GOE; see Fig. 1 (a). Note that we can see a small deviation from GOE at $N=4$. Thus the data strongly suggest that the level statistics of the finite-time Lyapunov spectrum agrees with that of GOE at any $t$, after taking the large-$N$ limit.

Next we add the mass term $\Delta L=-\frac{Nm^{2}}{4}{\rm Tr}\sum_{I}X_{I}^{2}$ to the D0-brane matrix model. The physically meaningful parameter is the dimensionless ratio $E/m$. Here we fix the energy to be $E=6(N^{2}-1)-27$ and change $m$. In the limit with an infinite mass, or equivalently the zero-energy limit, the theory becomes a free theory, which is not chaotic footnote_mass .

In Fig. 2 (a) the distribution of the unfolded level separations with $m=3$ is shown. Although it is linear in $s$ for small $s$, indicating level repulsion between Lyapunov exponents, the distribution disagrees with that of GOE, having a peak at smaller $s$ and a longer tail. However, as shown in Fig. 2 (b), the distribution goes close to GOE at $t>0$.

To make this observation more precise we calculated the difference, $\int ds|P(\lambda)-P_{\rm GOE}(\lambda)|$, of the distribution from that of GOE. The difference is plotted at $t=0$ for several values of $m$ in Fig. 3 (a). The spectrum disagrees with that of GOE at finite $m$, and the deviation is larger when $m$ is larger. In Fig. 3 (b), the time dependence is shown for $m=3$, $N=4,6,8$. The deviation from $P_{\rm GOE}(s)$ oscillates, and gradually decreases. This result strongly suggests that the distribution converges to $P_{\rm GOE}(s)$ when the limit $t\to\infty$ is taken after $N\to\infty$.

III.2 Beyond nearest neighbor

In order to see the agreement with RMT beyond the nearest-neighbor level correlation, we have compared the spectral form factor (SFF) defined by

and its RMT counterpart for Gaussian symmetric random matrices of the same dimension $K$,

The spectral form factor captures more information about the spectrum, the so-called spectral rigidity. The large $\tau$ behavior of the SFF reflects the fine grained structure of the energy spectrum. The small $\tau$ region is sensitive to the global shape of the spectrum, which is not expected to be universal.

We repeated the same analysis with a mass deformation. In Fig. 5, the SFFs $g(\tau)$ for the mass-deformed model with $N=8$ and $m=3$ for $t=1$ and $t=10$ are shown. The convergence to RMT at late time (large $t$) can be seen very clearly.

IV Product of random matrices

Let us consider a product of $t$ matrices randomly chosen from a certain ensemble (‘Random Matrix Product’, RMP),

We take the matrix size to be $K\times K$. The RMP has been studied as a toy model of the Lyapunov growth, by regarding $M_{i}$ to be an analogue of the transfer matrix at a short time separation. From the singular values $a_{i}(t)(i=1,2,\cdots,K)$, ordered as $a_{1}(t)\geq a_{2}(t)\geq\cdots\geq a_{K}(t)$, we define the finite-time Lyapunov exponents by $\lambda_{i}(t)=(\log a_{i}(t))/t$.

The RMP has also been considered in the study of quantum transport phenomena, such as the conduction of electrons in a disordered wire RefQT . Our analysis in this section is closely related to results in the literature of the quantum transport phenomena; our $K$ corresponds to the number of transport channels, and $t$ corresponds to the length of the disordered wire footnote_QTvsGOE . In quantum transport phenomena, the evolution is studied of the transmission eigenvalues when the length of the wire is changed EvolutionQuantumTransport . It would be interesting to consider the time evolution of Lyapunov spectrums of the classical (deterministic or non-deterministic) chaotic systems from a similar point of view.

If each $M_{i}$ is a real matrix (also a complex matrix) with the weight $e^{-K{\rm Tr}MM^{\dagger}}$, then the level spacing statics of Lyapunov exponents $\lambda_{i}(t)$ follow that of the standard GOE (GUE) for any fixed $t$. This is easily verified numerically, and for the complex matrices an analytic derivation can be found in ProductRandomMatrices . This is precisely analogous with the case of the massless D0-brane matrix model (1). Note that $t\to\infty$ with fixed $K$ is different from RMT Newman1986 footnote_larget .

One can also introduce a deformation of the RMP playing a role analogous to the mass deformation of the matrix model. We have numerically studied a product of real-valued random band matrices, whose $(i,j)$ components are set to zero unless $|i-j|<h$, with the periodic identification $i\sim i+K$. As shown in Fig. 6 (a), the deviation of $P(s)$ from GOE at $t=1$ converges to an $O(K^{0})$ value in the large-$K$ limit when $h/\sqrt{K}$ is fixed. In Fig. 6 (b), the results for the products with $h/\sqrt{K}=1/2$ are shown. At large $t$, the plot shows a clear tendency of the convergence to GOE.

We also calculate the average nearest neighbor gap, defined by

in which $s_{i}=\lambda_{i}-\lambda_{i+1}$ and the average $\langle\cdots\rangle$ is taken over $i=1,\ldots,K-2$ and all the samples. The average nearest neighbor gap characterizes the correlation between the neighboring gaps in the spectrum. In Fig. 7 we have plotted the value of $\langle r\rangle$, both for products of real and complex matrices, against the inverse of the number of multiplied matrices $t$, both for complex and real matrices with $K=900$ and $h=16,13,10$, along with the values for GOE and GUE matrices presented in Atas2013 . This is the evidence that the universality holds for next-to-next nearest neighboring levels.

V Discussions

In this paper we have suggested the existence of a new universality in the Lyapunov spectrum of the classical chaotic systems based on numerical evidence for the matrix models and random matrix products. The massless D0-brane matrix model and the product of un-banded Gaussian random matrices are special in that the universal behavior can be seen at any time scale. It is interesting to speculate that other Yang-Mills theories and/or quantum gravitational systems satisfy the same property. Classical field theory calculations which are useful for this direction can be found in e.g. Bolte:1999th ; Kunihiro:2010tg .

We have also studied several other systems, e.g. 3d Coulomb gas, coupled Lorenz attractors and coupled logistic maps, and observed qualitative evidence for the same universality HST_to_appear . In general, the scaling of $t$ and the number of degrees of freedom should be carefully studied. For example, although the random matrix product with fixed $h$ and fixed $t$ does not become RMT, it is likely that $h$ fixed and $t\sim K^{p}$, with a certain power $p>0$, can lead to RMT.

A possible path toward an understanding of the mechanism behind the universality is to see how the spectra of various systems converge to RMT. As we commented in section IV, the classical chaotic systems and quantum transport phenomena are mathematically closely related, and thus it may be possible to deepen understanding of existence of universalities by considering both phenomena together. It may also provide us with a new characterization of various chaotic systems; the amount of deviation from RMT may be reflecting the strength of chaos, and the special property in the D0-brane matrix model would be related to the fast scrambling Sekino:2008he ; Maldacena:2015waa . The generalization of this universality to the quantum chaos would be even more interesting. We hope that the study of the statistical properties of the Lyapunov exponents provides us with a new viewpoint for studying chaotic systems.

Acknowledgement: We would like to thank S. Aoki, P. Buividovich, P. Damgaard, E. Dyer, A. M. García-García, G. Gur-Ari, S. Hikami, J. Magan, S. Nishigaki, S. Sasa, A. Schäfer, S. Shenker, A. Streicher, K. Takeuchi, A. Ueda, P. Vranas and M. Walter for discussions.

This work was partially supported by JSPS KAKENHI Grant Numbers JP25287046 (M.H.), JP17K14285 (M.H.), JP15H05855 (M.T.), JP26870284 (M.T.), JP17K17822 (M.T.) and JP16H06490 (H.S.). Part of computation in this work was performed at Supercomputer Center, Institute for Solid State Physics, University of Tokyo.

References

Supplementary materials

We explain how we produced the plots in this paper. We take $W$ independent samples labelled by $w=1,2,\ldots,W$. Each sample consists of $K$ Lyapunov exponents $\lambda_{1}^{(w)}\geq\lambda_{2}^{(w)}\geq\ldots\geq\lambda_{K}^{(w)}$.

We first make a histogram with bins of width $\Delta\lambda$ using all $W$ samples. There are $WK$ exponents in total. We then normalize the histogram so that $\int\rho(\lambda)d\lambda=\sum_{i}\rho_{i}\Delta\lambda=1$, where $i$ is a label for the bins. For $\mathcal{O}(10^{7})$ exponents we use in the majority of our plots, we typically take $\mathcal{O}(10^{3})$ bins.

We plot the histogram of $s_{j}^{(w)}$. Namely, for each bin $[q\Delta s,(q+1)\Delta s)$, we count the number $n_{q}$ of $s_{j}^{(w)}$ within this bin, and take $P(s_{q}\equiv(q+\frac{1}{2})\Delta s)=n_{q}/(\Delta s\sum_{q}n_{q})$.

From the distribution $P(K,t)$ with given $(K,t)$, we define the deviation from the GOE distribution by

When the average separation is normalized to be 1, the GOE distribution is often approximated by Wigner’s surmise,

V.2 Error estimate

Firstly we separate the samples to $L$ groups. We used $L=4$. We prepare $L$ data sets, by excluding one of the $L$ groups. By using a certain bin size, we make a histogram for each data set, and determine the heights $P^{(l)}_{q}$, where $l=1,2,\cdots,L$ is the label for the data set, and $q$ is the label for the bin. The Jack-knife error is defined by

This error estimate is used for the error-bars in figures 1 and 2.

Let $P_{q}^{\rm max}\equiv P_{q}+\delta P_{q}$ and $P_{q}^{\rm min}\equiv P_{q}-\delta P_{q}$. We denote the bin width by $\epsilon$. We estimate the error-bar for $\Delta(K,t)$, which we denote by $\delta^{(\pm)}\left(\Delta(K,t)\right)$, as

and $\delta^{(-)}\left(\Delta(K,t)\right)_{q}=0$ if $P_{i}$ and $P_{\rm GOE}$ coincides within the error estimate explained above (i. e. if $P_{q}^{\rm min}\leq P_{{\rm GOE}.q}\leq P_{q}^{\rm max}$), otherwise

V.3 The Lyapunov spectrum for the D0-brane matrix model

In Figures 9 and 10 we plot the Lyapunov spectrum obtained for the D0-brane matrix model at $t=0$ and $t=10$, respectively. The plots are symmetric about $\lambda=0$, therefore we have plotted only the positive exponents. The data suggest that $\rho(\lambda)$ rapidly approaches the large-$N$ limit.