Spectral Form Factor in Non-Gaussian Random Matrix Theories

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 behaviour exactly at the critical points of each model. The dip-time of the SFF is estimated in each of these models. The change in the decay behaviour is explained by relating it to the change in the behvaiour of the mean level density at criticality, near the edges of its support, using the Paley-Wiener theorem. The late time behaviour 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.


Introduction
Quantum Chaos has attracted the attention of physicists 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. The first major steps in this direction were taken by F. Dyson in a series of papers ( [1][2][3][4]) in which he masterfully introduced Random Matrices to characterise chaotic quantum systems. According to this description, any arbitrarily complicated many-body Hamiltonian could be replaced by a matrix of random numbers which have been drawn from a Gaussian ensemble. One could create an ensemble of such matrices with the idea that all possible interactions between the energy levels of the many-body Hamiltonian (written in a fixed basis) would be accounted for by the various matrices in the ensemble. If the Hamiltonian were time-reversal symmetric, one would require the probability distribution to be invariant under orthogonal transformations; if not, they must be invariant under unitarity transformations. In the thermodynamic limit (N → ∞), the eigenvalue density of random matrices showed a universal behaviour given by the so called the Wigner's semi-circle formula. This result seemed to be applicable to a varied class of quantum systems displaying chaotic behaviour. However, chaos was also a hallmark of few-body Hamiltonians (N finite). A better diagnostic for chaotic quantum systems was devised in which one looked at the distribution of the nearest neighbour spacings (NNSD) of eigenvalues of the system. The system was termed chaotic if the distribution was of the Wigner-Dyson type, where (A β , B β ) are system dependent constants and β = 1(2) for GOE (GUE). Systems that have this distribution were said to be chaotic, irrespective of whether they have a classical counterpart.
In the past few years, a new diagnostic for quantum chaos in conformal field theories (CFTs) has been suggested [5][6][7][8] that uses out-of-time-ordered correlators (OTOC). The idea is to start with a Euclidean CFT in d-dimensions with a thermal circle β. The temperature of the system is given by β −1 . In this set-up, one considers two sets of operators V and W and places them on the thermal circle at fixed distances of β/4 3 . Keeping the W 's fixed, the V 's are analytically continued to Lorentzian time t along a complex time contour. The 4-pt. function of this OTOC W V (t)W V (t) β goes to zero at large times implying decaying overlap between the states |V W and |W V . This realises the classical intuition for sensitive dependence on initial condition at the quantum level and pitches it as a definiton for quantum chaos. One could also calculate the rate at which the correlator decays to zero and identify it with the Lyapunov exponent for the system. In [6], an upper bound for the Lyapunov exponent was proposed for large N holographic CFTs and the bound was shown to be saturated for CFTs dual to Black Holes in the AdS bulk implying them to be maximally chaotic.
On a different note, the Spectral Form Factor (SFF) was introduced long ago [9] and recently re-introduced [10-16] as a tool for probing the spectrum of quantum systems. The SFF is defined as, where Z(β) is the partition function of a quantum system and β is the inverse temperature. For β = 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 nearest neighbour spectral distribution (NNSD). In fact, the SFF when averaged over Gaussian Random Matrices, has a very particular behaviour 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. However, this way of diagnosing chaos is more along the first line of thought where the NNSD characterises the chaotic behaviour. This seems to have very little to do with the OTOC way of understanding chaos in CFTs. There is only one example where both these methods have been used to study chaos in the same system, that is the Sachdev-Ye-Kitaev (SYK) model. An effort to bridge this gap has been made in [17]. The SYK model was first shown to be chaotic in the OTOC way in a talk by Kitaev [5] and then by Maldacena and Stanford [18]. Later [12] calculated the SFF for the model and found that it displays the same late time behaviour as when it is averaged over a Gaussian random matrix ensemble. This equivalence suggests that both these models have the same structure for the NNSD, although no explicit calculation exists for the SYK to verify this claim.
In this paper, we consider non-Gaussian matrix models that have a richer phase structure than the Gaussian model 4 . We specifically work out the Quartic and the Sextic models. We average the SFF over each of these non-Gaussian ensembles and make the following observations, 1. The initial decay behaviour is the same as for the Gaussian ensemble, when the non-Gaussian ensemble is away from its critical points.
2. At the critical points, the decay behaviour of the SFF changes over to a faster decay, displaying a different power law behaviour. This gives an explicit example where the multi-crtical behaviour at equilibrium is extended to non-equilibrium dynamics.
3. The dip-time estimate changes with the change in the fall-off behaviour at criticality.
4. The late time behaviour of the SFF is still characterised by a linear growth followed by saturation.
The last observation is the most crucial and suggests one of the following consequences,

1.
Any even polynomial potential model has a characteristic eigenvalue spacing distribution given by the Wigner-Dyson distribution and hence is a candidate for describing chaotic behaviour in quantum systems.
2. The SFF is not a sensitive enough diagnostic for probing the nearest neighbour spectrum of quantum systems and hence displays a universal behaviour for all polynomial potential systems.
Like in the SYK model, in the absence of explicit calculations for the NNSD in these non-Gaussian models, it is hard to suggest either conclusion and shall have to await future endeavours on our part. This paper is organized as follows. We start by reviewing the basics of Random Matrix Theories in Section 2. In Section 3, we describe 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 densitites for these models. In Section 4, 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 behaviour of the disconnected part of the SFF exactly at the critical points of these models. In Section 5 we estimate the dip-time for these models. We also relate the change in the decay behaviour to the change in the form of the one-cut mean level density in these models. In Section 6, we talk about the implications of non-Gaussian Matrix Models on chaotic systems.

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 [21]. The joint probabitlity distribution of such matrices is, 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 (3) 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, 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, where V (M ) is the potential term. The ensemble corresponding to the Gaussian weights is V (M ) = 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 [22]. 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, where the action S(λ i ) is, The first term in (7) 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 √ N . This is done to ensure that the action doesn't scale with any factor of N . As a result, all eigenvalues are now O(1) numbers. The second term is new and has appeared as 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 (β = 1 or 2 accordingly). We shall only concern ourselves with GUE throughout the rest of the text and hence put β = 2 in the above expression. To find a solution, we need to extremise the action w.r.t. λ i , which gives us the following equation, The above integral equation was solved by [22] using the method of resolvents that we shall describe below.

The Method of Resolvents
Moving over to the continuum limit in the eigenvalues, we introduce a density of states ρ(λ) containing information about the number of eigenvalues between λ and λ + dλ. In the continuum limit, the saddle point equation (8) becomes, where P r refers to only the principal part of the integral. The solution of this equation would give us the eigenvalue density (or level density) ρ(µ) which would tell us about the distribution of the eigenvalues in the large N limit. To solve equation (9), we shall use the method of resolvents [22](see [23] for a review) 56 .
In the method of resolvents, we define the Stieltjes transform ω(λ) for the function ρ(λ) as, For random matrices corresponding to physical Hamiltonians, the eigenvalues are all real numbers. Thus, µ ∈ R. Outside the support of ρ(λ), the function ω(λ) has the folllowing asymptotic behaviour 7 , The Stieltjes transform extends the domain of the level density to complex values. Thus, λ ∈ C. It is then easy to check that, The property of the resolvent is also that, which means that it has a jump on the real line along the support of ρ(λ). Combining the above two properties, the resolvent assumes the following form, The way to obtain a solution for the level density from the above form has been described in Appendix A. The most general solution is, where M (λ) is some polynomial in λ and σ(λ) is a polynomial in λ of the following form, with n being the no. of intervals on which ρ(λ) is supported and (a 2i−1 , a 2i ) being the end-points of the intervals.

Matrix Models: Gaussian and Non-Gaussian
Having found the solution for the most general potential, 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 on chaotic behaviour, in particular the spectral form factor, which we shall define in the next section.

Gaussian Matrix Model
We begin by considering a Gaussian potential: V (M ) = 1 2 M 2 . The corresponding action in the continuum limit is, Extremising this action with respect to the level density ρ(λ) gives us, Taking a derivative of the overall equation w.r.t. λ gives (9), with V (λ) = λ. For a Gaussian potential in a symmetric interval − 1, 1 , the resolvent takes the following form, From the properties of the resolvent (stated in Appendix A), it follows that the level density takes the unique form, This is the Wigner's semi-circle law for Gaussian matrices (shown in Fig. 1) in the large N limit.

Non-Gaussian Matrix Models
Quartic Potential As the simplest example of a non-Gaussian matrix model, we consider a quartic potential of the form, where, in the second term, g is the coupling constant and the power of N has been adjusted to match the rest of the terms in the action. This model has been shown ( [22], [30], [31]) to have a critical point at g = g c , below which it is supposed to have two saddle point solutions. The point g c is the critical point separating the two phases of the model. It is our goal to study the behaviour of the spectral form factor in the one-cut phase and see how it changes at the critical point. It can be argued (Appendix A) that there is a unique resolvent corresponding to the saddle point equation in this phase. This unique form of the resolvent in the interval (−2a, 2a) is, with the constraint, The level density is found from (64) to be, The end-points of the interval (−2a, 2a) depend on the value of the coupling constant g through (23) as, It is easy to see from above that a 2 will have imaginary values for g ≤ −1/48 and hence, there are no solutions of (23) that exist below g c = −1/48. Till g = g c , the coupling constant has real values which correspond to the end points of the single interval. One shortcoming of this analysis is that the two-cut phase of the model (beyond g = g c ) cannot be reached from the one-cut phase. To remedy this, we need to go to a matrix model with a sextic potential [30].

Sextic Potential
The form of the potential in the sextic model is, The unique form of the resolvent in the one-cut phase is given by, with the constraint, The level density corresponding to (27) is, From (28) 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 again look for the points where the behaviour of the level density changes near the edges of the spectrum. This calculation has been carried out in Appendix B. In this model, there is a critical line separating the one-cut and two-cut phases that starts from (g, h) = (−1/48, 0) and ends at a tri-critical point. The values of g and h on this critical line is [30], For a 2 = 2, the above values reduce to the critical point of the quartic model. The tri-critical point appears at the following values of the couplings, In Sec. 4.1.2, we will find these critical points from the late time behaviour of the SFF in the sextic model.

Spectral Form Factor
The spectral form factor(SFF) was introduced as a probe for analysing the spectrum of quantum systems. It is defined in terms of the analytically continued partition function as, From the R.H.S. of the above expression, it is easy to see that at very high temperatures (β → 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 behaviour. This enables the SFF to play a crucial role in understanding the time dynamics of such quantum systems and also makes it a very useful tool for probing the discreteness in the energy spectrum of quantum systems. On the other hand, chaotic quantum systems have been found to satisfy the Wigner's surmise, which is a statement about the distribution of the nearest neighbour spacings of the energy levels. This makes the SFF a suitable observable for analysing chaos in quantum systems. As we learnt in the previous sections, a lot is known about the energy spectrum of Random Matrix Theories in the large N limit. When the SFF is averaged over an ensemble of random matrices, it has a very particular behaviour [12], [22]. In this light, we wish to understand the general dependence of the SFF on the NNSD of the RMTs and establish a connection with 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 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 quantum chaos. The quantity we wish to study is defined as, where D(λ) is the eigenvalue density and the integration is over the (unscaled) support of ρ(λ) in each case. For simplicity and clarity of the analysis, we wish to focus on the β = 0 (or T = ∞) window. In the large N limit, the two point function of the level densities in (33) can be divided into two parts; the connected G c part and the disconnected G dc part, which we define below.
, beyond which the plateau appears. The ramp and the plateau appear to be universal for all polynomial potentials V (λ) due to the universality of the short distance correlators of the level densities. The initial decay, however, changes at the critical points of these polynomial potentials.
We define the quantity in the large bracket 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, then the connected two point correlation function can be written as just δD(λ)δD(µ) . The connected part of SFF can then be re-written as, We wish to analyse the behaviour of these two parts separately in Random Matrix Models with non-Gaussian potentials. We shall consider the quartic and sextic potentials and analyse the behaviour of the SFF in these models 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 (24). 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.
where D(λ) counts the total number of eigenvalues between (λ, λ + dλ) and α ∼ O( √ N ). This is also the expression we mean to use for calculating the SFF (since Z(β + it) | β,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 O(1). The density we get from our action is normalized as, where the interval parameter a ∼ O(1). The two are related by a scaling of the eigenvalues in the form: λ = √ N σ. 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, 1. The intergration limits would correspond to a ∼ O(1).

The time shall be scaled accordingly as
Using the scaled mean level density we then proceed to calculate the SFF.

Quartic Model
To calculate the disconnected part of the SFF, we use the scaled mean level density for the one-cut potential (24) into the first term of (34).
The above expression 8 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 This suggests that it is only the Bessel function of the second kind that governs the fall off behaviour 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 than at any other value of g. To understand this better, we look at the large time asymptotics of the above expression and find, This clearly suggests that for all values other than the critical value, the power law assumes a τ −3/2 behaviour. However, at g = g c the coefficient of the τ −3/2 term vanishes indicating a crossover to a τ −5/2 power law decay behaviour. One way to explain this behaviour is to look at the form of the level density (ρ(λ)) at the critical point. Away from the critical point, the level density has the form (24) such that it goes to zero near the λ = 2a edge of the interval as (λ − 2a) 1/2 . This behaviour results in the τ −3 power law beahviour of the SFF. At the critical point g c , the SFF displays the following behaviour, This implies that at the λ = 2a edge, the level density now goes to zero as (λ − 2a) 3/2 . The change in the power law decay behaviour can then be completely attributed to this particular change in behaviour of the level density near the edge of the spectrum. We shall explain this relation using the Paley-Weiner theorem in Section 5.

Sextic Model
For the sextic model, the mean level density in the one-cut phase is given by (29). This model turns out to have a richer phase structure as compared to its lower degree cousins.
There are three possible phases in this model. The one and two cut phases are now separated by an entire critical line. This line ends at a tri-critical point where all the phases of the model meet. 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, This expression is subject to a constraint equation, It is easy to check that putting h = 0 in the above expression gives us back (38). To check our method for the prediction of the critical point in the quartic model, we shall try to predict the critical points of this model in a similar way. We shall look at the vanishing of the coefficients of the various powers of τ from the expression for the large time expansion of (42). The large time behaviour of the SFF is, (1 + 24a 2 g + 180a 4 h) τ 3/2 cos π 4 + 2at For arbitrary values of h and g, (44) will have a large time fall off behaviour governed by the τ −3/2 term. However, to see the transition to a faster decay we would expect the coefficient of τ −3/2 to vanish while satisfying (43). Solving the set of simultaneous equations for g and h, we find This, of course, is a one-parameter set of solutions. One has a solution for every value of which implies that ρ(λ) ∼ (λ − 2a) 3/2 at the λ = 2a edge 9 . Further, at the tri-critical point, the mean level density becomes This clearly implies that ρ(λ) ∼ (λ − 2a) 5/2 at the λ = 2a edge. This then is the change responsible for the change in the behaviour of the SFF from τ −5/2 to τ −7/2 .

Universality of short distance correlations
The connected part of the SFF is as defined in (35). 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 δρ(λ) around the mean-level densitites. A very standard result of RMT [21] is the exact form of this correlator near the centre of the spectrum of the eigenvalues, 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 should hold true for any polynomial potential measure (of single trace operators) from which the ensemble of matrices is chosen. This point was previously mentioned in [33]. In the context of the SYK model, it then suggests that all such polynomial potential models are equally good for describing its late time behaviour. A physical reason is that the various polynomial potentials significantly change only the eigenvalue spectrum near the edges of the distribution (Fig.2). The correlations of the level densitites near the centre of the spectrum (|λ| ∼ 0), however, are local observables that are oblivious to the change in the global structure of the eigenvalue density of the system and hence remain unchanged.

Fourier Transform
To get to the SFF from (49), we need to plug it into the expression for the connected part in (34) 10 . It is easy to see that with a change of variables from (λ, µ) to (E = λ+µ, ω = λ−µ), we can perform the integral. There are two parts to this integral, which we shall deal with separately.
1. The 1/N 2 part with the sine squared function that is responsible for the ramp. It is also the sub-dominant contribution.
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 (49) gives us the linear rise in time of the SFF, as mentioned earlier. We wish to show why this is the case. Performing the intergral in the new variables (E, ω), we notice that the integral over E just gives a trivial delta function. We shall choose to work near E = 0 since (49) is valid only in this region. The remaining integrand is a function of the variable ω which is just the difference of the eigenvalues. So we are left with, The above integral gives us, These expressions tell us about the very long time behaviour of the SFF at a fixed energy. At τ < 2πN , the SFF has a linear growth in time. At τ ≥ 2πN , this linear behaviour saturates abruptly to a plateau region. From a physical point of view, it is enlightening to understand how this behaviour arises. In (49), we can re-define the argument of the sine term in terms of, We would like to work in the limit where, The scaling parameter x in (51) can be made small or large (by changing the way ω goes to zero) depending on the region of the spectrum we wish to focus on.
1. For large x( 1), (sin(x)/x) → 0 in (49). Only the delta function term remains. This scaling focuses on the correlations between fluctuations that are far apart in the spectrum. For these terms, The vanishing of the sine term implies that these fluctuations do not contribute to the SFF. This behaviour is dubbed the spectral rigidity.
2. For small x( 1), (sin(x)/x) → 1 in (49). This scaling focuses on the correlations between fluctuations in the nearby energy levels. The integral in (50) gets its maximum contribution from near the ω = 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 behaviour of the SFF. Evidently, the ramp stops when we reach the ω corresponding to values lower than the smallest eigenvalue spacing in the system, i.e. ω < 1/N .

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 However, we should remember that τ = √ N t. Thus, in terms of the physical time parameter, which is the estimated dip-time for this model.

Non-Gaussian RMT
Quartic Model 1. 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 which is the same as in the Gaussian model. Hence, the estimated dip time here would also be t ∼ O(1)

Paley-Wiener Theorem
A direct relation between the fall-off behaviour of the SFF and the edge behaviour of the level density, at the critical points, can be established using the Paley-Wiener theorem (as pointed out in [15]). The statement of the theorem is as follows, For a function g(ζ) defined on a compact spatial support, its Fourier transform F (η) has a lower bound on the rate of decay given by where N is a rational number and γ N is a real constant. Let us prove that our initial fall of behaviour is in accordance with this theorem in each of the cases.

The Proof 11
To begin with, let us identify g(ζ) and F (η) in the above definition with the level density ρ(λ) and Z(β + it) respectively. Of course, we shall consider β = 0 for the proof. The expression for the spectral form factor is, Let us denote ρ(λ) =ρ(λ). Let us also multiply a (−it) n on either side of the above equation and then consider its modulus squared, On the R.H.S., all powers of t can be absorbed into the integrand as n derivatives of the exponential term w.r.t. λ. Then, Using integration by parts, all derivatives w.r.t. λ can be transferred from the exponential term toρ(λ). This process will generate surface terms containing derivatives ofρ(λ) less than n. We will only keep the term with n derivatives ofρ(λ) and shall explain the disappearence of the surface terms in a while. Thus, Using Cauchy-Schwarz inequality for integrals on the R.H.S., we get The second term can be explicitly shown to be bounded from above by (4a) 2 . Hence, Taking a square root on either side and considering only positive real roots, we get where 2 is the L 2 norm of an integrable function. However, one can show that the L 1 norm of an integrable function is always greater than or equal to its L 2 norm. In such a case, This result is what gives us a direct relation between the decay of the SFF and the edge effect of the mean level density. We shall consider each of the cases separately.

Gaussian
In case of a Gaussian potential, the mean level density is given by (20). Plugging this into (61), we find that the integral of the second derivative ofρ(λ) blows up at the edges. Hence, there is a bound that exists only for n = 1. This implies that the rate of decay of | Z(it) | is bounded below by a t −1 fall-off behaviour and that of the disconnected part for the SFF by a t −2 fall-off behaviour. From explicit calculations we see that the exact decay behaviour of SFF is governed by a t −3 fall-off behaviour which is in accordance with the bound. For higher order polynomial potentials, we shall see the bound transitions in accord with the changing decay behaviour of the SFF at the critical points.

Quartic Potential
For the quartic potential, the mean level density in the one-cut phase is given by (24). We would like to plug this into (61) and consider the two regions of interest separately .

Away from the critical point
Away from the critical point, the intergral of the second derivative ofρ(λ) (for n = 2) blows up again. This results in exactly the same fall-off behaviour for the SFF as in the Gaussian potential case with the bound being satisfied similarly.

At the critical point
At the critical point (g = −1/48, a = √ 2),ρ(λ) behaves as (41). This suggests that the edge behaviour of the level density transitions from (2 √ 2 − λ) 1/2 to (2 √ 2 − λ) 3/2 behaviour. Hence, it is now the integral of the third derivative (for n = 3) ofρ(λ) that blows up at the edges. This means that the bound on the fall-off behaviour for | Z(it) | changes to t −2 . Thus, the SFF gets a lower bound of t −4 . However, from explicit calculations, we see that the fall off behaviour of the SFF transitions to t −5 which is just about right to satisfy the bound.
The fact that the bound changes right at the critical point in accordance with the change in the behaviour of the mean level density displays a direct relation between the beahviour of the two. We shall elucidate the result for the sextic potential as well.

Sextic Potential
For the sextic potential, the mean level density in the one-cut phase is (29). We have three regions of interest in this particular model and we shall explain the behaviour of (61) in each of these separately.

Away from criticality
Away from criticality, it is again the second derivative ofρ(λ) that blows up. Hence, the SFF has the same t −3 fall-of behaviour as in the previous two cases, abiding by the same bound of t −2 obtained from (61).

On the critical line
On the critical line (45), the mean level density assumes the form specified in (47). The edge behaviour ofρ(λ) transitions from the ( causing the third derivative of the integral ofρ(λ) to blow up. The decay behaviour of | Z(it) | is then bounded by a t −2 fall-of behaviour from (61). From explicit calculations we see that | Z(it) | decays as t −5/2 which is within what the bound specifies. The behaviour of the disconnected part of SFF follows accordingly.

At the tri-critical point
At the tri-critical point (46), the mean level density is given by (48). The edge behaviour now transitions to (12 − λ) 5/2 . It is now the fourth derivative (for n = 4) of the integral ofρ(λ) that blows up near the edges. Hence, the bound on the decay of | Z(it) | is now t −3 . Again, from explicit calculations, we find the decay behvaiour to be t −7/2 which is just above the bound right at the tri-critical point. The same is true for the full SFF.
This change in the fall-off behaviour of the SFF exactly at criticality in direct correspondence with the change in the edge behaviour of the mean level density verifies the existence of a direct relation between the two quantities.

The Surface Terms
In trying to move all the derivatives from the exponential term to the density of states, one generates a number of intermediate total derivative terms, which we call the surface terms. All these terms depend on k ≤ m − 1 derivatives of the density of states, where m is the highest integer value of n for which the bound holds. This just implies that all these surface terms vanish when evaluated at the edges. As an example one may consider the quartic density at the critical point. The form of the density at the critical point is proportional to (2 √ 2 − λ) 3/2 . In this case, m = 2 and so the two surface terms generated depend onρ andρ respectively. However, since the power ofρ is 3/2, the one derivative still keeps the term in the numerator such that, along with theρ term, it vanishes at the edges. One can see that the same logic holds true for the vanishing of all the surface terms that are generated in each of the cases, and for every polynomial potential, making this a very general argument.

Implications on Quantum Chaos
The universality of the short distance correlations between the eigenvalues of large random matrices, near the centre of the spectrum, was found in [21,33]. In Section 4.2, we emphasized this result that was found to hold for all polynomial potentials V (λ). In this section, we wish to point out the implications of this result on chaotic quantum systems via the SFF.
The SFF has been developed as a tool for probing the nearest neighbour eigenvalue spacings in quantum systems. However, due to the universality of the short distance correlations between eigenvalues, the very late time behaviour of the SFF appears to be completely universal. Away from criticality, the complete SFF does not, in any way, distinguish between the Gaussian and non-Gaussian models. To place this observation in the context of chaotic systems, one must mention the SYK model. As mentioned in the introduction, the SYK model has been shown to be chaotic using the OTOC analysis [18]. Infact it saturates the bound for the Lyapunov exponent [6] rendering it to be maximally chaotic. Recently, the SFF was calculated for the SYK model [12], where it was shown to have a late time behaviour (later than the dip time), similar to the Gaussian Random Matrix Models. The initial decay behaviour, however, did not have an exact overlap. In this light, we wish to point out that our observation of the universality in the late time behaviour of the SFF for arbitrary non-Gaussian, polynomial potentials (as indicated in Sec. 4.2.1) suggests the following consequences, 1. 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 behaviour was also suggested in [34], where RMT with arbitrary polynomial potential was shown to describe the IR behaviour of QCD 12 .
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 behaviour in Appendix B.

Conclusions
We have studied the Spectral Form Factor (SFF) in non-Gaussian random matrix models. We found the initial decay behaviour 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 behaviour to the change in the behaviour of the mean eigenvalue density at the edges of its support, exactly at the critical points of the model. This behaviour was shown to be in accord with the Paley-Weiner theorem, which sets a lower bound on the decay rate of the Fourier transform of an analytic function defined on a finite support. A crucial observation was made for the late time behaviour of the SFF. Beyond the dip-time the SFF turned out to be universal for any polynomial potential. This behaviour was due to the universality of the short range correlations between the lowest eigenvalues of the matrix models. With regard to implications on chaotic quantum systems, we wish to point out that any of the random matrix models with a polynomial potential seem to be good enough to describe the chaotic behaviour in such models.
Due to this universal behaviour of the SFF, one might think that the nearest neighbour spacing distribution in all such models are well described by the Wigner's Surmise. We provided some numerical simulations to support this claim. However, it would be interesting to develop analytic techniques where one could explicitly calculate the NNSD of the non-Gaussian models and match them to the Gaussian behaviour. There are hints that this universal behaviour persists for non-polynomial potentials as well. One interesting form would be to include double-trace operators and calculate the late-time behaviour in such cases to check if the aforesaid universality still persists. One very important goal of this project was to see the connection between the two different definitions of quantum chaos. In this respect, it would be really interesting to see if one could find the (maximum) Lyapunov exponent of these models from the SFF. This would bridge the existing gap in the equivalence of the two definitions 13 .
In [35,36], 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 [37], 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 behaviour.

Acknowledgement
We would like to thank our supervisor, Gautam Mandal, for suggesting the problem to us and also for numerous discussions during the course of this work. We would like to extend our gratitude to Saumen Datta and Sambuddha Sanyal for helping us with the codes that were used for generating the plots. We would especially like to thank Geet Ghanshyam, Sounak Biswas and Sujoy Chakraborty for various illuminating discussions during the course of this work. 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. Most of all, we would like to thank the people of India for their constant and generous support towards research in String Theory.
The solution for ρ(λ) must conform to the analyticity properties for the resolvent stated above. In fact, the unique form of the function that does so is, where M and σ are polynomials in λ, with σ(λ) of the following form, where n is the no. of intervals on which ρ(λ) is supported and (a 2i−1 , a 2i ) are the end-points of the interval. However, the solution with two cuts has problems in the quartic model since below g = g c , we cannot find any real solutions for the intervals. Thus, the first solution, corresponding to the one-cut phase, is the only one that can be determined in the quartic model. We assume the following form for ρ(λ) in accord with (68), Putting this into (64) (on the UHP) and expanding around λ → ∞, we find (a 1 , a 2 ) = 2g, with the constraint, 12ga 4 + a 2 − 1 = 0 (72) The form of the one-cut level density in the quartic case is, ρ(λ) = 1 π 1 2 + 4a 2 g + 2gλ 2 √ 4a 2 − λ 2 (73)

Sextic Potential
For the sextic potential, P = 6. Now we have three different possibilities corresponding to the one-cut, two-cut and three-cut phases. For simplicity, we shall only display the method for finding the one-cut phase solution. The two and three cut phase solutions follow analogously.

Appendix B
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, where ω is the difference in the neighbouring eigenvalues. With the delta function, the integrals are easily evaluated to give, P (ω) = Aω 2 exp(−Bω 2 ) (82) where A and B are related by the normalization condition: ∞ 0 P (ω) = 1. The function P (ω) ∼ ω 2 near ω → 0 and has a Gaussian tail for ω → ∞. We wish to find similar expressions for the non-Gaussian ensembles. For that, we introduce the quartic terms in (81) and define the new NNSD as, P (ω) ∝ dλ 1 dλ 2 δ(ω − (λ 1 − λ 2 ))|λ 1 − λ 2 | 2 e − 1 2 (λ 2 1 +λ 2 2 )−g(λ 4 1 +λ 4 2 ) Integrating this, we get, (1 + gω 2 (2 + 7gω 2 )) ( where K1 4 (x) is the modified Bessel function of the second kind. To understand its relevance, one must look at Fig.6 . Analytically, one can understand the distribution better by looking at the two limits of ω, Near ω → 0, the distribution behaves analogous to the NNSD from the Gaussian ensemble, P (ω) ∼ ω 2 . However, in the other limit (ω → ∞), it has a non-Gaussian tail in contrast with its Gaussian counterpart. So, for g = 0, we see a deviation of P (ω) from the Wigner-Dyson distribution. However, from the form of P (ω) in Fig.6, added to the fact that the SFF diplays a universal behaviour due to the universal behaviour 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 [33]. 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.