Universality in Chaos: Lyapunov Spectrum and Random Matrix Theory

TL;DR

This paper reveals a new universality in large chaotic systems where the Lyapunov spectrum's statistical properties align with Random Matrix Theory, demonstrated through models related to black holes and matrix products.

Contribution

It introduces the concept of universality in Lyapunov spectra for large chaotic systems and connects it to Random Matrix Theory, with evidence from stringy black hole models and random matrix products.

Findings

Lyapunov spectra follow RMT statistics in large chaotic systems.

Universal behavior appears immediately in the massless limit.

Universal patterns emerge at late times in other models.

Abstract

We propose the existence of a new universality in classical chaotic systems when the number of degrees of freedom is large: the statistical property of the Lyapunov spectrum is described by Random Matrix Theory. We demonstrate it by studying the finite-time Lyapunov exponents of the matrix model of a stringy black hole and the mass deformed models. The massless limit, which has a dual string theory interpretation, is special in that the universal behavior can be seen already at t=0, while in other cases it sets in at late time. The same pattern is demonstrated also in the product of random matrices.

Figures

Equations

$$L=\frac{N}{2}\mathrm{Tr}\left(\sum _I(D_t{X}_I{)}^2+\frac{1}{2}\sum _{I\ne J}[X_I,X_J{]}^2\right),$$

$$\frac{d^2{X}_I}{d{t}^2}=\sum _J[X_J,[X_I,X_J]],$$

$$\sum _I\left[\frac{d{X}_I}{dt},X_I\right]=0.$$

$$\rho (\lambda ,t)=\frac{3}{4{\tilde{\lambda }}_{\max }^{3/2}}\sqrt{{\tilde{\lambda }}_{\max }-|\lambda |},$$

$$Z(\tau )=\sum _n{e}^{i{\lambda }_n\tau }$$

$$Z_{\mathrm{GOE}(K)}(\tau )=\sum _n{e}^{i{E}_n\tau }.$$

$$\mathcal{M}(t)=M_t{M}_{t-1}\cdots M_2{M}_1.$$

$$⟨r⟩={⟨\frac{\min (s_i,s_{i+1})}{\max (s_i,s_{i+1})}⟩}_i,$$

$$\Delta (K,t)$$

$$\simeq \sum _{q=0}^{q_{\max }}|P(s_q)-P_{\mathrm{GOE},q}|\Delta s,$$

$$P_{\mathrm{GOE}(\mathrm{Wigner})}(s)=\frac{\pi s}{2}{e}^{-\frac{\pi }{4}{s}^2}.$$

$$\delta P_q\equiv \sqrt{(L-1)\left(\frac{1}{L}\sum _{l=1}^L{\left(P_q^{(l)}\right)}^2-P_q^2\right)}.$$

$$\Delta (K,t)\pm {\delta }^{(\pm )}\left(\Delta (K,t)\right)=\sum _q{\delta }^{(\pm )}{\left(\Delta (K,t)\right)}_q\Delta s,$$

$${\delta }^{(+)}{\left(\Delta (K,t)\right)}_q$$

$$=$$

$${\delta }^{(-)}{\left(\Delta (K,t)\right)}_q$$

$$=$$

Full text

††preprint: YITP-17-17, KEK-TH 1957

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

[TABLE]

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

[TABLE]

supplemented with the Gauss’s law constraint

[TABLE]

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

[TABLE]

where $\tilde{\lambda}_{\rm max}$ is a time-dependent parameter which approximately equals to the largest Lyapunov exponent. Near the edge $|\lambda|\sim\tilde{\lambda}_{\rm max}$, this distribution is equivalent to the semi-circle, $\sqrt{\tilde{\lambda}_{\rm max}^{2}-\lambda^{2}}$. This is an indication of a possible connection to RMT.

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.

III.1 Mass deformation

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

[TABLE]

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

[TABLE]

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.

In Fig. 4 we have plotted $g(\tau)\equiv|Z(\tau)|^{2}/K^{2}$ calculated from the Lyapunov spectrum of the BFSS matrix model at $t=0$ and $g_{\mathrm{GOE}(K)}(\tau)\equiv|Z_{\mathrm{GOE}(K)}(\tau)|^{2}/K^{2}$. The agreement at large $\tau$ (the ramp $\sim\tau^{1}$ and the plateau $\sim\tau^{0}$) means the agreement of the Lyapunov spectrum and RMT energy spectrum beyond the nearest neighbor. Note that the disagreement in the small $\tau$ region is not a problem, it simply means the global shapes of the spectrum are different.

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),

[TABLE]

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

[TABLE]

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.

Furthermore, in order to see the correlation over even larger separations, in Fig. 8 (a) we have compared the SFFs for the product of real matrices, $|Z(\tau)|^{2}/|Z(\tau=0)|^{2}$, with that of GOE random matrices, $|Z_{\rm GOE}(\tau)|^{2}/|Z_{\rm GOE}(\tau=0)|^{2}$. We can see that $|Z(\tau)|^{2}$ approaches to $|Z_{\rm GOE}(\tau)|^{2}$ as $t$ increases. Also in Fig. 8 (b) we have plotted $g(\tau)$ for complex random matrix products against $g_{\mathrm{GUE}}(\tau)$ obtained from GUE random matrices. Here again, we can see the agreement between the finite-time Lyapunov exponents and RMT energy spectrum beyond the nearest neighbors.

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.

Supplementary materials

V.1 Details of the analysis of the unfolded spectrum:

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.

For Hamiltonian systems discussed in this paper, all exponents are paired with the exponent of the same absolute value and the opposite sign. Therefore we focus on positive Lyapunov exponents. We further omit both largest $5\%$ and smallest $5\%$ of the positive exponents, in order to avoid the exponents close to the edge affecting the fit discussed below. We denote the maximum and minimum of retained exponents by $\lambda^{(\mathrm{max})},\lambda^{(\mathrm{min})}$ respectively. For the bins containing retained exponents we fit the density of exponents $\rho(\lambda)$, by a polynomial $\tilde{\rho}(\lambda)=\sum_{k=0}^{k_{\mathrm{max}}}a_{k}(\lambda-\lambda_{0})^{k}$ of $\lambda$, for unfolding the spectrum. We typically choose $k_{\mathrm{max}}=10$. To reduce numerical error, $\lambda_{0}$ is chosen within the fitting range $[\lambda^{(\mathrm{min})},\lambda^{(\mathrm{max})}]$.

Then the spectrum is ‘unfolded’ by considering $s_{j}^{(w)}\equiv S(\tilde{R}(\lambda_{j}^{(w)})-\tilde{R}(\lambda_{j+1}^{(w)}))$, in which $\tilde{R}(\lambda)=\int_{\lambda_{0}}^{\lambda}\tilde{\rho}(\lambda^{\prime})d\lambda^{\prime}=\sum_{k=0}^{k_{\mathrm{max}}}\frac{a_{k}}{k+1}(\lambda-\lambda_{0})^{k+1}$ and $S\sim K$ is the normalizing factor chosen so that the average of $s_{j}^{(w)}$ is unity.

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

[TABLE]

in which we have defined $P_{\mathrm{GOE},q}\equiv P_{\mathrm{GOE}}(s_{q})$.

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

[TABLE]

However, for our purpose the Wigner’s surmise is not accurate enough. The correct distribution $P_{\rm GOE}(s)$ admits a Taylor series expansion and a Padé approximant, which are available in Dietz:1990 . In our analysis, it is sufficient to use the Taylor series expansion of $P_{\rm GOE}(s)$ as its approximation for $s\leq 3$. We use the upper limit, $s_{q_{\mathrm{max}}}\simeq 3$, in the summation (9).

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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.