Universality and Thouless energy in the supersymmetric Sachdev-Ye-Kitaev Model

We investigate the supersymmetric Sachdev-Ye-Kitaev (SYK) model, $N$ Majorana fermions with infinite range interactions in $0+1$ dimensions. We have found that, close to the ground state $E \approx 0$, discrete symmetries alter qualitatively the spectral properties with respect to the non-supersymmetric SYK model. The average spectral density at finite $N$, which we compute analytically and numerically, grows exponentially with $N$ for $E \approx 0$. However the chiral condensate, which is normalized with respect the total number of eigenvalues, vanishes in the thermodynamic limit. Slightly above $E \approx 0$, the spectral density grows exponential with the energy. Deep in the quantum regime, corresponding to the first $O(N)$ eigenvalues, the average spectral density is universal and well described by random matrix ensembles with chiral and superconducting discrete symmetries. The dynamics for $E \approx 0$ is investigated by level fluctuations. Also in this case we find excellent agreement with the prediction of chiral and superconducting random matrix ensembles for eigenvalues separations smaller than the Thouless energy, which seems to scale linearly with $N$. Deviations beyond the Thouless energy, which describes how ergodicity is approached, are universality characterized by a quadratic growth of the number variance. In the time domain, we have found analytically that the spectral form factor $g(t)$, obtained from the connected two-level correlation function of the unfolded spectrum, decays as $1/t^2$ for times shorter but comparable to the Thouless time with $g(0)$ related to the coefficient of the quadratic growth of the number variance. Our results provide further support that quantum black holes are ergodic and therefore can be classified by random matrix theory.


I. INTRODUCTION
Random matrix theory [1][2][3][4][5][6] is a powerful tool to explain universal features of complex quantum systems. In general, it is applicable for long times scales where the system has equilibrated and therefore its motion depends only on global symmetries, such as time reversal invariance and charge conjugation. Quantitative agreement with random matrix theory indeed have been reported in a variety of problems: from disordered and quantum chaotic systems in the limit of negligible localization effects [7] to the spectrum of highly excited nuclei and Lattice Quantum Chromodynamics (QCD) [8][9][10]. In the latter, it was found that both the density and level statistics of the low-lying eigenvalues of the QCD Dirac operator follow the predictions of random matrix theory where the choice of the ensemble, that labels the universality class, depends on the representation of the gauge group and the number and flavors. Despite its simplicity, random matrices can capture important dynamical features of realistic strongly correlated systems, for instance chiral symmetry breaking, a salient feature of non-perturbative QCD, related to the infrared limit of the spectral density [11] of the QCD Dirac operator.
Despite its success, the random matrix approach has obvious limitations as for sufficiently short times where the quantum dynamics is not universal. In certain cases, such as non-interacting disordered metals and quantum chaotic systems, where there is a good understanding of these corrections, it is possible to estimate the time or energy scale, usually called Thouless energy or time, where non-universal corrections are expected and to even get quantitative information on the dynamics in this intermediate scale by studying the corrections to the random matrix results [12]. For strongly interacting systems, the situation is less clear though several approaches have been developed to tackle this problem in the context of condensed matter [13], QCD [14] and nuclear physics [15][16][17][18][19].
In the latter, fermionic models with infinity range two-body interactions were introduced as a natural generalization of random matrix models with the expectation of a larger range of applicability while keeping some analytical control. More recently, research on this problem has received an important new impetus after the realization [20] that N fermions with infinite range two-body interactions, now called Sachdev-Ye- Kitaev (SYK) model [20][21][22], could be a toy model for holography, namely, the infrared limit be could dual to a certain quantum black hole in an anti-deSitter space in two bulk dimensions (AdS 2 ).
This raises several questions: do SYK models with different global symmetries still keep most of the features expected in a theory with gravity-dual such as quantum chaos, exponential increase of low energy excitations or a finite entropy at zero temperature?
What is the precise window of universality in which random matrix theory is applicable and how are corrections beyond this universality regime?
In this paper we aim to answer these questions. We study the spectral density, thermodynamic properties and level statistics in a supersymmetric SYK model by analytical and exact diagonalization techniques. Our main results are that this chiral SYK model still has all the expected features of a gravity-dual: finite entropy at zero temperature, exponential increases of low energy excitations and excellent agreement of level statistics with random matrix theory predictions for sufficiently long time or small energy separations. We also show that not only level statistics but also the spectral density in the infrared limit is universal and given by random matrix theory. We study quantitatively deviations from these universal results, observed for energies (times) larger (smaller) than the so called Thouless energy (time) [12]. The number variance grows quadratically for eigenvalue separations larger than the Thouless energy and the spectral form factor decays as a power-law for times shorter than the Thouless time.
The paper is organized as follows. In section II, we start with the definition of the model and then discuss analytical results for the spectral density and spectral resolvent of the supercharge and the thermodynamic properties of the Hamiltonian. In section III we compare numerical results for the microscopic spectral density with random matrix predictions. Level statistics and its comparison with random matrix theory, including deviations for (small) large (time) energies, are discussed in section IV. Concluding remarks are made in section V. In Appendix A we work out the large N limit of the Q-Hermite result for the spectral density while a simple analytical form for the resolvent of the SYK model is obtained in Appendix B.
where the supercharge is defined by, with γ i Majorana fermions defined by the following algebraic relation also verified by Euclidean Dirac γ matrices. The coupling J ijk is chosen to be a Gaussian random variable with probability distribution, where J sets the scale of the distribution. Unless specified otherwise J ≡ 1.
This Hamiltonian is supersymmetric with global symmetries that, depending on N , correspond to one of the chiral and superconducting random matrix ensembles [49,53]. In this paper we will study the thermodynamical properties of the Hamiltonian, but we will analyze the spectral properties of the supercharge for E ∼ 0 corresponding to the ground state and low energy excitations of the Hamiltonian.

B. The Spectral Density of the Supercharge
The average spectral density for even q, computed in Ref. [33] by an explicit evaluation of the moments, turned out to be given by the weight function of the Q-Hermite polynomials with Q = η where η(N, q) is a suppression factor related to the commutation of two products of q Majorana operators. A similar calculation can be carried out for odd q. The average density is still given by the same analytical expression but with a negative suppression parameter (for even q, η can also become negative for small values of N ), In terms of |η| the average density can be rewritten as Note that this expression is valid for both even and odd q. After a Poisson resummation, an explicit evaluation of the resulting integrals is possible (see Appendix A). In the large N limit, it simplifies to, provided that the absolute values of the energy is away from E 0 given by where with σ 0 the standard deviation of J. For even q, this asymptotic expression turned to be not only an excellent approximation of the Q-Hermite result Eq. (6) but also of the exact density obtained by numerical diagonalization of Eq. (2). Results depicted in Fig. 1 confirm that Eq. (7) is still an excellent approximation of the Q-Hermite results. However, the agreement with the spectral density from exact diagonalization for N = 34 is only qualitative. Since correction terms to the moments in Q-Hermite approximation Eq. (6) are of order q/N the excellent agreement for q even was unexpected.

C. The resolvent
In this section, we study the resolvent of the supercharge of the SYK model which is an alternative way to investigate its average spectral properties. It is defined by where the λ k are the eigenvalues of Q. The spectral density is given by the discontinuity of the resolvent across the real axis. Conversely, the resolvent Eq. (10) follows by integrating over the spectral density. Because the spectrum of Q in Eq. (2) is symmetric under λ k → −λ k , the resolvent on the imaginary axis is purely imaginary (so that iG(is) is real). According to the Banks-Casher relation [11], the resolvent in the s → 0 and thermodynamic limits is directly related to the spectral density at E = 0. In the context of QCD, a finite value in this limit signals the spontaneous breaking of chiral symmetry, one of its most relevant low energy features. In supersymmetric SYK model, it is unclear whether this interpretation is exactly applicable.
However, it is still of interest to study it as a finite spectral density at the origin is a strong indication of a highly entangled ground state. Moreover, we have derived a compact analytical expression based on the known [60,61] spectral density for Q-Hermite polynomials. We leave the details to the appendix B and state here the final result [60], In Fig. 2 we depict results for the resolvent as a function of s for different values of N . We find excellent agreement with the analytical prediction Eq. (12) except for small s, dominated by high moments, where we expect the Q-Hermite result to be less accurate. In order to explore numerically the limit s → 0, N → ∞ we carry a finite size scaling analysis in N for a small s ∼ 0.01. We fit the N dependence of the logarithm of the resolvent by a function a + bN where a and b are fitting parameters. The slope is equal to the nontrivial part of the entropy density. We find b ≈ −0.0695 which agrees with the analytical result −π 2 /(16q 2 ) to be discussed in the next subsection.

D. Thermodynamical properties
In this subsection we study the thermodynamical properties of the model in the large N and low temperature limit where the asymptotic average density Eq. (7) is expected to be a good approximation. We note that because the Hamiltonian H = Q 2 , where Q is the supercharge, the low temperature limit is controlled by states with positive energy close to E = 0. In contrast with even q where for large N , E 0 ∝ N , for q odd when η < 0, The region close to the origin for odd q case is therefore qualitative different from the even q case where the ground state energy scales with N as should be the case for a fermionic system.
The entropy at zero temperature s 0 follows from the the prefactor c N in Eq. (6).
Although this constant is known analytically, see Appendix A, its large N limit can simply be determined from the condition Evaluating the integral by a saddle point approximation and using that in this limit log |η| ∼ −2q 2 /N we find Since the Q-Hermite approximation only has a nontrivial large N limit when q 2 /N is kept fixed, it can only yield a 1/q 2 correction to the zero temperature entropy which is given by This result agrees up to order 1/q 2 with the exact result [49] of s 0 obtained by using path integral techniques. We note that s 0 is larger than the s 0 in the even q case, where the nontrivial contribution to the entropy is π 2 /4q 2 . This larger entropy, which suggests a highly entangled ground state, is likely a consequence of the special properties at E ≈ 0 due to the chiral symmetry of the model.
As a further confirmation of this, we calculate the low temperature limit of the free energy from the asymptotic spectral density Eq. (7). We show that it agrees with the large q limit of the free energy for odd q obtained in [49].
Since the Hamiltonian H = Q 2 and ρ asym (E) is the eigenvalue density of Q, the partition function is given by For large N , the integral can be evaluated by a saddle point approximation with saddle point equation given by The free energy is therefore given by, Substituting and using the saddle point equation the free energy can be rewritten as If we identify −E 2 0 log |η| = J this is exactly the saddle point equation obtained in [49]. Taking the large N limit, it gives the free energy obtained in [49] with the constant energy N J /4q 2 subtracted. In the low-temperature limit, the solution of the saddle point equation is given by v = π resulting in the free energy The conclusion of this analysis is that the Q-Hermite form of the spectral density which at large N only has a nontrivial limit when q 2 /N is kept fixed reproduces the leading 1/q 2 correction to the free energy. Contrary to the even q case, for odd q the low-temperature limit is determined by energies close to zero. In the following section, we focus on this region only.

III. UNIVERSAL MICROSCOPIC SPECTRAL DENSITY AND RESOLVENT
Depending on N , the supersymmetric SYK model [49,52] has additional chiral and discrete symmetries. The analysis of thermodynamical properties has revealed that the main differences with respect to the original SYK model is in the region close to the ground density ρ M (E) defined by with ∆ the level spacing near E = 0, is universal and given by random matrix theory. in Fig. 3, shows and excellent agreement with random matrix ensemble belonging to the same universality class. We note that this is a completely parameter-free comparison.
For the chiral symplectic ensemble, we observe a suppression of the oscillations for eigenvalue sufficiently far from the origin. This is consistent with results from other systems such as QCD [9] where dynamical features not present in the random matrix ensemble tend to soften the effect of discrete symmetries. The scale at which the oscillations disappear, known as the Thouless energy, increases with N , and although an estimate of the precise N dependence is hard, it seems that it scales linearly with N . We note the distribution of the first eigenvalue was already worked out in Ref. [53] where agreement with random matrix theory was also found.

ESTIMATION OF THE THOULESS ENERGY
We now proceed to the study of dynamical properties of the SYK model by studying level statistics beyond the Thouless energy. In the non-supersymetric SYK model, it was found [33,34] that in both the bulk of the spectrum and the region close to the ground state, spectral correlations for small eigenvalues separations are well described by random matrix theory. This is an indication that for sufficiently long times, of the order of the Heisenberg time, the SYK model is still quantum chaotic and has reached an ergodic state where the dynamics is universal since it depends only on the global symmetry of the system.
We focus again in the low energy region E ≈ 0 where we expect more differences with respect to the SYK model with q even. Since agreement with random matrix theory is expected for short-range spectral correlations, we investigate the number variance, a popular choice to characterize long range spectral correlations which will allow us a more systematic description of the corrections to the random matrix results. It is defined as the variance of the number of levelsÑ ( ) in an energy interval of width (in units of the mean level spacing): In order to proceed, we have to unfold the spectrum, namely, to rescale it so that the mean level spacing is energy independent. In addition we use the connected two-point function to study spectral correlations. The removal of the one point function contribution is crucial to relate results from level statistics to dynamical features of the system.
We employ two unfolding methods: a polynomial of fifth degree fitting for the whole averaged spectrum and the splines method where consecutive groups of several eigenvalues, at least of the order of the spectral window to be investigated, are fitted to a linear or quadratic polynomial. Results from both methods are similar so we stick to the latter which provides a slightly more accurate fitting for the eigenvalues closest to the origin.
In terms of the unfolded spectrum Ñ ( ) = , so Ñ ≡ L where L is the number of unfolded eigenvalues in the window of interest. The number variance can be expressed as the double integral of the connected two level correlation function whose analytical expression is known for all ensembles of random matrices [62]. However, it is often simpler to generate these curves from numerical exact diagonalization of random matrices, which we will do for comparison to the SYK model. The number variance for L < 10, starting from E = 0, is shown in Fig. 4 for different values of N , corresponding to different universality classes. In all cases we find excellent agreement with the random matrix prediction for the first few eigenvalues. As was expected, because the SYK model is quantum chaotic, the agreement extends to more eigenvalues as N increases. The point of departure from the random matrix prediction, usually called Thouless energy, seems to scale as N α with α ∼ 1 for all universality classes. This is consistent with previous results for the microscopic spectral density. However the value of α, and therefore of the Thouless energy, for the non-supersymmetric SYK models (q > 2 even) seems to be much larger (α ≈ 2) [37]. For q even, there is an analytical argument that indeed α = 2 [64]. A reason for that behavior is that, for smaller q, the Hamiltonian is sparser and therefore it takes more time, so less energy, to explore the full phase space available.
Beyond the Thouless energy, we observe a growth of the number variance much faster than the logarithmic growth predicted by random matrix theory. Our results, shown in where the (short-range) correction δρ 2 (λ 1 , λ 2 ) is chosen such that the exact sum-rule is satisfied dλ 1 ρ 2,c (λ 1 , λ 2 ) = 0.
Note that ρ RMT 2,c satisfies this sum-rule. The constant term in the two-point-function gives a quadratic contribution to the number variance. From Fig. 5 we find that the constant term in (26) is well approximated by The power law increase of the number variance is reminiscent of the result for weakly disordered metals ∼ L d/2 [12], where d is the spatial dimensionality of the system though it is not clear to us whether it is possible to push this analogy any further.

A. Spectral form factor
In a time representation, the Thouless energy, and the subsequent quadratic growth of the number variance, are interpreted as the minimum time after which the dynamics black curve represents the disconnected part of g(t), g d (t) = β 2 /(β 2 + t 2 ), the leading contribution to g(t) for short times, that must be subtracted from g(t) in order to extract the Thouless time t T .
Top-Right: Connected spectral form factor g c = g −g d . For sufficiently long times, we observe a dip (correlation hole) at t T followed by a ramp and eventual saturation, all typical features of quantum chaotic system and also predicted by random matrix theory. The peak for N = 32 and N = 34 is a known feature of quantum chaotic systems with symplectic symmetry. We have slightly shifted all curves by 0.0004 so that we can plot it in a log-scale. Bottom-Left: Zoom in of g c for short times.
The position of the correlation hole seems to roughly scale with N . Bottom-Right: the deviation from the random matrix theory result δg(t) for t < t T , before the dip is formed. Assuming fully chaotic eigenfunctions, δg(t) describes how the SYK model approaches equilibrium. We have found that the numerical result agrees well with the theoretical prediction Eq (32) (solid curve) which decays as t −2 , as t T is approached.
is universal and how this universal regime is approached respectively. In order to explore these time scales, we compute the spectral form factor [65][66][67][68] with the unfolded spectrum, with Z = i e iλ i t−βλ i t with λ i the unfolded eigenvalues and β the inverse temperature. We also remove the disconnected part g d (t) = β 2 /(β 2 + t 2 ) related to the one-point function.
Although not always employed in the literature [34], it is necessary to work with unfolded eigenvalues and connected two-point functions for a correct determination of the Thouless energy.
We illustrate this point in Fig. 6 where we compute the full spectral form factor, for β = 0.03, together with the connected part g c = g − g d . We observe the ramp and saturation expected from random matrix theory in both cases. We also observe a dip (correlation hole) for short times. Only for g c this is related to the Thouless time, a time scale the describes the limit of applicability of random matrix theory. The dip (correlation hole) occurs for N = 34 at t ∼ 0.1, about 1/10 of the Heisenberg time. This is is roughly consistent with the estimation of the Thouless energy from the number variance in Fig. 4 which was computed removing the double degeneracy (for g(t) we kept the double degeneracy) and also using the unfolded spectrum and the connected part of the two-level correlation function.
The scaling of the t T , see Fig. 6 (bottom left), seems to be linear with N , also consistent with the number variance estimation. For t < t T , the spectral form factor deviates from the random matrix theory prediction. Simple perturbation theory [66,67], valid for t t T , shows that Similarly, for t t T , the constant contribution c M in Eq. 26 controls δg(t) so valid for any t < t T . The agreement is excellent until times very close to t T .
As for the number variance, our results are reminiscent to those for non-interacting weakly disordered metals in more than two spatial dimensions where a power-law decay was observed [12,69]. However there are important differences. In our case, diffusion is in Fock space and it is unclear whether the techniques and ideas employed for a noninteracting problem can be translated to the interacting SYK model.

V. OUTLOOK, HOLOGRAPHIC INTERPRETATION AND CONCLUSIONS
We have studied spectral and thermodynamical properties of the supersymmetric SYK model (q odd). We have found that the microscopic spectral density corresponding to the O(N ) smallest eigenvalues is universal, and is well described by chiral or superconducting random matrix ensembles depending on the value of N . The average spectral density, computed analytically and numerically, for E ≈ 0 grows exponentially with N , though the chiral condensate, obtained from the spectral resolvent, which is normalized with respect to the total number of eigenvalues, vanishes in the thermodynamical limit.
For slightly larger energies, the average spectral density grows exponentially. Spectral correlations in the region E ≈ 0 also agree with the random matrix prediction for sufficiently short eigenvalue separations. For eigenvalue separation larger than the Thouless energy, which roughly scales with N , we observe deviations from random matrix theory which are characterized by a quadratic growth of the number variance. We have also found a simple analytical expression for the spectral form factor, valid for times shorter than the Thouless time, that provides useful information about how the SYK model approaches ergodicity.
One of the main motivations to study the SYK model is its relevance in holography.
The observation of quantum chaos for long times and the exponential growth of the low energy excitations together with a finite zero temperature entropy strongly suggest that this supersymmetric extension could also have a gravity dual. Moreover it also raises some interesting questions. For instance, it would be interesting to understand in more detail the interpretation in the quantum gravity dual of the observed universal microscopic spectral density or the existence of a larger entropy at zero temperature with respect to the non-supersymmetric case. In this Appendix we calculate the large N limit of the spectral density for odd q.
For negative η the spectral density can be written as With the normalization constant c N given by the spectral density is normalized to 2 N/2 . After a Poisson resummation, the expression can be rewritten as Both integrals can be evaluated analytically and 1 − cosh 2π|n| log |η| arcsin(E/E 0 ) n sinh(nπ 2 / log |η|) .
For large N away from the edge of the spectrum the leading contribution to the spectral density is dominated by the n = 0 term resulting in the asymptotic form for the spectral density ρ asym (E) = c N cosh π arcsin(E/E 0 ) log |η| exp 2 arcsin 2 (E/E 0 ) log |η| . (A9) The large N limit of the normalization constant is determined by dEρ asym (E)dE = 2 N/2 .
For large N the integral can be evaluated by a saddle point approximation. Using that log |η| ∼ −2q 2 /N in this limit we find which give exactly the leading order 1/q 2 correction to the zero temperature entropy [49]. The resolvent of a Hamiltonian H is defined by where N is the dimension of the Hilbert space. If |z/λ k | > 1 for all eigenvalues λ k of H it can be expressed in terms of the moments of H. For a Hamiltonian for which the odd moments vanish, such as the SYK model, we obtain This sum has a finite radius of convergence due to the fact that it has a cut along the support of the spectrum (for a Gaussian the radius of convergence is zero). However, it may be possible to analytically continue G(z) inside the radius of convergence. This can be done using the resummation The sum is convergent on the imaginary axis. To prove this identity we start from the binomial expansion Defining z as The second line of (B9) also converges when s is small.
In the Q-Hermite case the moments are given by we will show next that the result (B9) can be simplified to To prove this, we first shift s → s √ 1 − η which cancels against the prefactor of the Q- to rewrite (B9) as a double sum which can be rewritten as a sum over p ≡ k + j fixed and a sum over j which runs from 0 to p. By splitting the sum over k form −p to p in the expression for the Q-Hermite moments into a sum from 0 to p and a sum form −p to −1 we see that the two expressions for the resolvent are the same.
A scale can be re-introduced by the substitution iG(is) → 1 σ iG(is/σ).