From Integrable to Chaotic Systems: Universal Local Statistics of Lyapunov exponents

Introduction

Random matrices have a long tradition in describing universal aspects of spectral statistics, typically of static Hamiltonians in quantum physics . Today’s applications of random matrix theory (RMT) cover all areas of physics and include other sciences, cf. for a recent compilation. The random matrix description is usually based on global symmetries and not on dynamical principles. However, already in early studies Dyson analysed dynamical aspects of Brownian motion (BM) in matrix space, especially the dynamics which are formed by sums of matrices. This work initiated a discussion also on random matrix dynamics, stochastic processes in matrix space and the underlying dynamical principles.

A different and a priori unrelated idea to include dynamics is to study systems where the time dependence is multiplicative. This is the main objective of the present work. Such processes can be found in many physical systems e.g. in the form of transfer matrices, leading to the Dorokhov–Mello–Pereyra–Kumar (DMPK) equation , through propagation kernels , or more recently in quantum entanglement . Long time ago Furstenberg and Kesten have proposed to multiply $M$ random matrices as a toy model for chaotic dynamical systems. The corresponding Lyapunov exponents become deterministic in the large-$M$ limit , and we refer to the reviews for applications in disordered systems and to for the mathematical statements. It was conjectured early in that the Lyapunov spectrum and the universal statistics of a single random matrix are related. The claim was substantiated in a stability analysis of complex systems using stochastic processes . Different recent applications of such multiplicative processes include non-perturbative quantum gravity with space-time foliation , or the propagation of information in telecommunication . Multiplying $M$ random $N\times N$ matrices applies further to multi-layered complex networks with $M$ layers and $N$ degrees of freedom per layer , to the complexity of random maps , and in computer science (cf. ) through machine learning . All these examples can be viewed as certain realisations of progressive scattering, where a distinguished (time) direction exists. Therefore, the product matrix can be often interpreted as a transfer matrix. We will stick to this interpretation in our discussion.

How can one relate time dependent processes and universal equilibrium statistics of RMT? The standard dichotomy considered in RMT for the spacing distribution of consecutive energy levels of a static Hamiltonian is between Poisson statistics characterising integrable systems - the Berry–Tabor conjecture - and the (approximate) Wigner surmise in RMT characterising chaotic quantum systems - the Bohigas–Giannoni–Schmit conjecture . Such transitions have been an intense object of study .

Fixing all quantum numbers but one in an integrable Hamiltonian quantum system, like the hydrogen atom or harmonic oscillator, the respective spectrum can be unfolded to become equidistant. It is the incommensurable superposition of several such picket fence spectra that then becomes Poisson . In contrast, our results only exhibit a single picket fence structure, realised by the individual Lyapunov exponents, in the limit where they become deterministic. When increasing the system size compared to the number of time steps, the Lyapunov exponents start to interact with each other. We will find that their statistics can then be described by Dyson’s BM with fixed initial conditions.

Rather than focussing on the level spacing distribution which is a classical measure for random matrix statistics, we focus on the correlation kernel, encoding all correlation functions in the underlying determinantal point process . For recent work on the spacing of Lyapunov exponents, cf. .

Let us review some recent developments in products of random matrices, cf. . Multiplying $M$ independent complex Ginibre matrices having $N^{2}$ complex normal elements, exact analytical expressions for the joint probability distribution functions of the eigenvalues and singular values have been derived for finite $M$ and $N$. They are given by a determinantal point process which means that all $k$-point correlation functions can be written as determinants,

of a correlation kernel $K(x_{j},x_{l})$. The kernel is expressed in terms of Meijer G-functions . The level density $R_{1}(x)=K(x,x)$ is the simplest example. These findings have initiated considerable activity for finite $M$ and $N\rightarrow\infty$, see e.g. and references in the review for more general products.

Based on these findings the full statistics for finite Lyapunov exponents was derived in at fixed $M$ and $N$. Therein, it was argued that for the local statistics a non-trivial double scaling limit $M,N\rightarrow\infty$ should exist, leading away from the deterministic limit $M\to\infty$, for $N$ fixed (or $N\to\infty$) . In contrast, for fixed $M$ the universal local statistics of a single Gaussian Unitary Ensemble (GUE) was found in the bulk and at the soft edge of the spectrum of Lyapunov exponents . A similar transition from deterministic to GUE statistics was also seen in for the solution of the DMPK equation, where the authors pointed out a connection to Dyson’s BM. Let us emphasise that the matrix products taken for DMPK are perturbations of the unit matrix while in the present work we consider more general products of complex matrices, for which this transition has not been studied before, and point out the general mechanism behind this transition.

Qualitative Discussion

We consider $M+1$ time steps (layers) of the system, where each step is described by an $N$-dimensional vector $\vec{v}_{j}$, $j=0,\ldots,M$. In the simplest situation $\vec{v}_{j}=X_{j}\vec{v}_{j-1}$, the evolution from an initial condition $\vec{v}_{0}$ to time $M$ follows the product matrix $Y=X_{M}\cdots X_{1}$. We call $Y$ the transfer matrix. Much information is encoded in the Lyapunov exponents, being eigenvalues of the Lyapunov matrix

where the invertibility of the product is assumed throughout this work.

Let us consider an arbitrary product of operators, without prior assumption about their dependence. While our explicit calculations are for Ginibre matrices, our discussion should apply to more general situations, as well. As an example, in Fig. 2, we compare the microscopic level density in the bulk of the spectrum of the product of Ginibre matrices, complex Bernoulli matrices (with entries $0,\pm 1,\pm i,\pm 1\pm i$), and sums of Ginibre matrices $A_{j}$, $X_{j}=A_{j}+A_{j-1}$, with $A_{0}=0$. The latter choice can be understood as a model for progressive scattering with a short memory, making it non-Markovian. For $M=2$ it was considered in . Indeed, we could also choose $X_{j}=1+\gamma\delta t+A_{j}\sqrt{\delta t}$ with $A_{j}$ independent Ginibre matrices, the scalar $\gamma$ being proportional to the variance of $A_{j}$, and the time increment $\delta t\propto 1/M$. This model leads to the DMPK equation and in the limit $M\to\infty$ it models the Brownian stochastic process

In the present work, we will not go into details of this process and refer to .

First, we want to recall what is known. In several works, e.g. , it has been shown for various kinds of $X_{j}$ being i.i.d. that the distributions of the Lyapunov exponents $\lambda_{j}$ become asymptotically Gaussian when choosing $M\gg N$. The level density is then given by

with deterministic means $\bar{\lambda}_{j}$ and standard deviations $\sigma_{j}$. For $M$ Ginibre matrices they are given by

where $\psi(z)=(\log\Gamma(z))^{\prime}$ is the digamma function . When $N/M\to 0$ the level density becomes a sum of Dirac delta functions, corresponding to picket fence statistics after proper unfolding.

Clearly, the approximation (4) holds as long as the overlap between individual Lyapunov exponents is small, meaning that the $j$th width-to-spacing-ratio (WSR) defined by

is small. In the case of $M$ Ginibre matrices this is always satisfied when $j$ is fixed while $N,M\to\infty$ (hard edge scaling, smallest eigenvalue), regardless of any relation between $N$ and $M$.

In contrast, for the largest, $N$th Lyapunov exponent the ${\rm WSR}_{N}\approx\sqrt{N/M}$ is only small when $M\gg N$. For $M\ll N$ we cannot expect the approximation (4) to hold; the overlap between individual eigenvalue distributions becomes large, their interactions become important, and the standard GUE results should follow. The critical regime is when $M\propto N$ or more precisely $M\propto j$. In the ensuing sections, we will concentrate on this regime.

Summarizing the discussion above, we obtain the following picture illustrated in Fig. 1 for the product of $M$ Ginibre matrices. When studying the double scaling limit $M,N\rightarrow\infty$ with the relation $a=\lim_{N\rightarrow\infty}N/M$, we find three distinct asymptotic regimes: (i) $M=M(N)$ increasing super-linearly ($a=0$), (ii) linearly ($0<a<\infty$), or (iii) sub-linearly ($a=\infty$) with $N$. This identification of scales is our first main result. It shows us that there is always a critical spectral regime for $M\ll N$ in contrast to standard random matrix models and more natural for physical operators. Below we will only consider the bulk and soft edge (largest eigenvalue) of the Lyapunov exponents as they display new features. The local WSR is parametrised by $a$ as ${\rm WSR}_{j=Np}=\sqrt{ap}$, see (6). It increases when moving from left to right in the spectrum, cf. Fig. 1. A similar picture is expected for a general product of operators, with the WSR as a parameter, cf. Fig. 2.

Unfolding for Ginibre

To compare microscopic properties of different spectra one needs to unfold them via the averaged cumulative density

$\bar{R}_{1}(\lambda)$ is the non-fluctuating part of the level density, approaching the macroscopic level density at $N\to\infty$. When no analytical formulas for $\bar{R}_{1}(\lambda)$ and $\bar{N}(\lambda)$ are available, e.g. for comparison with experimental data, one has to perform a fit. Fortunately, $\bar{R}_{1}(\lambda)$ is known analytically for products of Ginibre matrices . In the bulk of the spectrum, taking $M\to\infty$ regardless of whether and how $N$ approaches infinity, this density is $\bar{R}_{1}(\lambda)=2e^{2\lambda}\Theta(N-e^{2\lambda})$, with $\Theta$ the Heaviside step function. This mapping will be explained in detail below. At the soft edge where $\bar{R}_{1}(\lambda)$ vanishes as a square root, we will specify the fit for the function $\bar{N}(\lambda)$.

The eigenvalue statistics of the random matrix $\exp[2ML]$ follows a determinantal point process (1) at fixed $M$ and $N$ with the kernel

Here, $G_{j}$ is essentially a Meijer G-function, and Eq. (8) is equivalent to the kernel derived in . The normalized level density of $L$ is

The prefactor stems from the Jacobian of the unfolding. Similarly, the $k$-point correlation functions (1) of $L$ are given by the kernel $2Me^{2M\lambda_{j}}K(e^{2M\lambda_{j}},e^{2M\lambda_{l}})$.

Due to the explicit form of the kernel (8) we can take various limits. For fixed $N$ and sufficiently large $M$ one reproduces the Gaussian result (4): For $y=e^{2M\lambda}$, the integral (9) picks up its main contribution from the saddle point at $t=0$. When $M\rightarrow\infty$ with finite $j$, this yields a Dirac delta at $\bar{\lambda}_{j}=\psi(j)/2$. Thence, the spectrum at the hard edge is always discrete and deterministic on a microscopic scale and, after proper unfolding, gives picket fence statistics.

The unfolding in the bulk follows from the asymptotic behaviour of the eigenvalues of the matrix $e^{2L}/N$ which are approximately $e^{\psi(j)}/N$. For

being of order $N$, the limiting form is $e^{\psi(j)}/N\approx p$ for $N\rightarrow\infty$, following from the asymptotic formula $\psi(Np)\approx\log(N)+\log(p)+\mathcal{O}(1/N)$. Thus, the eigenvalues $p=e^{2\lambda}/N$ of the matrix $\exp[2L]/N$ are uniformly distributed on $(0,1)$. This only holds when being away from the edges at $p=0,1$. The map $\lambda\mapsto p$ is just the cumulative distribution (7) (quantile) and therefore the proper unfolding.

Bulk Statistics

Zooming into the microscopic scale in the bulk at $p\in(0,1)$, we consider two neighbouring points

in the kernel (8), with $\xi$ and $\zeta$ of order unity. Multiplied by the Jacobian $M(Np+\xi)^{M-1}$, the kernel (8) in the bulk becomes

for $a=\lim_{N,M\to\infty}N/M\in(0,\infty)$. Here, $g(\xi)$ is an appropriate function to guarantee the existence of the limit, in this case $g(\xi)=(pN+\xi)^{Nj_{0}/a}$, with $j_{0}=[Np]$. Because it drops out from the determiant that yields the $k$-point correlation functions (1), we will not specify it any more below.

In , we give details of the derivation of (13). To sketch the idea, we apply the saddle point approximation and expand the summand and integrand in (8) and (9) about the summation index $j=Np+\delta j$ and the integration variable $t=i\delta t$, with $\delta j$ and $\delta t$ being of order one. This expansion is cumbersome but it is inevitable to go up to subleading orders since the leading orders cancel.

We would like to point out that the limit (13) can be also found in the spectrum when $M$ behaves sub-linearly with $N$ ($a=\infty$). There is always a small transitional regime close to the hard edge limit where the index $j$ of the Lyapunov exponent is of the order $M$, see Fig. 1.

Let us reformulate (13) by the Poisson summation formula to

The prefactor $\exp[(\xi^{2}-\zeta^{2})/(2ap)]$ can be skipped in the last formula since it drops out in all correlation functions (1), too. Though the second line of (13) has the advantage that each summand can be associated to a single Lyapunov exponent, the representation (14) reveals the intimate relation with the kernel from [45, Theorem 2.5], where Dyson’s BM with equidistant initial conditions has been studied, which is also a determinantal point process. This agreement is quite surprising, and analytically it shows the universality of our limit (13). In this particular BM has been identified with the solution of the DMPK equation. Thus, it can be expected that our result (13) corresponds to the local spectral statistics of the Anderson transition .

Let us discuss now the microscopic level density on the scale of the local mean level spacing, given by the relation

The behaviour of $\rho_{\rm bulk}(\xi;a)$ is shown in Fig. 3 for various values of $a$. In order to further corroborate the universality found above, a comparison to Monte Carlo simulations of three different products of random matrices is made in Fig. 2.

The interpolating microscopic density in the bulk (15) as well as the corresponding kernel (13) enjoy a discrete periodicity $\rho_{\rm bulk}(\xi;a)=\rho_{\rm bulk}(\xi+1;a)$ on a local scale. A full continuous translation invariance is restored only in the limit $a\rightarrow\infty$ where we analytically recover the sine-kernel

This follows from approximating the sum (13) by an integral that can be performed . Yet, there is always a small region between the hard edge and the bulk, even for $a\to 0$, where the transition kernel (15) can be still found because ${\rm WSR}_{j}$, see (6), or equivalently $ap=j/M={\rm WSR}_{j}^{2}$ enters the kernel (15).

For $a\rightarrow 0$ the density (15) reduces to a sum of Dirac delta functions,

since the erfi-function narrows to a Gaussian. The same applies to any $k$-point correlation function, yielding picket fence statistics .

Soft Edge Statistics

At the soft edge ($p=1$) the mean level spacing has a different scaling behaviour. Here, we adopt the finite $M$ scaling of and find the position of the upper spectral edge of $\exp[2ML]$ at $N^{M}(M+1)^{M+1}/M^{M}$, when $N\rightarrow\infty$. The probability of finding an eigenvalue above the edge drops exponentially. Therefore, we zoom into the spectrum as

The prefactor $a^{-2/3}$ in front of $\xi$ is reminiscent to the Airy-kernel scaling . We insert this scaling into (8) and fix the limit $a=\lim_{N,M\to\infty}N/M$. After relabelling the summation index $j\to N-j$ and expanding the summand as well as integrand in $1/N\sim 1/M$ in (8) and (9), the double scaling limit $M,N\to\infty$ (18) leads to the interpolating kernel,

our third main result. Again we postpone the detailed derivation to . The sum could be evaluated explicitly after the expansion.

An equivalent result was derived independently in . The interpolating microscopic level density at the soft-edge reads $\rho_{\rm soft}(\xi;a)=K_{\rm soft}\left(\xi,\xi;a\right)$. We also find that the soft edge kernel (19) appears in Dyson’s BM , and via the identification also for the DMPK equation. Therefore, also the interpolating kernel at the soft edge is universal and should hold for more general products of operators, as long as a soft edge is present.

For $a\to\infty$, we recover the GUE Airy-kernel ,

depending on the Airy function. For the interpolating microscopic density at the soft edge this implies

As for the Airy-kernel (20), the spectrum of (19) is not unfolded. For example, for large negative argument $\xi$, the limiting Airy-density (21) increases as $\sqrt{|\xi|}$ to the left, originating from the macroscopic semi-circular law of the GUE. In Fig. 4 we compare the unfolded density for the interpolating kernel (19) and the Airy-density. Our unfolding, fitting(7) by $\bar{N}(x)=ax+bx^{3/2}+cx^{2}$, is only valid inside the support of the macroscopic level density which vanishes as a square root, cf. . For the tails, which include the Tracy–Widom distribution for the GUE, a different means of comparison has to be sought.

In the opposite limit, for $a\to 0$, we rescale $\xi=\hat{\xi}/a^{1/3}$,

This is the picket fence density with a lower bound and agrees with the spectrum of a quantum harmonic oscillator. The rescaling by $a^{-1/3}$ is essential to get a normalized mean level spacing.

Comparing the observations above with the largest eigenvalue distribution, it is well known that the one corresponding to the Airy-density (21) is given by the Tracy–Widom distribution . In contrast, in the picket fence limit it is approximately Gaussian, see (4). This was also found in via a Fokker-Planck equation approach, and in . The interpolating kernel (19) thus describes an interpolation between the two, demanding further investigations.

Conclusions

The local statistical properties of Lyapunov exponents of transfer matrices, modelled by products of random matrices, have been discussed in the limit $N,M\to\infty$, where $M$ typically represents the number of time steps or layers, and $N$ the degrees of freedom and thus the complexity of the underlying system. The critical scaling $a\propto N/M$ separates two different phases that can be associated with integrable or chaotic behaviour. For $M\gg N$ the entire Lyapunov spectrum is deterministic. In the opposite limit $M\ll N$, the bulk and the largest Lyapunov exponents at the soft edge follow GUE statistics, associated with quantum chaotic behaviour. We filled the gap in analytically describing the local spectral statistics in the bulk and at the edge of the spectrum in the critical regime $M\propto N$, in deriving two limiting kernels that interpolate between the deterministic and GUE regimes. Interestingly, the transition between these two phases still exists for $M\ll N$ albeit not the whole spectrum lies in one or the other phase. The critical regime shrinks to a narrow scale about the hard edge.

Our findings agree with Dyson’s BM with picket fence initial conditions , showing the universality of our results. Here, we want to underline that this agreement is far from obvious since products of matrices cannot be easily expressed into sums of other independent objects like the sums of their logarithms. Indeed, this essential difference to scalars can be noticed by the validity of the map to Dyson’s BM, which only holds locally and not for the whole spectrum. For instance at the hard edge, the deterministic behaviour for the smallest Lyapunov exponents persists, regardless of what relation $N$ and $M$ have.

The relation between transfer matrices and the Dyson BM studied in has been already pointed out in for a specific choice of products of random matrices modelling the stochastic process (3). Therein, the product of matrices perturbing the unit matrix has been analysed which yields the DMPK equation in the limit of large $M$. Therefore, we believe that our interpolating regime is related to the Anderson transition, also resulting from the DMPK equation .

Since our model is analytically solvable for any finite $N$ and $M$, one can also address the nearest neighbour spacing distribution and the distribution of the largest Lyapunov exponent. The latter has important consequences for the stability of the system and its Kaplan-Yorke dimension . Another direction of investigation would be the analysis of the microscopic statistics of the complex eigenvalues of the product matrix. For the spectral radius first investigations were done , finding a transition, too.

References