Spectral Form Factor in Non-Gaussian Random Matrix Theories

TL;DR

This paper analyzes the spectral form factor in non-Gaussian random matrix models with rich phase structures, revealing multi-critical behavior and universality at large times, with implications for chaotic quantum systems.

Contribution

It provides explicit calculations of the spectral form factor in non-Gaussian matrix models, highlighting multi-criticality and universality in their spectral statistics.

Findings

Spectral form factor exhibits critical behavior at phase transition points.

Late time spectral correlations show universality across polynomial potentials.

Estimated dip-time varies with potential, indicating phase-dependent dynamics.

Abstract

We consider Random Matrix Theories with non-Gaussian potentials that have a rich phase structure in the large $N$ limit. We calculate the Spectral Form Factor (SFF) in such models and present them as interesting examples of dynamical models that display multi-criticality at short time-scales and universality at large time scales. The models with quartic and sextic potentials are explicitly worked out. The disconnected part of the Spectral Form Factor (SFF) shows a change in its decay behavior exactly at the critical points of each model. The dip-time of the SFF is estimated in each of these models. The late time behavior of all polynomial potential matrix models is shown to display a certain universality. This is related to the universality in the short distance correlations of the mean-level densities. We speculate on the implications of such universality for chaotic quantum systems including the SYK model.

Equations

$$|Z(\beta +it){|}^2=\sum _{m,n}{e}^{-\beta (E_m+E_n)}{e}^{-it(E_m-E_n)}$$

$$\mathcal{Z}=\left(\prod _{i=1}^N\int d{\lambda }_i\right)\exp \left[-\frac{N^2}{2}\left(\int d\lambda \rho (\lambda ){\lambda }^2-2\int d\lambda d\mu \rho (\lambda )\rho (\mu )\log (|\lambda -\mu |)\right)\right]$$

$$\frac{1}{2}{\lambda }^2=2{\int }_{\mathbb{R}}d\mu \rho (\mu )\log |\lambda -\mu |$$

$$\omega (\lambda \pm i0)=\frac{\lambda }{2}\pm \frac{1}{2}\sqrt{{\lambda }^2-4}$$

$$\rho (\lambda )=\frac{1}{2\pi }\sqrt{4-{\lambda }^2}$$

$$V(M)=\frac{1}{2}{M}^2+\frac{g}{N}{M}^4$$

$$\omega (\lambda )=\frac{\lambda }{2}+2g{\lambda }^3-\left(\frac{1}{2}+4g{a}^2+2g{\lambda }^2\right)\sqrt{{\lambda }^2-4a^2}$$

$$12g{a}^4+a^2-1=0$$

$$\rho (\lambda )=\frac{1}{\pi }\left(\frac{1}{2}+4g{a}^2+2g{\lambda }^2\right)\sqrt{4a^2-{\lambda }^2}$$

$$a^2=\frac{1}{24g}(\sqrt{1+48g}-1)$$

$$V(M)=\frac{1}{2}{M}^2+\frac{g}{N}{M}^4+\frac{h}{N^2}{M}^6$$

$$\omega(\lambda)=\frac{1}{2}\lambda+2g\lambda^{3}+3h\lambda^{5}-\bigg{(}3h\lambda^{4}+(2g+6ha^{2})\lambda^{2}+(\frac{1}{2}+4ga^{2}+18ha^{4})\bigg{)}\sqrt{\lambda^{2}-4a^{2}}$$

$$60h{a}^6+12g{a}^4+a^2-1=0$$

$$\rho(\lambda)=\frac{1}{\pi}\bigg{(}3h\lambda^{4}+(2g+6ha^{2})\lambda^{2}\hskip 85.35826pt\\ \hskip 42.67912pt+\left(\frac{1}{2}+4ga^{2}+18ha^{4}\right)\bigg{)}\sqrt{4a^{2}-\lambda^{2}}$$

$$g=\frac{1}{12a^4}(3-2a^2)h=\frac{1}{60a^6}(a^2-2)$$

$$g=-1/36h=1/1620a^2=3$$

$$|Z(\beta +it){|}^2=\sum _{m,n}{e}^{-\beta (E_m+E_n)}{e}^{-it(E_m-E_n)}$$

$$G(\beta,t)=\frac{\langle|Z(\beta+it)|^{2}\rangle_{\text{GUE}}}{\langle Z(\beta)\rangle_{\text{GUE}}^{2}}\hskip 91.04872pt\\ =\frac{\int\,d\lambda\,d\mu\,e^{-\beta(\lambda+\mu)}e^{-it(\lambda-\mu)}\big{\langle}\mathcal{D}(\lambda)\mathcal{D}(\mu)\big{\rangle}_{\text{GUE}}}{\int\,d\lambda\,d\mu\,e^{-\beta(\lambda+\mu)}\big{\langle}\mathcal{D}(\lambda)\rangle\langle\mathcal{D}(\mu)\big{\rangle}_{\text{GUE}}}$$

$$G_{dc}=\frac{\langle Z(\beta+it)\rangle\langle Z(\beta-it)\rangle}{\langle Z(\beta)\rangle^{2}}\bigg{|}_{\beta=0}\hskip 56.9055pt\\ \hskip 56.9055pt=\frac{\int\,d\lambda\,d\mu\,e^{-it(\lambda-\mu)}\big{\langle}\mathcal{D}(\lambda)\big{\rangle}\big{\langle}\mathcal{D}(\mu)\big{\rangle}}{\int\,d\lambda\,d\mu\big{\langle}\mathcal{D}(\lambda)\big{\rangle}\big{\langle}\mathcal{D}(\mu)\big{\rangle}}$$

$$G_c=(G-G_{dc}){|}_{\beta =0}$$

$$\mathcal{D}(\lambda )=\bar{\mathcal{D}}(\lambda )+\delta \mathcal{D}(\lambda )$$

$$G_{c}=\frac{\int\,d\lambda\,d\mu\,e^{-it(\lambda-\mu)}\langle\delta\mathcal{D}(\lambda)\delta\mathcal{D}(\mu)\rangle}{{\int\,d\lambda\,d\mu\big{\langle}\mathcal{D}(\lambda)\big{\rangle}\big{\langle}\mathcal{D}(\mu)\big{\rangle}}}$$

$${\int }_{-2\alpha }^{2\alpha }\mathcal{D}(\lambda )d\lambda =N$$

$${\int }_{-2a}^{2a}\rho (\lambda )d\lambda =1$$

$$\langle Z(\beta+it)\rangle_{\text{nGUE}}|_{\beta=0}=\int\,d\lambda\,e^{-i\tau\lambda}\big{\langle}\rho(\lambda)\big{\rangle}_{\text{nGUE}}\\ =\int_{-2a}^{2a}\,d\lambda\,e^{-i\tau\lambda}\frac{1}{\pi}\bigg{(}\frac{1}{2}+2g\lambda^{2}+4ga^{2}\bigg{)}\sqrt{4a^{2}-\lambda^{2}}\\ \\ =\frac{1}{\tau^{2}}\bigg{(}a\,\tau\,(1+24a^{2}g)\,J_{1}(2a\tau)-24a^{2}g\,J_{2}(2a\tau)\bigg{)}$$

$$a^2=\frac{1}{24g}(\sqrt{1+48g}-1){|}_{g=g_c}=-\frac{1}{24g_c}$$

$$⟨Z(\beta +it){⟩}_{\text{nGUE}}{|}_{\beta =0}\stackrel{t\to \infty }{\to }-\sqrt{\frac{a}{\pi }}\frac{(1+24a^{2}g)}{{\tau }^{3/2}}\cos \left(\frac{\pi }{4}+2a\tau \right)-\frac{1}{\sqrt{a\pi }}\frac{(1+152a^{2}g)}{{\tau }^{5/2}}\sin \left(\frac{\pi }{4}+2a\tau \right)+\mathcal{O}\left(\frac{1}{{\tau }^{7/2}}\right)$$

$$\rho (\lambda )=\frac{1}{12\pi a^2}(4a^2-{\lambda }^2{)}^{3/2}$$

$$\displaystyle\langle Z(\beta+it)\rangle_{\text{nGUE}}|_{\beta=0}=\int\,d\lambda\,e^{-i\tau\lambda}\big{\langle}\rho(\lambda)\big{\rangle}_{\text{nGUE}}$$

$$\displaystyle=\int_{-2a}^{2a}\,d\lambda\,e^{-i\tau\lambda}\frac{1}{\pi}\bigg{(}3h\lambda^{4}+(2g+6ha^{2})\lambda^{2}+\bigg{(}\frac{1}{2}+4ga^{2}+18ha^{4}\bigg{)}\bigg{)}\sqrt{4a^{2}-\lambda^{2}}$$

$$\displaystyle=\frac{1}{\tau^{4}}\bigg{(}a(-360a^{2}h\tau+(1+24a^{2}g+180a^{4}h)\tau^{3})\,J_{1}(2a\tau)-24a^{2}(-30h+(g+15a^{2}h)\tau^{2})J_{2}(2a\tau)\bigg{)}$$

Full text

††preprint: TIFR/TH/17-24

I Introduction

Quantum chaos in many body systems has attracted the attention of physicists across various disciplines for over many decades now. Given that classically chaotic systems could be characterised using only a few stringent conditions, one hoped that the same would be true for chaotic quantum systems. However, it turns out that chaos in many body quantum systems is a much more complicated phenomenon. A key role to characterise chaos in quantum systems is played by Random Matrix Theories (RMTs). The idea (first proposed in Dyson (1962)) is to look at the energy spectrum of a given quantum system. Typically, for a few body system, if the levels follow the Wigner’s semi-circle law, then such systems are very well approximated by RMTs. For many body quantum systems, a sharp diagnostic of chaos is obtained by looking at the distribution of the nearest neighbour spacings (NNSD) of the energy levels of the system. The system is termed chaotic if the distribution is of the Wigner-Dyson type Bohigas et al. (1984)111See Garcia-Garcia and Verbaarschot (2016) for a discussion of the SYK model.. Such a behaviour is considered to be the hallmark of chaos in quantum many-body systems. Interestingly, the above expression has been shown to hold analytically only for few-body Hamiltonians ($N$ finite). For many-body quantum systems, the result has been established only numerically. Nonetheless, the remarkable validity of this result over a wide range of systems leads one to believe that chaos at the quantum level is an extremely rich phenomenon, that may exist irrespective of whether the system has a classical counterpart.

That being said, classically chaotic systems can still teach us a thing or two about quantum chaos. An important diagnosis for chaos in classical systems is the sensitive dependence on initial conditions, characterised by a positive Lyapunov exponent. An analog of this definition for a few-body quantum system was outlined in Larkin and Ovchinnikov (1969) using out-of-time-ordered correlators (OTOCs). Recently, the OTOCs have been widely adopted for the analysis of chaos in quantum many body systems, including large $N$ CFTs and even Black Holes Kitaev (2014); Maldacena et al. (2016); Shenker and Stanford (2014); Roberts and Stanford (2015).

The novelty of OTOCs is in that it helps us explore an interesting timescale, called scrambling time, in systems with a large number of degrees of freedom. First discussed in the context of black holes in Hayden and Preskill (2007); Sekino and Susskind (2008), this timescale characterises how long it effectively takes for an initial localized perturbation to spread over all degrees of freedom in a system. Usually, for chaotic many body systems, such a timescale is parametrically larger than the thermalisation timescale in the system. In order to make any connection between the two afore-mentioned definitions of quantum chaos, it would be worthwhile to understand whether and how scrambling of degrees of freedom ultimately leads to the onset of RMT behaviour in the system. In fact, some concrete steps along this direction have already been undertaken in the context of the SYK model. The SYK model was first shown to be chaotic using OTOCs in Kitaev (2014); Maldacena and Stanford (2016). More recently, the model was shown to have a late-time behaviour similar to Gaussian RMTs Cotler et al. (2016). However, the analysis in Cotler et al. (2016) was mostly numerical, and a more analytic understanding of this cross-over remains an open problem.

An intriguing aspect of RMTs is an inherent timescale, called the ramp time Gharibyan et al. (2018). It is the timescale below which no universal correlations between the energy levels exists. It usually indicates the onset of Random Matrix behaviour in the system. Above the ramp time, the universal correlations between the nearby energy levels starts to manifest themselves. A good diagnostic for exploring such timescales is an observable called the Spectral Form Factor (SFF). The SFF was introduced long ago Brezin and Hikami (1996a) and recently re-introduced Garcia-Garcia and Verbaarschot (2016); Dyer and Gur-Ari (2016); Cotler et al. (2016); Krishnan et al. (2017a); Balasubramanian et al. (2016); del Campo et al. (2017); Krishnan et al. (2017b); Mehta (2004); Papadodimas and Raju (2015) as a tool for probing the spectra of quantum systems. It is defined as,

[TABLE]

where $Z(\beta)$ is the partition function of a quantum system and $\beta$ is the inverse temperature. For $\beta=0$, it is easy to see that the above expression picks out contributions only from the differences between the nearest neighbour energy eigenvalues at very late times. This makes it a good probe for understanding not only the discreteness of the energy spectrum, but also gives a way to characterise a chaotic system based on its NNSD. In fact, the SFF (when averaged over Gaussian Random Matrices) has a very particular behavior at large $N$ characterised by an intial decay, followed by a linear rise and saturation. This indicates how the SFF perceives the nearest neighbour energy spectrum. In fact in Cotler et al. (2016), the SFF was used to show the similarities between the late time behaviour of the SYK model and Gaussian RMTs.

The goal of this paper is to understand whether and how the structure of the SFF is modified under a controlled deformation of the RMT from a Gaussian to a non-Gaussian model. For this, we will consider the SFF corresponding to polynomial non-Gaussian matrix models, that have a richer phase structure as compared to the Gaussian model222Phase transitions, corresponding to such a phase structure, were first shown to exist in Unitary Matrix Models and are termed the Gross-Witten-Wadia (GWW) phase transitions Wadia (1980); Gross and Witten (1980).. Simultaneously, we will try to understand how such deformations affect the chaotic structure of the energy eigenvalues. This shall be reflected in any changes in the onset time for the ramp in the SFF. We specifically work out the Quartic and the Sextic models, where we average the SFF over each of these non-Gaussian ensembles and make the following observations,

The initial decay behavior is the same as for the Gaussian ensemble, when the non-Gaussian ensemble is away from its critical points. 2. 2.

At the critical points, the decay behavior of the SFF changes over to a faster decay, displaying a different power law behavior. This gives an explicit example where the multi-critical behavior at equilibrium is extended to non-equilibrium dynamics. 3. 3.

The dip-time estimate changes with the change in the fall-off behavior at criticality, providing an upper bound on the ramp-time. This in turn provides us with an idea of when RMT behaviour sets into the system. 4. 4.

The late time behavior of the SFF continues to be characterised by a linear growth followed by saturation, indicating no change in the NNSD.

We finish by speculating on how to obtain a possible analytic understanding of the crossover between the two different descriptions of chaotic behaviour in the same system.

This paper is organized as follows. In II, we start with describing the Gaussian Matrix Model and calculate its mean level density in the large $N$ limit. That is followed by a description of two specific non-Gaussian Matrix Models, namely the quartic and the sextic models. We then proceed to calculate the mean-level densities for these models. In III, we introduce the Spectral Form Factor and the method for calculating it in the large $N$ limit, by separating it out into a connected and a disconnected part. In this section, we state our main results regarding the change in the decay behavior of the disconnected part of the SFF exactly at the critical points of these models. In IV we estimate the dip-time (indicating the onset of ramp) for these models. In V, we talk about the implications of non-Gaussian Matrix Models on chaotic systems.

II Matrix Models: Gaussian and Non-Gaussian

The approach to the complete non-perturbative solution of RMTs with the most general polynomial potential (at large $N$) is well establishedMehta (2004)333See Appendix A and B for a comprehensive review of the solution.. Equipped with this solution, we shall now consider specific forms of the potential and study the corresponding matrix models case by case to understand their properties. The idea would be to understand their implications in describing chaotic behavior, in particular, the spectral form factor, which we shall define in the next section.

II.1 Gaussian Matrix Model

We begin by considering a Gaussian potential: $V(M)=\frac{1}{2}M^{2}$. The corresponding action in the continuum limit is,

[TABLE]

Extremising this action with respect to the level density $\rho(\lambda)$ gives us,

[TABLE]

Taking a derivative of the overall equation w.r.t. $\lambda$ gives (44) with $V^{\prime}(\lambda)=\lambda$. For a Gaussian potential in a symmetric interval \big{(}-1,1\big{)}, the resolvent takes the following form,

[TABLE]

From the properties of the resolvent (stated in Appendix B), it follows that the level density takes the unique form,

[TABLE]

This is the Wigner’s semi-circle law for Gaussian matrices (shown in Fig. 1) in the large $N$ limit.

II.2 Non-Gaussian Matrix Models

Quartic Potential

As the simplest example of a non-Gaussian matrix model, we consider a quartic potential of the form,

[TABLE]

where, in the second term, $g$ is the coupling constant and the power of $N$ is so adjusted as to make the quartic interaction strength comparable to the quadratic term in the action at large $N$. This model has been shown (Brezin et al. (1978); Bhanot et al. (1990); Mandal (1990)) to have a critical point at $g=g_{c}=-1/48$, below which it is supposed to have no solutions. At $g_{c}$ the free-energy of the model has a non-analyticity. It is our goal to study the behavior of the spectral form factor computed using one-cut solution444An $n$-cut solution is a function which has support on $n$ disjoint intervals. Throughout this paper, we will work with one cut solutions only. and see how it changes at the critical point. It can be argued (Appendix B) that there is a unique resolvent corresponding to the saddle point equation in this phase. This unique form of the resolvent is,

[TABLE]

where $a$ is defined through the constraint,

[TABLE]

The level density is found from (52) to be,

[TABLE]

The end-points of the interval $(-2a,2a)$ (on which $\rho(\lambda)$ is real) depend on the value of the coupling constant $g$ through (8) as,

[TABLE]

It is easy to see from above that $a^{2}$ will have imaginary values for $g\leq-1/48$ and hence, there are no real solutions of (8) that exist below $g_{c}=-1/48$.

One shortcoming of this analysis is the absence of two cut solutions, making it impossible to access the phase below $g=g_{c}$. To remedy this, we need to go to a matrix model with a sextic potential Bhanot et al. (1990). Note that at $g=g_{c}$, the behaviour of the level density, near the edge, changes to (23) 555We will consider the change in the edge behaviour of the level density as a sign of criticality..

Sextic Potential

The form of the potential in the sextic model is,

[TABLE]

The unique form of the resolvent for the one-cut solution with the above potential is given by,

[TABLE]

with $h$ being defined through the constraint,

[TABLE]

The level density, in the one-cut phase, corresponding to (12) is,

[TABLE]

From (13) we see that the end points of the interval now depend on both the coupling constants $g$ and $h$. To find the critical points in this model, we shall look for the points where the behavior of the level density changes near the edges of the spectrum. This calculation has been carried out in Appendix C. In this model, there is a line of critical points in the $(g,h)$ plane, which starts from $(g,h)=(-1/48,0)$. The values of $g$ and $h$ on this critical line are Bhanot et al. (1990),

[TABLE]

For $a^{2}=2$, the above values reduce to the critical point of the quartic model. The critical line encounters a tri-critical point for the following values of the couplings,

[TABLE]

Eliminating $a$ from (14), we can construct a function $g_{c}\equiv g(h)$ along the critical line. Now, keeping $h$ fixed, if we start decreasing the value of $g$ starting from the origin, we find a number of solutions for the level density organised in the following way. For $g>g_{c}(h)$, and a fixed value of $h$, there exist all one, two and three-cut solutions. However, for $g<g_{c}(h)$, the one-cut solution vanishes, leaving behind the two and three-cut solutions which co-exist. From a naive analysis of the free-energy, it can be shown that in the region $g>g_{c}(h)$, the one-cut solution is the most stable, thereby dominating the thermodynamic behaviour. The nomenclature of the one-cut phase for $g>g_{c}(h)$ and the two-cut phase for $g<g_{c}(h)$ is, thereby derived from such a behvaiour.

III Spectral Form Factor

The spectral form factor(SFF) was originally introduced as a probe for analysing the spectrum of quantum systems. It is defined in terms of the analytically continued partition function as,

[TABLE]

From the R.H.S. of the above expression, it is easy to see that at very high temperatures $(\beta\rightarrow 0)$, the SFF is highly sensitive to the differences in the energy levels of the system. Interestingly, only the nearest neighbour energy spacings contribute to the very late time behavior. This enables the SFF to play a crucial role in understanding the time dynamics of chaotic quantum systems and also makes it a very useful tool for probing the discreteness in their energy spectrum. A property of chaotic quantum systems is that they satisfy the Wigner’s surmise, which is a statement about the distribution of the nearest neighbour spacings of the energy levels. This qualifies the SFF as a suitable observable for analysing chaos in such quantum systems.

As we learnt in the previous sections, a lot is known about the energy spectrum of RMTs in the large $N$ limit. When the SFF is averaged over an ensemble of random matrices, it has a very particular behavior, especially at late times Cotler et al. (2016); Brezin et al. (1978). In this light, we wish to understand the general dependence of the SFF on the specific nature of the RMTs and understand the time-scales for the onset of RMT behaviour in quantum chaotic systems. The idea is to consider the various non-Gaussian ensembles with polynomial potentials discussed in the previous section, and average the time-dependent SFF over all these ensembles. All such ensembles have a distribution of eigenvalues that are different from each other, except at small energy scales. We wish to understand in what sense the SFF is sensitive to the difference in the ensembles, which would tell us its effectiveness in analysing the onset of quantum chaos. The quantity we wish to study is defined as,

[TABLE]

where $\mathcal{D}(\lambda)$ is the eigenvalue density and the integration is over the (unscaled) support of $\rho(\lambda)$ in each case. For simplicity and clarity of the analysis, we wish to focus on the $\beta=0$ (or $T=\infty$) window. In the large $N$ limit, the two point function of the level densities in (III) can be divided into two parts; the connected part ($G_{c}$) and the disconnected part ($G_{dc}$), which we define below.

[TABLE]

[TABLE]

We define the quantity (17) as the connected two-point correlation function. In fact, it is easy to show that if we consider fluctuations around the extremised value of the level density,

[TABLE]

then the connected two point correlation function can be written as just $\langle\delta\mathcal{D}(\lambda)\delta\mathcal{D}(\mu)\rangle$. The connected part of SFF can then be re-written as,

[TABLE]

We wish to analyse the behavior of these two parts separately in RMTs with non-Gaussian potentials. We shall consider the quartic and sextic potentials and analyse the behavior of the SFF, when averaged over these ensembles, using expressions for the level density we found in the previous sections.

Convention

For calculating the SFF, we shall use the mean level density in the one cut phase of the quartic model, given by (9). Before proceeding, we shall put in place some comments regarding the conventions we will be using, in order to connect up to results in the rest of the literature. The mean level density can be normalised in two ways.

[TABLE]

where $\mathcal{D}(\lambda)$ counts the total number of eigenvalues between $(\lambda,\lambda+d\lambda)$ and $\alpha\sim\mathcal{O}(\sqrt{N})$. This is also the expression we mean to use for calculating the SFF (since $\langle Z(\beta+it)\rangle|_{\beta,t=0}=N$). However, this is not the mean level density we get from extremising the action since all factors of $N$ have been scaled out from our action and all eigenvalues are just $\mathcal{O}(1)$. The density we get from our action is normalized as,

[TABLE]

where the interval parameter $a\sim\mathcal{O}(1)$. The two are related by a $\sqrt{N}$ scaling of the eigenvalues. If we rewrite $G_{c}$ and $G_{dc}$ in terms of this scaled density, all expressions would remain the same (with appropriate scalings for the coupling constants wherever necessary) except for two changes,

The intergration limits would correspond to $a\sim\mathcal{O}(1)$. 2. 2.

The time shall be scaled accordingly as $\tau=\sqrt{N}t$.

Using the scaled mean level density we then proceed to calculate the SFF.

III.1 The Disconnected Part : $G_{dc}$

Quartic Model

To calculate the disconnected part of the SFF, we use the scaled mean level density of the one-cut potential (9) into the first term of (III).

[TABLE]

The above expression 666Note that the expression for SFF on the L.H.S. is in terms of $t$ but all expressions on the R.H.S. are in terms of $\tau$. has a very intriguing form. The reason being that the coefficient of the first term vanishes exactly at the thermodynamic critical point, $g=g_{c}(=-1/48)$, of the quartic model, since

[TABLE]

This suggests that it is only the Bessel function of the second kind that governs the fall off behavior of the SFF at the critical point. This has the non-trivial implication that the SFF has a faster power law decay at the critical point as compared to any other value of $g$. To understand this better, we look at the large time asymptotics of the above expression and find,

[TABLE]

This clearly suggests that for all values other than the critical value, the power law assumes a $\tau^{-3/2}$ behavior. However, at $g=g_{c}$ the coefficient of the $\tau^{-3/2}$ term vanishes indicating a crossover to a $\tau^{-5/2}$ power law decay behavior.

One way to explain this behavior is to look at the form of the level density $(\rho(\lambda))$ at the critical point. Away from the critical point, the level density has the form (9) such that it goes to zero near the $\lambda=2a$ edge of the interval as $(\lambda-2a)^{1/2}$. This behavior results in the $\tau^{-3}$ power law beahviour of the SFF. At the critical point $g_{c}$, the SFF displays the following behavior,

[TABLE]

This implies that at the $\lambda=2a$ edge, the level density now goes to zero as $(\lambda-2a)^{3/2}$. The change in the power law decay behavior can then be completely attributed to this particular change in behavior of the level density near the edge of the spectrum. One can check this behaviour to be in accord with the Paley-Weiner theorem.

Sextic Model

For the sextic model, the mean level density in the one-cut phase is given by (II). We shall work with the expression for the level density in the one cut phase and predict all critical points in the theory accessible from this phase. Calculating the SFF in this model, we get,

[TABLE]

This expression is subject to a constraint equation,

[TABLE]

It is easy to check that putting $h=0$ in the above expression gives us back (III). To check that our method correctly predicts the critical points (as it did in the quartic model), we shall try to predict the critical points of the sextic model in a similar way. We shall look at the vanishing of the coefficients of the various powers of $\tau$ from the expression for the large time expansion of (24).

The large time behavior of the SFF is,

[TABLE]

For arbitrary values of $h$ and $g$, (26) will have a large time fall off behavior governed by the $\tau^{-3/2}$ term. However, to see the transition to a faster decay we would expect the coefficient of $\tau^{-3/2}$ to vanish while satisfying (25). Solving the set of simultaneous equations for $g$ and $h$, we find

[TABLE]

This is exactly the same as (14). This, of course, is a one-parameter set of solutions. One has a solution for every value of $a\geq\sqrt{2}$. The lower bound is the value where, $h=0$ and $g=-1/48$ which are the values for the critical point of the quartic model. The vanishing of this coefficient suggests that there is a transition of the power of the decay from $\tau^{-3/2}$ to $\tau^{-5/2}$. The solution in terms of an unfixed parameter suggests that we have a critical line separating the one and two cut phases of the sextic model. There is no upper bound on the value of $a$ yet. That can be obtained by analysing the point where the coefficient of the next highest power of $\tau$, namely $\tau^{-5/2}$, vanishes. Solving for the value of $a$ by plugging in the values of $g$ and $h$ from (27) into the coefficient of the $\tau^{-5/2}$ term, we get $a=\sqrt{3}$ and correspondingly,

[TABLE]

This is, again, the exact same point as in (15). This indicates the place where the critical line, separating the one and two cut phases, encounters a tri-critical point. It turns out that exactly at this point, the long time fall off behavior of the SFF changes over from $\tau^{-5/2}$ to $\tau^{-7/2}$.

Again, this change in the fall-off behavior can be attributed to the change in the manner in which the mean level density approaches the edge of one of its supports. Away from the critical point, $\rho(\lambda)\sim(\lambda-2a)^{1/2}$ near the $\lambda=2a$ edge. Along the critical line and upto the tri-critical point,

[TABLE]

which implies that $\rho(\lambda)\sim(\lambda-2a)^{3/2}$ at the $\lambda=2a$ edge777The transition in the power-law behavior of the SFF, to a decay faster than $1/\tau^{3}$near criticality, as a result of the change in the edge behavior of the density of states was hinted at in Torres-Herrera et al. (2018). We would like to thank Antonio M. Garcia-Garcia for pointing this out to us. . Further, at the tri-critical point, the mean level density becomes

[TABLE]

This clearly implies that $\rho(\lambda)\sim(\lambda-2a)^{5/2}$ at the $\lambda=2a$ edge. This then is the change responsible for the change in the behavior of the SFF from $\tau^{-5/2}$ to $\tau^{-7/2}$.

III.2 Universality of the connected part: $G_{c}$

III.2.1 Universality of $G_{c}$ under the choice of potential

The connected part of the SFF is as defined in (18). It depends on the connected part of the two-point function of the mean level densities. In fact, the connected part depends only on the correlations of the fluctuations $\delta\rho(\lambda)$ around the mean-level densities. Let us use this fact to calculate some universal properties of $G_{c}$. Consider the continuum action (2) but with a general potential $V(\lambda)$ (replace $\lambda^{2}$ with $2V(\lambda)$ in (2)). Let us calculate correlations of the fluctuations $\delta\rho(\lambda)$ around the mean-level density perturbatively. Let $\rho(\lambda)=\rho_{0}+\delta\rho(\lambda)$. Here $\rho_{0}$ is the large $N$ mean-level density. Plugging this into (2) we get

[TABLE]

where $S_{0}$ will be canceled by denominator in the definition of $G_{c}$ (18), $S_{1}(\delta\rho)=0$ as $\rho_{0}$ is a large $N$ saddle-point solution, and

[TABLE]

Going over to Fourier space, we get

[TABLE]

With this expression, we identify (18) to be $|\tau|/(2\pi N^{2})$. This result then explains the ramp in SFF. In order to analyse very large time behavior of the SFF, we take a long time average of SFF and get

[TABLE]

where $g(\lambda)$ is the degeneracy in $\lambda$. For a non-degenrate spectrum, $g(\lambda)=1\ \forall\ \lambda$, which then gives

[TABLE]

The RHS, in the above equation, is independent of time. This result then explains the origin of a time-independent plateau in SFF.

Note:* The ramp and plateau behavior of the SFF is independent of the choice of potential. The connected part of the SFF displays a universal behavior under the choice of any polynomial potential, as long as the potential has the symmetries of the respective ensemble.*

III.2.2 Universality of short distance correlations

There is a more sophisticated way to understand this result. A very standard result of RMT Mehta (2004) is the exact form of this correlator near the centre of the spectrum of the eigenvalues,

[TABLE]

This result is dubbed the sine kernel and was originally derived using the method of Orthogonal polynomials for the Gaussian ensembles. However, we wish to emphasize that this result holds true for any polynomial potential measure (of single trace operators) from which the ensemble of matrices is chosen. This point was previously mentioned in Brezin and Zee (1993). In the context of the SYK model, it then suggests that all such polynomial potential models are equally good for describing its late time behavior. A physical reason is that the various polynomial potentials significantly change the eigenvalue spectrum only near the edges(Fig.2 and 3). The correlations of the level densities near the center of the spectrum ($|\lambda|\sim 0$), however, are observables that are oblivious to the change in the global structure of the eigenvalue density of the system, especially near the edges, and hence remain unchanged.

III.2.3 Fourier Transform

To get to the SFF from (35), we need to plug it into the expression for the connected part in (18)888It should be noted that the scaling of $\mathcal{D}(\lambda)$ has been cancelled out between the numerator and the denominator, such that the same expression is valid with $\mathcal{D}(\lambda)$ replaced by $\rho({\lambda})$.. It is easy to see that with a change of variables from $(\lambda,\mu)$ to $(E=\lambda+\mu,\omega=\lambda-\mu)$, we can perform the integration. There are two parts to this integral, which we shall deal with separately.

The $1/N^{2}$ part with the sine squared function that is responsible for the ramp. It is also the sub-dominant contribution. 2. 2.

The $1/N$ part with the Dirac-delta function (the dominant part) which is reponsible for the plateau.

The Ramp and Plateau

The integral involving the first term in (35) is responsible for the linear rise in time of the SFF, as mentioned earlier. We wish to show why this is the case. The expression for the SFF is given by,

[TABLE]

where $\omega$ denotes the difference in energy levels. Evaluating the above integral gives us,

[TABLE]

These expressions tell us about the very long time behavior of the SFF at a fixed energy. At $\tau<2N$, the SFF has a linear growth in time. At $\tau\geq 2N$, this linear behavior saturates abruptly to a plateau region.

From a physical point of view, it is enlightening to understand how this behavior arises. In (35), we can re-define the argument of the sine term in terms of,

[TABLE]

We would like to work in the limit where,

[TABLE]

The scaling parameter $x$ in (37) can be made small or large (by changing the way $\omega$ goes to zero) depending on the region of the spectrum we wish to focus on.

For large $x(\gg 1)$, $(\sin(x)/x)\sim 1/x$ in (35). This behaviour implies that correlations between energy levels that are far apart, are suppressed as $1/N^{2}$. These fluctuations, in turn, do not contribute to the SFF. This behavior is dubbed the spectral rigidity. 2. 2.

For small $x(\ll 1)$, $(\sin(x)/x)\rightarrow 1$ in (35). This scaling focuses on the correlations between fluctuations in the nearby energy levels. The integral in (36) gets its maximum contribution from near the $\omega=0$ region. The non-zero value of the $(\sin(x)/x)$ part is responsible for the time dynamics of the SFF. These differences in the nearest neighbour energy eigenvalues are therefore the ones that cause the ramp behavior of the SFF. Evidently, the ramp stops when we reach the $\omega$ corresponding to values lower than the smallest eigenvalue spacing in the system, i.e. $\omega<\mathcal{O}(1/N)$.

IV Estimation of dip time

One of the non-trivial effects of considering the SFF in a non-Gaussian RMT is the change in its initial fall-off behavior exactly at the critical points. This change also changes the dip-time for the various models we consider. As a corollary, this behaviour helps us uncover the timescale when the ramp starts, thereby giving us an upper bound on the ramp time Gharibyan et al. (2018). The dip-time is estimated by comparing the initial fall-off behavior with the late-time behavior of the curve where it starts to rise linearly. We shall start by estimating the dip-time in the Gaussian model just to jog our memory, and quickly proceed to estimate it at the various critical points in the quartic and the sextic models.

IV.1 Gaussian RMT

Equating the decay time from the disconnected part to the linear rise time from the connected part of the SFF in the Gaussian model, we find

[TABLE]

However, we should remember that $\tau=\sqrt{N}t$. Thus, in terms of the physical time parameter,

[TABLE]

which is the estimated dip-time for this model.

IV.2 Non-Gaussian RMT

Quartic Model

Away from the critical point: Equating the decay rate with the linear rise, away from the critical point of the quartic model, we see that

[TABLE]

which is the same as in the Gaussian model. Hence, the estimated dip time here would also be

[TABLE] 2. 2.

At the Critical Point: Equating the decay time with the linear rise at the critical point, we find

[TABLE]

In terms of the physical time parameter, the estimated dip-time at the critical point is

[TABLE]

Sextic Model

Away from criticality: Away from criticality, the estimated dip-time for this model is the same as that for the previous models, namely

[TABLE] 2. 2.

On the critical line: On the critical line for this model, we get

[TABLE]

This behavior is the same as that at the critical point of the Quartic model. The estimated dip-time in terms of the physical parameter is then,

[TABLE] 3. 3.

At the tri-critical point: At the tri-critical point,

[TABLE]

In terms of the physical time parameter, the estimated dip time is,

[TABLE]

V Onset of Random Matrix Behaviour and Quantum Chaos

The SFF has been developed as a tool for probing the nearest neighbour eigenvalue spacings in quantum systems. Due to the universality of the short distance correlations between eigenvalues, the very late time behavior of the SFF appears to be completely universal Mehta (2004); Brezin and Zee (1993). In Section III.B.1, we claim that this universality holds for all polynomial potentials $V(\lambda)$. In this section, we wish to point out the implications of this result on chaotic quantum systems via the SFF.

One very good example of a chaotic system is the SYK model. The SFF of this model, averaged over a Gaussian ensemble of the couplings, has been shown to behave like a Gaussian RMT at large times Cotler et al. (2016). On another instance, the authors in Krajewski et al. (2018) evaluated the large $N$ effective action for the SYK model averaged over a non-Gaussian ensemble of the couplings. The effect of such an averaging was found to shift the variance (proportional to $J$) by a term proportional to the non-Gaussian coupling. For the late time behaviour of the SFF, this would imply that the time-scale for the onset of ramp in the SFF will be affected (see Fig. 7 in Cotler et al. (2016)), without affecting the plateau in anyway. This is precisely our conclusion from the general analysis we carried out. To summarise our observations of universality in the late time behavior of the SFF averaged over arbitrary non-Gaussian ensembles with polynomial potentials,

Any RMT with a polynomial potential can describe a chaotic quantum system (including the SYK model) so long as we concern ourselves with the low energy eigenvalues and their correlations. This behavior was also suggested in Verbaarschot and Wettig (2000), where RMT with arbitrary polynomial potential was shown to describe the IR behavior of QCD 999We thank Nikhil Karthik for bringing this work to our attention. 2. 2.

The NNSD structure of the polynomial potential models seems to be similar to the Gaussian model, implying that they too are well described by the Wigner’s Surmise. In the absence of analytic calculations, we shall present some numerical evidence for the suggested behavior in Appendix D.

Another interesting aspect of the non-Gaussian averaging is the early time behaviour and actual onset of the ramp in chaotic systems. Away from criticality, the complete SFF does not, in any way, distinguish between the Gaussian and non-Gaussian ensembles. However, at the critical points, the early time decay of the SFF switches over to a faster fall-off behaviour. This change shows that the time scale at which ramp begins is much earlier than the dip time. In other words, it demonstrates a straightforward manner in which we can explore the actual ramp time in matrix models (see Gharibyan et al. (2018) for other ways to so).

VI Conclusions

We have studied the Spectral Form Factor (SFF) in non-Gaussian random matrix models. We found the initial decay behavior of the SFF to be the same in both the Gaussian and non-Gaussian models except at the critical points of the latter class, where it was found to display a faster decay. We related this change in behavior to the change in the behavior of the mean eigenvalue density at the edges of its support, exactly at the critical points of the model. Moreover, at criticality, this faster decay helps us get a better upper bound on the ramp time. This is important for the understanding of the onset of random matrix behaviour in a quantum chaotic system.

A crucial observation was made for the late time behavior of the SFF. Beyond the dip-time the SFF turned out to be universal for any polynomial potential. This behavior was due to spectral rigidity. With regard to implications on chaotic quantum systems, we emphasised that any random matrix models with a polynomial potential seems good enough to describe the late time chaotic behavior in such models.

Due to this universal behavior of the SFF, one might think that the nearest neighbour spacing distribution in all such models with polynomial potentials are well described by the Wigner-Dyson distribution. We provided some numerical simulations to support this claim. However, analytic techniques to calculate the NNSD of the non-Gaussian models remains an open problem. There are hints that this universal behavior persists for non-polynomial potentials as well. One interesting future direction would be to include double-trace operators and calculate the late-time behavior in such cases, to check if the aforesaid universality still persists. One very important goal of this project was to shed some light on how a chaotic quantum system displaying a typical scrambling behaviour at early times, gives way for the onset of a random matrix-like behaviour at late times. In this regard, we have been able to explore the various time-scales when RMT like behaviour might set into the system. However, in order to truly understand the cross-over, one would require an understanding of the phenomena of scrambling between the scrambling time and the ramp time, probably in a system with finite number of degrees of freedom. Scrambling on these time-scales would probably set the stage for the onset of an RMT description of the system beyond ramp time.

In Dijkgraaf and Vafa (2002a, b), a direct connection was made between matrix models (with non-Gaussian terms) and supersymmetric gauge theories, by identifying the superpotential of the gauge theory with the potential of the matrix model. One has to work with the gauged version of the matrix model to use this equivalence. However, it would be rather interesting to explore whether, in the spirit of Dijkgraaf et al. (2003), it is possible to directly explore the non-perturbative spectrum of the gauged SUSY theory from the perturbative dynamics of the matrix model, and thereby try to shed some light on their chaotic behavior.

Acknowledgements.

We thank Gautam Mandal for suggesting the problem and also for numerous important discussions during the course of this work. We also thank Saumen Datta, Sambuddha Sanyal, Geet Ghanshyam, Sounak Biswas and Sujoy Chakraborty for discussions on various aspects of the problem. Our computational work relied exclusively on the computational resources of the Department of Theoretical Physics at the Tata Institute of Fundamental Research (TIFR). R.S. would like to thank the Galileo Galilei Institute for Theoretical Physics for the hospitality and the INFN for partial support during the completion of this work. This work was also partly supported by Infosys Endowment for the study of the Quantum Structure of Space Time.

Appendix A Overview of RMT

Gaussian Matrix Ensembles are created by considering a large number of matrices, each of which are filled with random numbers drawn from a Gaussian distribution. Depending on whether the matrix elements are real, complex or quaternion numbers, we can have an orthogonal (GOE), unitary (GUE) or symplectic (GSE) ensemble respectively. While describing real many-body systems, time reversal symmetric Hamiltonians can be described by a GOE, whereas Hamiltonians with a broken time-reversal symmetry (for e.g. in a spin system with an external magnetic field) can be described by a GUE Mehta (2004). The joint probabitlity distribution of such matrices is,

[TABLE]

where $N$ is the rank of the matrix and the product in the measure is over all independent and identically distributed (i.i.d.) variables. However, the L.H.S. of (A) has a more generic meaning. It tells us that we can consider any ensemble of matrices that keeps the measure invariant under a similarity transformation,

[TABLE]

where $U$ could be an orthogonal or a unitary matrix. This implies that the random elements can be drawn from any ensemble that satisfies the above condition. The most general ensemble for the time-independent RMT can be written as,

[TABLE]

where $V(M)$ is the potential term. The ensemble corresponding to the Gaussian weights is $V(M)=\frac{1}{2}M^{2}$ . Let us work with this most general case and sketch out a solution in the large N limit along the lines of Brezin et al. (1978).

We start by diagonalizing the matrix via the above similarity transformation, $M=U^{-1}DU$. The ensemble, written in the basis of the eigenvalues of the matrix, looks as follows,

[TABLE]

where the action $S({\lambda_{i}})$ is,

[TABLE]

The first term in (42) is just the potential written as a function of the eigenvalues. The additional $N$ has appeared in front of this term as a result of the scaling of the eigenvalues by a factor of $\sqrt{N}$. This is done to ensure that the action doesn’t scale with any factor of $N$. As a result, all eigenvalues are now $\mathcal{O}(1)$ numbers. The second term is new and has appeared in the diagonalization process as a part of the measure. It is the Vandermonde determinant and it depends on the modulus of the difference between any pair of eigenvalues. It also depends on whether the matrices are orthogonal or unitary ($\beta=1$ or $2$ accordingly). We shall only concern ourselves with GUE here and hence put $\beta=2$ in the above expression. To find a solution, we need to extremise the action w.r.t. $\lambda_{i}$, which gives us the following equation,

[TABLE]

The above integral equation was solved by Brezin et al. (1978) using the method of resolvents that we shall describe below.

Appendix B The Method of Resolvents

Moving over to the continuum limit in the eigenvalues, we introduce a density of states $\rho(\lambda)$ containing information about the number of eigenvalues between $\lambda$ and $\lambda+d\lambda$. In the continuum limit, the saddle point equation (43) becomes,

[TABLE]

where $Pr$ refers to only the principal part of the integral.

The solution of this equation would give us the eigenvalue density (or level density) $\rho(\mu)$ which would tell us about the distribution of the eigenvalues in the large $N$ limit. To solve equation (44), we shall use the method of resolvents Brezin et al. (1978)(see Eynard et al. (2015) for a review)101010To find other ways to calculate the density of states, see Brezin and Hikami (1996b); Alexandrov et al. (2007) and references within Eynard et al. (2015).111111This method only helps us find the leading large $N$ results. To get to the subleading $1/N$ corrections, see Migdal (1983); David (1990); Kazakov (1991); Ambjorn et al. (1993). We thank Andrei Mironov for pointing out these references to us..

In the method of resolvents, we define the Stieltjes transform $\omega(\lambda)$ for the function $\rho(\lambda)$ as,

[TABLE]

For random matrices corresponding to real physical Hamiltonians, the eigenvalues are all real numbers. Thus, $\mu\in\mathbb{R}$. Outside the support of $\rho(\lambda)$, the function $\omega(\lambda)$ has the folllowing asymptotic behavior121212Assuming, $\int_{\text{supp}\mu}d\mu\,\rho(\mu)=1$,

[TABLE]

The Stieltjes transform extends the domain of the level density to complex values. Thus, $\lambda\in\mathbb{C}$. It is then easy to check that,

[TABLE]

The property of the resolvent is also that,

[TABLE]

which means that it has a jump on the real line along the support of $\rho(\lambda)$.

Combining the above two properties, the resolvent assumes the following form,

[TABLE]

The most general solution is,

[TABLE]

where $M(\lambda)$ is some polynomial in $\lambda$ and $\sigma(\lambda)$ is a polynomial in $\lambda$ of the following form,

[TABLE]

with $n$ being the no. of intervals on which $\rho(\lambda)$ is supported and $(a_{2i-1},a_{2i})$ being the end-points of the intervals.

Appendix C Level density for Polynomial Potentials

In case of the quartic potential, Brezin et al. (1978) provided us with a solution to the integral equation using the resolvent method. We wish to review that method here. The conditions satisfied by the resolvent are as follows,

it is analytic in the entire complex plane $\mathbb{C}$ except on the interval $(-2a,2a)$ along the real line. 2. 2.

it is real everywhere outside the interval $(-2a,2a)$. 3. 3.

it goes as $1/\lambda$ as $\lambda\rightarrow\infty$, with the eigenvalue density appropriately normalised. 4. 4.

Along the interval $(-2a,2a)$, the resolvent has a discontinuity of the following form,

[TABLE] 5. 5.

From its definition, the resolvent satisfies the following property,

[TABLE]

The last two properties together suggest the following form for the resolvent,

[TABLE]

The solution for $\rho(\lambda)$ must conform to the analyticity properties for the resolvent stated above. In fact, the unique form of the function that does so is,

[TABLE]

where $M$ and $\sigma$ are polynomials in $\lambda$, with $\sigma(\lambda)$ of the following form,

[TABLE]

where $n$ is the no. of intervals on which $\rho(\lambda)$ is supported and $(a_{2i-1},a_{2i})$ are the end-points of the interval.

Suppose we have an even polynomial potential of degree $deg[V]=P$. Let $n$ be the no. of intervals on the real line that form the support of $\rho(\lambda)$. Then, for a solution of $\rho(\lambda)$, it can be seen that,

[TABLE]

Quartic Potential

For $P=4$ and $n=1$, this tells us that,

[TABLE]

For $n=2$,

[TABLE]

However, the solution with two cuts in the quartic model is not physical, since it has complex solutions below $g=g_{c}$.

Thus, the one-cut solution is the only one that can be determined in the quartic model.

We assume the following form for $\rho(\lambda)$ in accord with (56),

[TABLE]

Putting this into (52) (on the UHP) and expanding around $\lambda\rightarrow\infty$, we find

[TABLE]

with the constraint,

[TABLE]

The form of the one-cut level density in the quartic case is,

[TABLE]

Sextic Potential

For the sextic potential, $P=6$. Now we have three different possibilities corresponding to the one, two and three-cuts.

For $n=1$,

[TABLE]

For $n=2$,

[TABLE]

For $n=3$,

[TABLE]

For simplicity, we shall only display the method for finding the one-cut solution. The two and three cut solutions follow analogously.

The unique form of the one-cut solution in accord with (62) is,

[TABLE]

Putting this form into the resolvent for the sextic model given by (12) and expanding around $\lambda\rightarrow\infty$, we find that

[TABLE]

with the constraint,

[TABLE]

Thus, the unique form of the one-cut level density in the sextic case is,

[TABLE]

Appendix D NNSD in Non-Gaussian Models

The nearest neighbour spacing distribution (NNSD) has been the hallmark of quantum chaos for a wide range of quantum mechanical systems. It is easily calculated for a large ensemble of two-level systems, where the different copies provide instances of neighbouring eigenvalues with varying differences. In the eigenvalue basis, it can be written as,

[TABLE]

where $\omega$ is the difference in the neighbouring eigenvalues. With the delta function, the integrals are easily evaluated to give,

[TABLE]

where $A$ and $B$ are related by the normalization condition: $\int_{0}^{\infty}P(\omega)=1$. The function $P(\omega)\sim\omega^{2}$ near $\omega\rightarrow 0$ and has a Gaussian tail for $\omega\rightarrow\infty$. We wish to find similar expressions for the non-Gaussian ensembles. For that, we introduce the quartic terms in (69) and define the new NNSD as,

[TABLE]

Integrating this, we get,

[TABLE]

where $K_{\frac{1}{4}}(x)$ is the modified Bessel function of the second kind. To understand its relevance, one must look at Fig.8 . Analytically, one can understand the distribution better by looking at the two limits of $\omega$,

[TABLE]

Near $\omega\rightarrow 0$, the distribution behaves analogous to the NNSD from the Gaussian ensemble, $P(\omega)\sim\omega^{2}$. However, in the other limit ($\omega\rightarrow\infty$), it has a non-Gaussian tail in contrast with its Gaussian counterpart. So, for $g\neq 0$, we see a deviation of $P(\omega)$ from the Wigner-Dyson distribution131313Note that this function has no special behaviour at $g_{c}=-1/48$.. However, from the form of $P(\omega)$ in Fig.8, added to the fact that the SFF diplays a universal behavior due to the universal behavior of the short distance correlations, one might speculate that the NNSD structure is still such that the system is chaotic. Our speculation is also partly based on an RG analysis that was carried out in Brezin and Zee (1993). According to their argument, the universality of the short distance correlations as well as the chaotic nature of the system has much to do with the presence of a trivial stable fixed point at the origin in the space of couplings.