Supersymmetric SYK model and random matrix theory
SSupersymmetric SYK model and random matrix theory
Tianlin Li a Junyu Liu b Yuan Xin c Yehao Zhou d a Department of Physics and Astronomy, University of Nebraska,1400 R Street, Lincoln, Nebraska 68588, USA b Walter Burke Institute for Theoretical Physics, California Institute of Technology,1200 east California Boulevard, Pasadena, California 91125, USA c Department of Physics, Boston University,590 Commonwealth Avenue, Boston, MA 02215, USA d Perimeter Institute for Theoretical Physics,31 Caroline Street North, Waterloo, Ontario N2L 2Y5, Canada
E-mail: [email protected] , [email protected] , [email protected] , [email protected] Abstract:
In this paper, we investigate the effect of supersymmetry on the symmetry clas-sification of random matrix theory ensembles. We mainly consider the random matrix be-haviors in the N = 1 supersymmetric generalization of Sachdev-Ye-Kitaev (SYK) model, atoy model for two-dimensional quantum black hole with supersymmetric constraint. Someanalytical arguments and numerical results are given to show that the statistics of the su-persymmetric SYK model could be interpreted as random matrix theory ensembles, with adifferent eight-fold classification from the original SYK model and some new features. Thetime-dependent evolution of the spectral form factor is also investigated, where predictionsfrom random matrix theory are governing the late time behavior of the chaotic hamiltonianwith supersymmetry. a r X i v : . [ h e p - t h ] J un ontents N = 1 supersymmetric extension 5 N = 1 supersymmetric classification 103.2.1 Supercharge in N = 1 SYK 113.2.2 Hamiltonian in N = 1 theory 12 Physical systems with some stochastic or chaotic properties have some randomness in thesetup of the fundamental hamiltonian, which could be effectively simulated in the context ofrandom matrix theory. When choosing an ensemble from random matrix theory for a chaotichamiltonian, we often need to consider the symmetries in the dynamics of the related physicalsystem. The choice of standard matrix ensembles from symmetries, historically comes fromthe invention of Dyson [1], which is called three-fold way when classifying Gaussian UnitaryEnsemble (GUE), the Gaussian Orthogonal Ensemble (GOE), and the Gaussian SymplecticEnsemble (GSE). For more general symmetry discussion of interaction systems, the Altland-Zirnbauer theory gives a more complete description as a ten-fold classification [2, 3]. In thepractical usage, one of the most celebrated works would be the classification of interactioninside topological insulators and topological phases in a ten-fold way [4, 5].– 1 –n the recent study, the rising interests of studies on Sachdev-Ye-Kitaev (SYK) model givesanother profound application in the random matrix theory classification. SYK model [6, 7]is a microscopic quantum hamiltonian with random Gaussian non-local couplings among ma-jonara fermions. As is maximally chaotic and nearly conformal, this model could be treatedas a holographic dual of quantum black hole with AdS horizon through the (near) AdS/CFTcorrespondence [8–17]. In the recent research people have also discussed several generaliza-tions of the SYK model [18–21], such as higher dimensional generalizations and supersymmet-ric constraints. Some other related issues and similar models are discussed in [22–50]. In therecent discussions, people have discovered that the SYK hamiltonian has a clear correspon-dence with the categories of the three fold standard Dyson ensembles, unitary, orthogonaland sympletic ensembles, in the random matrix theory [51–54]. In the recent work, [53, 54],it is understood that the time-dependent quantum dynamics of the temperature-dependentspectral form factor, namely, the combinations of partition functions with a special analyticcontinuation in SYK model, is computable in the late time by form factors in the randommatrix theory with the same analytic continuation, as a probe of the discrete nature of theenergy spectrum in a quantum black hole, and also a solid confirmation on the three-foldclassification [54].In the route towards Dyson’s classification, one only considers the set of simple unitary oranti-unitary operators as symmetries when commuting or anticomuting with the hamiltonian.An interesting question would be, what is the influence of supersymmetry, the symmetry be-tween fermions and bosons in the spectrum, in the classification of symmetry class?As is illuminated by research in the past, supersymmetry [55] has several crucial influencesin the study of disorder system and statistical physics [56], and could be emergent from con-densed matter theory models [57]. Originating from particle physics, supersymmetry willenlarge the global symmetry group in the theory, has fruitful algebras and strong mathemat-ical power used in several models in quantum mechanics and quantum field theory, and isextremely useful to simplify and clarify classical or quantized theories. In the recent studyof SYK model, the supersymmetric generalization for the original SYK has been discussedin detail in [21], which shows several different behaviors through supersymmetric extensions.This model might give some implications in the quantum gravity structure of black hole intwo dimension in a supersymmetric theory, and also a related conjecture in [54] for spectralform factor and correlation functions in super Yang-Mills theory.In order to explore the supersymmetric constraints on the random matrix theory classifi-cation, in this paper we will study the symmetry classification and random matrix behaviorof the N = 1 supersymmetric extension of SYK model by Fu-Gaiotto-Maldacena-Sachdev’sprescription [21]. The effect of supersymmetry in the symmetry classification could be sum-marized in the following aspects, – 2 – Supersymmetry will cause the hamiltonian to show a quadratic expression. Namely, wecould write H as the square of Q . This condition will greatly change the distributionof the eigenvalues. From random matrix language [58], if Q is a Gaussian randommatrix, then H should be in a Wishart-Laguerre random matrix, with the eigenvaluedistribution changing from Wigner’s semi-circle to the Marchenko-Pastur distribution.In another sense, the quadratic structure will fold the eigenvalues of Q and cause apositivity condition for all eigenvalues. Namely, if Q has the eigenvalue distributionthat eigenvalues come in pair with positive and negative signs, the squaring Q willcause larger degeneracies and a folded structure in eigenvalues of energy. More over,the coupling degree might be changed when considering Q instead of H . For instance,in the N = 1 extended SYK model, Q is a non-local three point coupling, which isnot even. This will change the previous classification in the hamiltonian based on therepresentation of Clifford algebra from mathematical point of view. • We find the Witten index or Witten parity operator ( − F , which is well-known as acriterion for supersymmetry breaking [55, 59–61], is crucial in classifying the symmetryclass for supercharge Q . Some evidence of this point also could be found in some othermodels or setups. For instance, Witten parity is the Klein operator which separates thebosonic and fermionic sectors in the N = 2 supersymmetric systems [62, 63]. [64] pro-vides a more nontrivial example, where the odd parity operators are used to move statesalong a chain of different fermion sectors. Reversely, in some systems where one candefine a graded algebra, Klein operator serves as a key factor in realizing supersymme-try, which is helpful in models of bosonization and higher spin theories, etc.[65–68]. Forexample, [67] constructs the bosonized Witten supersymmetric quantum mechanics byrealizing the Klein operator as a parity operator. [68] realize a Bose-Fermi transforma-tion with the help of the deformed Heisenberg algebra which involves a Klein operator.Another interesting application of Witten operator is [69], where the author argues thatincorporating the Witten operator is crucial in some computation in supersymmetricsystems with finite temperature. In the supersymmetric SYK model we are consider-ing, Witten parity and the anti-unitary operator together become a new anti-unitaryoperator, which will significantly enlarge the set of symmetries in the hamiltonian, andchange the eight-fold story for supercharge Q and hamiltonian H .These aspects will be investigated in a clearer and more detailed way in the paper.This paper will be organized as the following. In Section 2 we will review the model construc-tion and thermodynamics of SYK model and its supersymmetric extensions. In Section 3 wewill discuss the random matrix classification for models, especially supersymmetric extensionsof the SYK model. In Section 4 we will present our numerical confirmation for symmetryclassifications from the exact diagonalization, including the computation of the density ofstates and spectral form factors. In Section 5, we will arrive at a conclusion and discuss thedirections for future works. In the appendix, we will review some knowledge to make this– 3 –aper self-contained, including basics on Altland-Zirnbauer theory and a calculation on therandom matrix theory measure. In this paper, we will mostly focus on SYK models and their extensions. Thus before themain investigation, we will provide a simple introduction on the necessary knowledge of relatedmodels to be self-contained.
In this part, we will simply review the SYK model mainly following [12]. The SYK model isa microscopic model with some properties of quantum black hole. The hamiltonian is givenby H = (cid:88) i 275 (2.17)Moreover, the one-loop correction from Schwarzian action is different. As a result of super-symmetry constraint, the one-loop factor is ( βJ N =1 ) − / Z Sch ( β ) ∼ βJ N =1 ) / e Ns + cN/ β (2.18)which predicts a different behavior for the density of states ρ ( E ) ∼ EJ N =1 ) / e Ns +2 cNE (2.19) For the generic positive integer ˆ q we can also define the N = 1 supersymmetric extension with non-localinteraction of 2ˆ q − Q = i ˆ q − (cid:88) i
It is established that SYK model is classified by random matrix theory in that the ran-dom interacting SYK hamiltonian fall into one of the three standard Dyson ensembles inthe eight-fold way [51–54]. It is natural to believe that the supersymmetric extension canalso be described by random matrix theory. To sharpen the argument, we derive the exactcorrespondence between each SYK hamiltonian and some random matrix ensembles, in otherwords, the eight-fold rule for supersymmetric case. A priori, the supersymmetric SYK hamil-tonian should lead to a different random matrix theory description than the original case.Superficially, the original SYK theory and its supersymmetric cousin are different have twomajor differences, which have been also mentioned in the previous discussions. • The degeneracy of the two hamiltonian matrices are different. The degeneracy of su-persymmetric SYK model is also investigated by [21], which we derive again using somedifferent discussion in Section 3.2.2. The degeneracy space is enlarged by supersym-metry. Generally, the energy level distribution of random matrices is sensitive to thedegeneracy and is thus sensitive to the supersymmetric extension. • Another difference is the apparent positive semidefiniteness of the hamiltonian beingthe square of the supercharge. We will see later that the positive constraint leads to anew eigenvalue distribution different from those of Gaussian ensembles.Symmetry analysis is crucial in classifying the random matrix statistics of hamiltonian matri-ces. [51, 54] argue that the particle-hole symmetry operator determines the class of randommatrix theory statistics. The random matrix classification dictionary is determined by the de-generacy and the special relations required by having the symmetry. The systematic methodof random matrix classification is established as the Atland-Zirnbauer theory [2, 3], reviewedin appendix A. The anti-unitary operators play a central role in the classifications. TheAtland-Zirnbauer also applies to extended ensembles different from the three standard Dysonensembles, which we find useful in classifying the supersymmetric SYK theory. In Section 3.1we derive again the eight-fold way classification of original SYK hamiltonian using Atland-Zirnbauer theory and find unambiguously the matrix representations of hamiltonian in eachmod eight sectors. We notice that the matrix representation of hamiltonian takes block di-agonal form with each block being a random matrix from a certain ensemble. This blockdiagonal form is also found by [51] in a different version.Naively one would apply the same argument to the supersymmetric hamiltonian, since italso enjoys the particle-hole symmetry. But this is not the full picture. First, one need totake into account of hamiltonian being the square of the supercharge and is thus not totallyrandom. In Section 3.2.1 we argue that the supercharge Q has a random matrix descriptionwhich falls into one of the extended ensembles. Using the Atland-Zirnbauer theory on Q weobtain its matrix representation in block diagonal form and use it to determine the matrix– 7 –epresentation of the hamiltonian in Section 3.2.2. Second, in order to obtain the correctclassification one needs to consider the full set of symmetry operators. Apparently particle-hole is not enough since supersymmetry enlarges the SYK degeneracy space. We argue thatthe Witten index operator, ( − F , is crucial in the symmetry analysis of any system withsupersymmetry. Incorporating ( − F we obtain the full set of symmetry operators. Finally,the squaring operation, will change the properties of the random matrix theory distribution ofsupercharge Q , from Gaussian to Wishart-Laguerre. The quantum mechanics and statisticsin supersymmetric SYK models, based on the main investigation in this paper, might be anon-trivial and compelling example of supersymmetric symmetry class. Now we apply the Altland-Zirnbauer classification theory (see appendix A for some necessaryknowledges) to the original SYK model [51–54]. This is accomplished by finding the symmetryof the theory (and has been already discussed in other works, see [51, 54]). First, one canchange the majonara fermion operators to creation annihilation operators c α and ¯ c α by ψ α = c α + ¯ c α √ ψ α − = i ( c α − ¯ c α ) √ α = 1 , · · · , N d = N/ 2. The fermionic number operator F = (cid:80) α ¯ c α c α divides the totalHilbert space with two different charge parities. One can define the particle-hole operator P = K N d (cid:89) α =1 ( c α + ¯ c α ) (3.2)where K is the complex conjugate operator ( c α and ¯ c α are real). The operation of P onfermionic operators is given by P c α P = ηc α P ¯ c α P = η ¯ c α P ψ i P = ηψ i (3.3)where η = ( − [3 N d / − (3.4)From these commutation relation we can show that[ H, P ] = 0 (3.5)To compare with the Altland-Zirnbauer classification, we need to know the square of P andthis is done by direct calculation P = ( − [ N d / = +1 N mod 8 = 0+1 N mod 8 = 2 − N mod 8 = 4 − N mod 8 = 6 (3.6)– 8 –ow we discover that P can be treated as a T + operator and it completely determines theclass of the hamiltonian. Before we list the result, it should be mentioned that the degeneracyof hamiltonian can be seen from the properties of P : • N mod 8 = 2 or 6:The symmetry P exchanges the parity sector of a state, so there is a two-fold degeneracy.However, there is no further symmetries caused by P in each block, Thus it is given asa combination of two GUEs, where two copies of GUEs are degenerated. • N mod 8 = 4:The symmetry P is a parity-invariant mapping and P = − 1, so there is a two-folddegeneracy. There is no further independent symmetries. From Altland-Zirnbauertheory we know that in each parity block there is a GSE matrix. Also, where two copiesof GSEs are independent. • N mod 8 = 0:The symmetry P is a parity-invariant mapping and P = 1. There is no further sym-metries so the degeneracy is one. From Altland-Zirnbauer theory we know that in eachparity block there is a GOE matrix. Also, two copies of GOEs are independent.We summarize these information in the following table as a summary of SYK model, N mod 8 Deg. RMT Block Type Level stat.0 1 GOE (cid:32) A B (cid:33) A, B real symmetric R GOE2 2 GUE (cid:32) A 00 ¯ A (cid:33) A Hermitian C GUE4 2 GSE (cid:32) A B (cid:33) A, B Hermitian quaternion H GSE6 2 GUE (cid:32) A 00 ¯ A (cid:33) A Hermitian C GUEwhere the level statistics means some typical numerical evidence of random matrix, for in-stance, Wigner surmise, number variance, or ∆ statistics, etc. Although the SYK hamilto-nian can be decomposed as two different parity sectors, we can treat them as standard Dysonrandom matrix as a whole because these two sectors are either independent or degenerated(The only subtleties will be investigating the level statistics when considering two independentsectors, where two mixed sectors will show a many-body localized phase statistics instead ofa chaotic phase statistics, which has been discussed originally in [51].) In the following wewill also numerically test the random matrix behavior, and based on the numerical testingrange of N we can summarize the following table for practical usage. N 10 12 14 16 18 20 22 24 26 28Ensemble GUE GSE GUE GOE GUE GSE GUE GOE GUE GSE– 9 – .2 N = 1 supersymmetric classification Supersymmetry algebra is a Z -graded algebra, where states and operators are subdividedinto two distinct parity sectors. In such an algebra there may exist a Klein operator [70]which anti-commutes with any operators with odd parity and commutes with any operatorswith even parity. The Klein operator of supersymmetry algebra is naturally the Witten indexoperator.Witten index might plays a role in the symmetric structure and block decomposition inthe supersymmetric quantum mechanics. A simple example is [70], in N = 2 supersymmetryalgebra, Define W be the Witten operator. The Witten operator has eigenvalue ± H = H + ⊕ H − . (3.7)We can also define projection operators P ± = 1 / ± W ). In the parity representation theoperators take 2 × W = (cid:32) − (cid:33) , P + = (cid:32) (cid:33) , P − = (cid:32) (cid:33) . (3.8)Because of Q = 0 and { Q, W } = 0 the complex supercharges are necessarily of the form Q = (cid:32) A (cid:33) , Q † = (cid:32) A † (cid:33) , (3.9)which imply Q = 1 √ (cid:32) AA † (cid:33) , Q = i √ (cid:32) − AA † (cid:33) . (3.10)In the above equation, A takes H − → H + and its adjoint A † takes H + → H − . The super-symmetric hamiltonian becomes diagonal in this representation H = (cid:32) AA † A † A (cid:33) . (3.11)In this construction, the Hilbert space is divided by Witten parity operator. The hamiltonianis shown to take the block diagonal positive semidefinite form without even referring to theexplicit construction of the hamiltonian. It is remarkable that the above computation is verysimilar to our work from Section 3.2.1 to 3.2.2. Applications of this property can be foundin[62, 63]. They describe a supersymmetric Quantum Mechanics system where fermions scat-ter off domain walls. The supercharges are defined as a differential operator and its adjoint.From (3.11) the number of ground states of each Z sector is simply the kernel of the differ-ential operator and the Witten index is computed. A more non-trivial example is provided– 10 –y [64]. In this work, the Hilbert space is divided into an N fermions Fock space. Thus theHamiltonian can be expressed as the direct sum of all fermion sectors. The ladder operators Q and Q † are odd operators and move states between different sectors.The argument can also work in reversive way. Hidden supersymmetry can be found in abosonic system such as a Calogero-like model [71], a system of one dimensional Harmonicoscillators with inverse square interactions and extensions. What makes supersymmetry man-ifest is the Klein operator. The model and its various extensions are studied in [65–68, 72]. Atrivial simple Harmonic operator has algebra [ a − , a + ] = 1. The algebra describes a bosonicsystem. Z grading is realized by introducing an operator K = cos( πa + a − ). The new oper-ator anti-commutes with a − and a + thus is a Klein operator. Based on the Klein operatorone can construct the projection operators on both sectors and also the supercharge. In thisway the simple harmonic oscillator is “promoted” to have supersymmetry. A generalizationto simple hamornic oscillator is the deformed Heisenberg algebra, [ a − , a + ] = 1 + νK . Thecorresponding system is an N = 2 supersymmetric extension of the 2-body Calogero. Themodel is also used in considerably simplifying Calogero model.These evidences strongly support the argument that supersymmetry will change the clas-sification of symmetry class in quantum mechanical models. In the following work, we willshow that supersymmetric SYK model symmetry class can be explicitly constructed andchange the classification of random matrix theory ensembles. N = 1 SYK In the N = 1 supersymmetric model, it should be more convenient to consider the spectrumof Q instead of H , because H is the square of Q . Although Q is not a hamiltonian, since weonly care about its matrix type, and the Altland-Zirnbauer theory is purly mathematical, Q can be treated as a hamiltonian. Similiar to the original SYK model, we are concerned aboutthe symmtry of the theory. We notice that the Witten index ( − F is( − F = ( − i ) N d N (cid:89) i =1 ψ i = N d (cid:89) α =1 (1 − c α c α ) (3.12)which is the fermionic parity operator up to a sign ( − N d . Witten index and particle-holesymmetry have the following commutation relation: P ( − F = ( − N d ( − F P (3.13)Now we define a new operator, R = P ( − F . It has a compact form R = K N d (cid:89) α =1 ( c α − ¯ c α ) (3.14) R and P are both anti-unitary symmetries of Q , with commutation relations:– 11 – mod 8 P R P, Q ] = 0 { R, Q } = 02 { P, Q } = 0 [ R, Q ] = 04 [ P, Q ] = 0 { R, Q } = 06 { P, Q } = 0 [ R, Q ] = 0and squares P = ( − [ N d / , R = ( − [ N d / N d (3.15)Thus, in different values of N , the two operators P and R behave different and replace therole in T + and T − in the Altland-Zirnbauer theory. Now we can list the classification for thematrix ensemble of N = 1 supersymmetric SYK model N mod 8 T T − Λ Cartan Label Type0 P = 1 R = 1 1 BDI (chGOE) R R = − P = 1 1 DIII (BdG) H P = − R = − H R = 1 P = − R One can also write down the block representation of Q . Notice that the basis of block decom-position is based on the ± N = 1 theory Now we already obtain the random matrix type of the supercharge. Thus the structure ofthe square of Q could be considered case by case. Before that, we can notice one generalproperty, that unlike the GOE or GSE group in SYK, in the supersymmetric model there isa supercharge Q contains odd number of Dirac fermions as a symmetry of H , thus it alwayschanges the parity. Thus the spectrum of H is always decomposed to two degenerated blocks.Another general property is that the spectrum of H is always positive because Q is Hermitianand H = Q > 0. Thus the random matrix class of N = 1 will be some classes up to positivityconstraint. • N = 0 mod 8: In this case Q is a BDI (chGOE) matrix. Thus we can write down theblock decomposition as Q = (cid:32) AA T (cid:33) (3.16)where A is a real matrix. Thus the hamiltonian is obtained by H = (cid:32) AA T A T A (cid:33) (3.17)– 12 –ince AA T and A T A share the same eigenvalues ( { R, Q } = 0 thus R flips the signof eigenvalues of Q , but after squaring these two eigenvalues with opposite signaturesbecome the same), and there is no internal structure in A (in this case P is a symmetryof Q , [ P, Q ] = 0, but P = 1, thus P cannot provide any further degeneracy), we obtainthat H has a two-fold degeneracy. Moreover, because AA T and A T A are both realpositive-definite symmetric matrix without any further structure, it is nothing but thesubset of GOE symmetry class with positivity condition. These two sectors will beexactly degenerated. • N = 4 mod 8: In this case Q is a CII (chGSE) matrix. Thus we can write down theblock decomposition as Q = (cid:32) BB † (cid:33) (3.18)where B is a quaternion Hermitian matrix. Thus after squaring we obtain H = (cid:32) BB † B † B (cid:33) (3.19)Since BB † and B † B share the same eigenvalues, and each block has a natural two-folddegeneracy by the property of quaternion (Physically it is because { R, Q } = 0 thus R flips the sign of eigenvalues of Q , but after squaring these two eigenvalues with oppositesignatures become the same. Also, in this case P is a symmetry of Q , [ P, Q ] = 0, and P = − H . Because BB † and B † B are quaternion Hermitian matrices when B is quaternion Hermitian , BB † = B † B areboth quaternion Hermitian positive-definite matrix without any further structure. Asa result, it is nothing but the subset of GSE symmetry class with positivity condition.These two sectors will be exactly degenerated. • N = 2 mod 8: In this case Q is a DIII (BdG) matrix. Thus we can write down theblock decomposition as Q = (cid:32) Y − ¯ Y (cid:33) (3.20)where Y is a complex, skew-symmetric matrix. Thus after squaring we obtain H = (cid:32) − Y ¯ Y − ¯ Y Y (cid:33) (3.21)Firstly let us take a look at the degeneracy. Since − Y ¯ Y and − ¯ Y Y share the same eigen-values and each block has a natural two-fold degeneracy because in skew-symmetricmatrix the eigenvalues come in pair and after squaring pairs coincide (Physically it is– 13 –ecause { P, Q } = 0 thus P flips the sign of eigenvalues of Q , but after squaring thesetwo eigenvalues with opposite signatures become the same. Also, in this case R is asymmetry of Q , [ R, Q ] = 0, and R = − H .Now take the operator Q as a whole, from the previous discussion, we may note that itis quaternion Hermitian because it could be easily verified that Q Ω = Ω Q and Q † = Q .Thus Q = H must be a quaternion Hermitian matrix (there is another way to see that,which is taking the block decomposition by another definition of quarternion Hermitian,squaring it and check the definition again). Moreover, H has a two-fold degenerated par-ity decomposition thus in each part it is also a quarternion Hermitian matrix. Becausein the total matrix it is a subset of GSE symmetry class (with positivity constraint), ineach degenerated parity sector it is also in a subset of positive definite GSE symmetryclass (one can see this by applying the total measure in the two different, degeneratedpart). • N = 6 mod 8: In this case Q is a CI (BdG) matrix. Thus we can write down the blockdecomposition as Q = (cid:32) Z ¯ Z (cid:33) (3.22)where Z is a complex symmetric matrix. Thus after squaring we obtain H = (cid:32) Z ¯ Z 00 ¯ ZZ (cid:33) (3.23)Since Z ¯ Z and ¯ ZZ share the same eigenvalues ( { P, Q } = 0 thus P flips the sign of eigen-values of Q , but after squaring these two eigenvalues with opposite signatures becomethe same), and there is no internal structure in Z (in this case R is a symmetry of Q ,[ R, Q ] = 0, but R = 1, thus R cannot provide any further degeneracy), we obtain that H has a two-fold degeneracy.Similar with the previous N mod 8 = 2 case, we can take the operator Q and H asthe whole matrices instead of blocks. For H we notice that the transposing operationsmake the exchange of these two sectors. However, the symmetric matrix statement isbasis-dependent. Formally, similar with the quarternion Hermitian case, we can extendthe definition of symmetric matrix by the following. DefineΩ (cid:48) = (cid:32) (cid:33) (3.24)and we could see that a matrix M is symmetric real (or symmetric Hermitian) if andonly if M † = M and M T Ω (cid:48) = Ω (cid:48) M (where Ω (cid:48) means the basis changing over two– 14 –ectors). We can check easily that Q satisfies this condition, thus Q = H must satisfy.Thus we conclude that the total matrix H in a subset of GOE symmetry class (withpositivity constraint).Although from symmetric point of view, the hamiltonian of N = 1 model should be classifiedin the subsets of standard Dyson ensembles. But what the subset exactly is? In fact, the spe-cial structure of the squaring from Q to H will change the distribution of the eigenvalues fromGaussian to Wishart-Laguerre [58, 73, 74] (Although there are some differences in the powersof terms in the eigenvalue distributions.) We will roughly called them as LOE/LUE/LSE, ashas been used in the random matrix theory research. Some more details will be summarizedin the appendix B.However, the difference in the details of the distribution, beyond numerical tests of the dis-tribution function of the one point-eigenvalues, will not be important in some physical tests,such as spectral form factors and level statistics (eg. Wigner surmise). The reason could begiven as follows. From the supercharge point of view, because Q is in the Altland-Zirnbauerdistribution with non-trivial ˜ α (see appendix B), the squaring operation will not change thelevel statistics such as Wigner surmise and spectral form factors (which could also be verifiedby numerics later). From the physical point, as is explained in [51], the details of distribution(even if not Gaussian), cannot change the universal properties of symmetries.Finally, we can summarize these statements in the following classification table (the degen-eracies have been already calculated in [21]), N mod 8 Deg. RMT Block Type Level stat.0 2 LOE (cid:32) AA T A T A (cid:33) A real R GOE2 4 LSE (cid:32) − Y ¯ Y − ¯ Y Y (cid:33) Y complex skew-symmetric H GSE4 4 LSE (cid:32) BB † B † B (cid:33) B Hermitian quaternion H GSE6 2 LOE (cid:32) Z ¯ Z 00 ¯ ZZ (cid:33) Z complex symmetric R GOEFor our further practical computational usage, we may summarize the following table fordifferent N s in the supersymmetric SYK random matrix correspondence. As we show in thenext section, for N ≥ 14, these theoretical consideration perfectly fits the level statistics. N 10 12 14 16 18 20 22 24 26 28RMT LSE LSE LOE LOE LSE LSE LOE LOE LSE LSEUniversal Stat. GSE GSE GOE GOE GSE GSE GOE GOE GSE GSE– 15 – Exact Diagonalization In this part, we will present the main results from numerics to test the random matrix theoryclassification in the previous investigations. One can diagonalize the hamiltonian exactly withthe representation of the Clifford algebra by the following. For operators acting on N d = N/ γ ζ − = 1 √ N d − (cid:89) p =1 σ zp σ xN d γ ζ = 1 √ N d − (cid:89) p =1 σ zp σ yN d (4.1)where σ p means standard Pauli matrices acting on the p -th qubit, tensor producting theidentity matrix on the other parts, and ζ = 1 , , ......, N d . This construction is a representationof the Clifford algebra { γ a , γ b } = δ ab (4.2)And one can exactly diagonalize the hamiltonian by replacing the majonara fermions withgamma matrices to find the energy eigenvalues. Thus, all quantities are computable by bruteforce in the energy eigenstate basis.The main results of the following investigation would be the following. In the density ofsupercharge eigenstates and energy eigenstates in the supersymmetric SYK model, the be-havior is quite different, but coincides with our estimations from the random matrix theoryclassification: the spectral density of supercharge Q shows clearly the information about ex-tended ensembles from Altland-Zirnbauer theory, and the spectral density of energy H showsa clear Marchenko-Pastur distribution from the statistics of Wishart-Laguerre. Moreover,because both Q and H both belongs to the universal level statistical class for GOE, GUEand GSE, the numerics from Wigner surmise and spectral form factor will show directly theseeight-fold features. We say a matrix M is a quaternion Hermitian matrix if and only if M = (cid:32) A + iB C + iD − C + iD A − iB (cid:33) for some real A, B, C, D in a basis, and A is symmetric while B, C, D is skew-symmetric. There is an equivalentdefinition that, defining Ω = (cid:32) − (cid:33) thus M is a quaternion Hermitian matrix if and only if M † = M and M Ω = Ω M . Thus it is shown directlythat if M is quaternion Hermitian then ( MM † ) † = MM † and MM † Ω = M ( M Ω) = M Ω M = Ω M = Ω MM † ,thus MM † = M = M † M is still a quaternion Hermitian matrix. – 16 – / N N ρ ( E ) / N N ρ ( E ) - - / N N ρ ( E ) Figure 1 . The density of states for original SYK model Hamiltonian (left), supersymmetric SYKHamiltonian (middle) and SUSY SYK supercharge operators treated as Hamiltonian (right) by exactdiagonalization. Density of states from N = 10 to N = 28 are plotted in colors from light blue todark blue. The eigenvalues have been rescaled by E ( Q ) /N J while the density of states has been alsorescaled to match the normalization that the integration should be 1. The plots for density of states in SYK model and its supersymmetric extension are shownin Figure 1 for comparison. For each realization of random hamiltonian, we compute alleigenvalues. After collecting large number of samples one can plot the histograms for allsamples as the function ρ ( E ). For density of states in SYK model, in small N tiny vibrationsare contained, while in the large N the distribution will converge to a Gaussian distributionbesides the small tails. However, in the supersymmetric SYK model the energy eigenvaluestructure is totally different. All energy eigenvalues are larger than zero because H = Q > N , andthe figure will come to a convergent distribution. The shape of this distribution matches theeigenvalue distribution of Wishart-Laguerre, which is the Marchenko-Pastur distribution [75]in the large N limit. For the supercharge matrices, as N becomes larger the curve acquiresa dip at zero, which is a clear feature for extended ensembles and could match the averageddensity of eigenvalues of random matrices in CI, DIII [3] and chiral [76] ensembles at large N .For numerical details, we compute N = 10 (40000 samples), N = 12 (25600 samples), N = 14(12800 samples), N = 16 (6400 samples), N = 18 (3200 samples), N = 20 (1600 samples), N = 22 (800 samples), N = 24 (400 samples), N = 26 (200 samples), and N = 28 (100samples). The results for original SYK model perfectly match the density of states obtainedin previous works (eg. [12, 54]). There exists a practical way to test if random matrices from a theory are from some specificensembles. For a random realization of the hamiltonian, we have a collection of energyeigenvalues E n . If we arrange them in ascending order E n < E n +1 , we define, ∆ E n = E n − E n − to be the level spacing, and we compute the ratio for the nearest neighbourhood spacingas r n = ∆ E n / ∆ E n +1 . For matrices from the standard Dyson ensemble, the distribution– 17 – igure 2 . The theoretical Wigner surmises for three different standard ensembles. The lower (blue),middle (red) and higher (green) curves are corresponding to GOE, GUE and GSE universal classrespectively. of level spacing ratio satisfies the Wigner-Dyson statistics[77]) (which is called the Wignersurmise p ( r ) = 1 Z ( r + r ) ˜ β (1 + r + r ) β/ (4.3)for GOE universal class, ˜ β = 1, Z = 8 / 27; for GUE universal class, ˜ β = 2, Z = 4 π/ (81 √ β = 4, Z = 4 π/ (729 √ 3) (In fact, these are labels for the field ofrepresentation. See appendices for more details). Practically we often change r to log r , andthe new distribution after the transformation is P (log r ) = rp ( r ). Standard Wigner surmisesare shown in the Figure 2. [51] has computed the nearest-neighbor level spacing distributionof the SYK model, which perfectly matches the prediction from the eight-fold classification.What is the story for the N = 1 supersymmetric SYK model? A numerical investigationshows a different correspondence for the eight-fold classification, which is given by Figure 3.One can clearly see the new correspondence in the eight-fold classification for supersymmetricSYK models, as has been predicted in the previous discussions.Some comments should be given in this prediction. Firstly, one have some subtleties in obtain-ing correct r s. Considering there are two different parities in the SYK hamiltonian ( F mod 2),each group of parity should only appear once in the statistics of r n . For N mod 8 = 0 , P maps each sector to itself, thus if we take all r n the distribu-tion will be ruined, serving as a many-body-localized distribution (the Poisson distribution).For N mod 8 = 2 , P maps even and odd parities toeach other, and one can take all possible r s in the distribution because all fermionic paritysectors are degenerated. Similar things are observed for all even N in the supersymmetric– 18 – Level StatN = 10 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 12 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 14 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 16 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 18 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 20 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 22 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 24 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 26 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) H Level StatN = 28 H MatrixGUEGOEGSE - P ( l og ( r )) log ( r ) Figure 3 . The nearest-neighbor level spacing distribution forhamiltonian of N = 1 supersymmetric SYK model for different N . The lower (blue), middle (red) and higher (green) curves aretheoretical predictions of Wigner surmises from GOE, GUE andGSE respectively. The black dashed curves are distributions forall r s from a large number of samples. SYK model. As we mentioned before, the reason is that the supercharge Q is a symmetryof H , which always changes the particle number because it is an odd-point coupling term.Moreover, the standard ensemble behavior is only observed for N ≥ 14, and for small enough N s we have no clear correspondence. Similar things happen for original SYK model, wherethe correspondence works only for N ≥ 5, because there is no thermalization if N is toosmall [51]. However, the threshold for obtaining a standard random matrix from N = 1supersymmetric extension is much larger.In Section 3.2.1, we argued that the supercharge operator Q in N = 1 supersymmetric SYKtheory are also random matrices in some extended ensembles [2, 3]. We compute the levelstatistics of Q and compare it with the Wigner surmises of three standard Dyson ensemblesin cases with different N . The result is presented in Figure 4. We see the level statistics of Q matrices match the same ensembles as the corresponding hamiltonian. This result confirmsthe relationship between Q ’s random matrix ensemble and that of the corresponding H . Thatwe do not see extended ensemble in the Q ’s level statistics because the level statistic does not– 19 – Level StatN = 10 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 12 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 14 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 16 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 18 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 20 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 22 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 24 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 26 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Q Level StatN = 28 Q MatrixGUEGOEGSE - - P ( l og ( r )) log ( r ) Figure 4 . The nearest-neighbour level spacing distributionfor the supercharge matrix Q of N = 1 supersymmetric SYKmodel for different N . The lower (blue), middle (red) and higher(green) curves are theoretical prediction of Wigner surmisesfrom GOE, GUE and GSE, respectively. The black dashedcurves are distributions for all r s from a large number of sam-ples. see all the information in the ensembles. Before presenting the numeric results of spectral form factors, we will review the discretenessof spectrum and the spectral form factor following [54]. For a quantum mechanical system,the partition function Z ( β ) = Tr( e − βH ) (4.4)could be continued as Z ( β, t ) = Z ( β + it ) = Tr( e − βH − iHt ) (4.5)The analytically continued partition function Z ( β, t ) is an important quantity to understanda discrete energy spectrum. Typically, people will compute the time average to understandthe late time behavior, but for Z ( β, t ), it vibrates near zero at late time and the time average– 20 – - - - - - - - t g ( t ) - - - - - - - t g c ( t ) - - - - - - - - - t g d ( t ) - - - - - - - t g ( t ) - - - - - - - t g c ( t ) - - - - - - - - - t g d ( t ) - - - - - - - t g ( t ) - - - - - - - t g c ( t ) - - - - - - - - - t g d ( t ) N = 10 N = 12 N = 14 N = 16 N = 18 N = 20 N = 22 N = 24 N = 26 N = Figure 5 . The spectral form factors g ( t ), g c ( t ) and g d ( t ) in the supersymmetric SYK model with J N =1 = 1, β = 0 , , 10 respectively. should be zero. Thus, we often compute (cid:12)(cid:12)(cid:12) Z ( β,t ) Z ( β ) (cid:12)(cid:12)(cid:12) . For a discrete energy eigenvalue spectrum,we have (cid:12)(cid:12)(cid:12)(cid:12) Z ( β, t ) Z ( β ) (cid:12)(cid:12)(cid:12)(cid:12) = 1 Z ( β ) (cid:88) m,n e − β ( E m + E n ) e i ( E m − E n ) t (4.6)It’s hard to say anything general directly for a general spectrum, but one can use the long-termaverage 1 T (cid:90) T (cid:12)(cid:12)(cid:12)(cid:12) Z ( β, t ) Z ( β ) (cid:12)(cid:12)(cid:12)(cid:12) dt = 1 Z ( β ) (cid:88) E n E e − β E (4.7)– 21 – - - - - - - - t g ( t ) - - - - - - - t g c ( t ) - - - - - - - - - t g d ( t ) - - - - - - - t g ( t ) - - - - - - - t g c ( t ) - - - - - - - - - t g d ( t ) - - - - - - - t g ( t ) - - - - - - - t g c ( t ) - - - - - - - - - t g d ( t ) N = 10 N = 12 N = 14 N = 16 N = 18 N = 20 N = 22 N = 24 N = 26 N = Figure 6 . The “spectral form factors” g ( t ), g c ( t ) and g d ( t ) in the supersymmetric SYK model,treating the supercharge matrix as the Hamiltonian, with J N =1 = 1, β = 0 , , 10 respectively. for large enough T ( n E means the degeneracy). For a non-degenerated spectrum, it shouldhave a simple formula (cid:12)(cid:12)(cid:12)(cid:12) Z ( β, t ) Z ( β ) (cid:12)(cid:12)(cid:12)(cid:12) = Z (2 β ) Z ( β ) (4.8)However, for a continuous spectrum, the quantity has vanishing long-term average. Thus, thequantity should be an important criterion to detect the discreteness. In this paper, we will– 22 –se a similar quantity, which is called the spectral form factor g ( t, β ) = (cid:104) Z ( β + it ) Z ( β − it ) (cid:105)(cid:104) Z ( β ) (cid:105) g d ( t, β ) = (cid:104) Z ( β + it ) (cid:105) (cid:104) Z ( β − it ) (cid:105)(cid:104) Z ( β ) (cid:105) g c ( t, β ) = g ( t, β ) − g d ( t, β ) = (cid:104) Z ( β + it ) Z ( β − it ) (cid:105) − (cid:104) Z ( β + it ) (cid:105) (cid:104) Z ( β − it ) (cid:105)(cid:104) Z ( β ) (cid:105) (4.9)In the SYK model, these quantities will have similar predictions with the hamiltonian replacedby random matrix from some specific given Dyson ensembles. For example, for a givenrealization M from a random matrix ensemble with large L , we have the analytically continuedpartition function Z rmt ( β, t ) = 1 Z rmt (cid:90) dM ij exp (cid:18) − L M ) (cid:19) Tr( e − βM − iMt ) (4.10)where Z rmt = (cid:90) dM ij exp (cid:18) − L M ) (cid:19) (4.11)The properties of spectral form factors given by random matrix theory, g rmt ( t ), have beenstudied in [54]. There are three specific periods in g rmt ( t ). In the first period, the spectralform factor will quickly decay to a minimal until dip time t d . Then after a short increasing(the ramp ) towards a plataeu time t p , g rmt ( t ) will arrive at a constant plataeu. This patternis extremely similar with SYK model. Theoretically, in the early time (before t d ), g ( t ) shouldnot obtained by g rmt ( t ) because of different initial dependence on energy, while in the latetime these two systems are conjectured to be coincide [54].With the data of energy eigenvalues one could compute the spectral form factors, whichhave been shown in Figure 5 for supersymmetric SYK model. We perform the calculation forthree different functions g ( t ), g d ( t ) and g c ( t ) with β = 0 , , 10 and several N s. Clear patternssimilar with random matrix theory predictions are shown in these numerical simulations. Onecould directly see the dip, ramp and plateau periods. For small β s there exist some small vi-brations in the early time, while for large β this effect disappears. The function g d is stronglyvibrating because we have only finite number of samples. One could believe that the infinitenumber of samples will cancel the noisy randomness of the curves.A clear eight-fold correspondence has been shown in the spectral form factor. Near theplateau time of g ( t ) one should expect roughly a smooth corner for GOE-type, a kink forGUE-type, and a sharp peak for GSE-type. Thus, we observe roughly the smooth cornersfor N = 14 , , , 24, while the sharp peaks for N = 18 , , , 28 (although the peaks looknot very clear because of finite sample size). For N = 10 , 12, as shown in Figure 3 there is– 23 – ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ 10 12 14 16 18 20 22 24 26 2810 10 12 14 16 18 20 22 24 26 28N d i p t i m e t d ● β = ■ β = ◆ β = ● ● ● ● ● ● ● ● ● ● 10 12 14 16 18 20 22 24 26 2810050010005000 10 12 14 16 18 20 22 24 26 28N D i p t i m e t d J ● β = t d ⅇ + ■ ■ ■ ■ 10 12 14 16 18 20 22 24 26 2810050010005000 10 12 14 16 18 20 22 24 26 28N D i p t i m e t d J β = t d , N = ( mod8 ) ⅇ + ■ β = t d , N = ( mod8 ) ⅇ + Figure 7 . The dip time t d for supersymmetric SYK model. In the left figure, we evaluate threedifferent temperatures and compute the dip time with respect to N , where the error bar is given asthe standard deviation when evaluating t d because of large noise is around the minimal point of g ( t ).In the middle figure we fit the dip time by polynomials and exponential functions for t d ( N ) at thetemperature β = 5. In the right figure we separately fit the dip time for two different random matrixclasses with the same temperature β = 5 and the same fitting functions. no clear random matrix correspondence because N is too small, thus we only observe somevibrations near the plateau time.We also perform a similar test on the supercharge Q , plotted in Figure 6. In Section 4.2,we numerically tested the nearest neighbour level statistics of Q which matches perfectly thestatistics of the corresponding H . The spectral form factors of Q are slightly different fromthose of H , yet they show exactly the same eight fold behavior. More quantitative data could be read off from the spectral form factors. In Figure 7, Figure8 and Figure 9 we present our numerical results for dip time t d of g ( t ), plateau time t p of g ( t ), and plateau height g d of g c ( t ) respectively. For numerical technics, we choose the linearfitting in the ramp period, and the plateau is fitted by a straight line parallel to the time axis.The dip time is read off as the averaged minimal point at the end of the dip period, and theerror bar could be computed as the standard deviation.It is claimed in [54] that polynomial and exponential fitting could be used to interpret thedip time as a function of N with fixed temperature. We apply the same method to thesupersymmetric extension. However, we find that in the supersymmetric extension, the fit-ting is much better if we fit the GOE-type group ( N mod 8 = 0 , 6) and the GSE-type group( N mod 8 = 2 , 4) separately. On the other hand, although we cannot rule out the polyno-mial fitting, the fitting effect of exponential function is relatively better. On the exponentialfittings with respect to different degeneracy groups, the coefficients before N are roughly thesame (0 . N for β = 5) while the overall constants are different. That indicates that theeight-fold degeneracy class or random matrix class might influence the overall factors in thedip time exponential expressions.One could also read off the plateau time and exponentially fit the data. Similar with dip– 24 – ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ 10 12 14 16 18 20 22 24 26 2810 N P l a t eau t i m e t b J β = t p ,allN ⅇ + ■ β = t p ,allN ⅇ + ◆ β = t p ,allN ⅇ + ■ ■ ■ ■ ■ ■□ □ □ □▲ ▲ ▲ ▲ ▲ ▲△ △ △ △● ● ● ● ● ●○ ○ ○ ○ 10 12 14 16 18 20 22 24 26 2810 N P l a t eau t i m e t b J ■ β = t p ,N = ( mod8 ) ⅇ + □ β = t p ,N = ( mod8 ) ⅇ + ▲ β = t p ,N = ( mod8 ) ⅇ + △ β = t p ,N = ( mod8 ) ⅇ + ● β = t p ,N = ( mod8 ) ⅇ + ○ β = t p ,N = ( mod8 ) ⅇ + Figure 8 . The plateau time t p for supersymmetric SYK model. We choose three different temper-atures and evaluate the plateau time with respect to N , and we use the exponential function to fit t p ( N ). In the left figure we use all N s, while in the right figure we separately fit two different randommatrix classes. time, we could also separately fit the plateau time with respect to two different random ma-trix classes, and one could find a difference in the overall factors of these two groups, whilethe coefficients before N are the same. There is a non-trivial check here. Theoretically fromrandom matrix theory one can predict that the plateau time scales like t p ∼ e S (2 β ) [54]. Inthe large β limit, the entropy should be roughly the ground state entropy. Analytically, theentropy is predicted by S ( β = ∞ ) = N s = 0 . N . Now check the largest β we take( β = 10), we can read off the entropy by 0 . N (GSE-type), 0 . N (GOE-type), or 0 . N (two groups together), which perfectly matches our expectation. ● ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○□ □ □ □ □ □ □ □ □ □ □ □◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ △ △ △ △ △ △ △ △ △ △ △ △▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ ▽ - β p l a t eauhe i gh t ● N = ■ N = ◆ N = ▲ N = ▼ N = ○ N = □ N = ◇ N = △ N = ▽ N = ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ 10 15 20 251. × - × - p l a t eauhe i gh t ● β = ■ β = ◆ β = ● ( β )/ Z ( β ) ^2, β = ■ ( β )/ Z ( β ) ^2, β = ▲ ( β )/ Z ( β ) ^2, β = ● ( β )/ Z ( β ) ^2, β = ■ ( β )/ Z ( β ) ^2, β = ▲ ( β )/ Z ( β ) ^2, β = Figure 9 . The plateau height g p for supersymmetric SYK model. In the left figure we choose severaltemperatures and fix N in each curve, while in the right we fix β and evaluate g p ( N ). For the plateau height, one can clearly see an eight-fold structure. From the previous dis-cussion we obtain that the plateau height should equals to Z (2 β ) /Z ( β ) times a contributionfrom the degeneracy, which is clearly shown in the figure. For N = 14 , , , 24 (GOE-type),the degeneracy is two thus points should be on the lower line, while for N = 18 , , , Conclusion and outlook In this paper, we use analytic arguments and numerical evidence to explore the supersym-metric constraints on the random matrix theory symmetry class. We focus on the N = 1supersymmetric SYK model, a supersymmetric generalization of nonlocal-coupled majonarafermions with similar chaotic behavior for a two dimensional quantum black hole.Use the direct classification from random matrix theory, we show that for N = 1 super-symmetric SYK model has a different behavior for N mod 8 structure. These argumentsmight be made to be more general: supersymmetry could directly change the universal classof Hamiltonian (GOE/GUE/GSE) by classifying the symmetry class of supercharge, wherecombinations of Witten index and antiunitary operators will make some new anti-unitaries;On the other hand, the quadratic structure of the Hamiltonian will change the original typeof distribution from Gaussian to Wishart-Laguerre. These points may happen for genericsupersymmetric statistical physics models.We also use numerical method, exact diagonalization to confirm the random matrix theoryclassification on the Hamiltonian and the supercharge of the supersymmetric SYK model.It is clear that if we check the spectrum density, the supercharge Q shows a clear featurefrom one-point function of extended random matrix theory ensembles, while the Hamiltonianshows a feature of quadratic semi-circle (Marchenko-Pastur). However, for level statistics(eg. Wigner surmise and spectral form factor), the universal class GSE/GOE could captureimportant physical features, and the new eight-fold rule could be verified.Several future directions could be investigated. Firstly, one may consider higher supersym-metry constraints on the SYK model, such as N = 2 generalization. Many thermodynamicaland field theory properties of higher supersymmetric SYK theory are non-trivial, and it mightbe interesting to connect these properties to random matrix theory. Moreover, to understandthe spectral form factor with supersymmetric constraints, one could also try to study super-conformal field theory partition functions at late time. Finally, introducing supersymmetriesin the symmetry classification of phases in the condensed matter theory will bring more un-derstanding at the boundary of condensed matter and high energy physics. We leave theseinteresting possibilities to future works. Acknowledgments We thank Xie Chen, Kevin Costello, Liam Fitzpatrick, Davide Gaiotto, Yingfei Gu, NicholasHunter-Jones, Alexei Kitaev, Andreas Ludwig, Evgeny Mozgunov, Alexandre Streicher forvaluable discussions. We thank Takuya Kanazawa for comments on the draft. JL is deeplygrateful to Guy Gur-Ari for communications on the symmetry of the original and supersym-– 26 –etric SYK models. TL, JL, YX and YZ are supported by graduate student programs ofUniversity of Nebraska, Caltech, Boston University and Perimeter Institute. A Review on Altland-Zirnbauer theory In this appendix we make a brief review the Altland-Zirnbauer theory (eg., see [2, 3]) thatbrings hamiltonians to ten different random matrix classes. In a physical system, symmtriescan appear and they consist a group G , then the space of physical states is a projective rep-resentation of the symmetry group. A fundamental question we can ask is, what is the mostgeneral type of hamiltonian the system can have.We may visit the simplest example to get some intuitions. The action of an element of G on the Hilbert space V can be either unitary or antiunitary, thus there is a homomorphismfrom group G to Z which labels unitarity of operators. Let G be the subgroup of unitaryoperators, then V splitts into irreps of G : V = (cid:77) i V i ⊗ C m i (A.1)where V i are irreps and m i are their multiplicities in V . If there is no antiunitary operatorsthen followed by Schur’s lemma, the most general Hamitonians are those belong to the set (cid:77) i End G ( V i ⊗ C m i ) = (cid:77) i End( C m i ) (A.2)plus Hermicity. This is called Type A in the Altland-Zirnbauer theory, without any anti-unitary operators. The case with the presence of antiunitary operators is more complicated.Let T be an antiunitary operator, then the conjugation by T , i.e. U (cid:55)→ T U T − , is anautomorphism of G , thus T maps a component V i ⊗ C m i to another V j ⊗ C m j . A simple caseis when i (cid:54) = j , which is easy to see that the most general hamiltonian is of form [2, 3]( H, T HT − ) (A.3)where H is an Hermitian operator in component i and T HT − acts on component j . Thusit’s also of Type A.The Type A is the simplest structure without any further symmetries. However, if we con-sider i = j , and consider more anti-unitary operators, the situation is much more technical. Itturns out that possible hamiltonians with specific symmetric structures can be classified intoten classes. Here we skip the detailed analysis and directly present the final results. Theseclasses are classified by the following three different operators, • T + , antiunitary, commutes with hamiltonian, and T = ± • T − , antiunitary, anticommutes with hamiltonian, and T − = ± 1– 27 – Λ, unitary, anticommutes with hamiltonian, and Λ = 1If two of these three operators exist, the third will be determined by the following identity,Λ = T + T − (A.4)The properties of these three operators can classify the hamiltonian in the following tenclasses, T T − Λ Cartan label Block TypeA (GUE) M complex: M † = M C M real: M T = M R − M quaternion: M † = M H (cid:32) ZZ † (cid:33) Z complex C − (cid:32) A B ¯ B − ¯ A (cid:33) A Hermitian B complex symmetric C M pure imaginary, skew-symmetric C (cid:32) AA T (cid:33) A real R − (cid:32) Z ¯ Z (cid:33) Z complex symmetric R − (cid:32) Y − ¯ Y (cid:33) Y complex, skew-symmetric H − − (cid:32) BB † (cid:33) B quaternion H where there are no values in some corresponding operators we mean that there is no such asymmetry in the system. We also present the block representation in this table, where blocksare classified by the ± B Eigenvalue distribution This appendix is a simple introduction on the random matrix theory eigenvalue distribution(for instance, see [73, 74]), the measure in the eigenvalue basis. For Wigner-Dyson ensemble,this is given by the formula P ( λ ) dλ = C ( N, ˜ β ) | ∆( λ ) | ˜ β (cid:89) k e − N ˜ β λ k dλ k (B.1)– 28 –here λ = ( λ , · · · , λ N ) is the set of eigenvalues, ∆( λ ) is the Vandermont determinant definedby ∆( λ ) = (cid:89) k>l ( λ k − λ l ) (B.2)and C ( N, ˜ β ) is a normalization constant depending on ˜ β and N . For different ensembles, ˜ β is defined as RMT ˜ β AI(GOE) 1A(GUE) 2AII(GSE) 4For the remaining ensembles, the eigenvalues occur in pairs (because the T − operator intro-duced in the last appendix anticommutes with Q ), and the eigenvalues probability distributionis given by P ( λ ) dλ = C ( N, ˜ β, ˜ α ) | ∆( λ ) | ˜ β (cid:89) k λ ˜ αk e − N ˜ β λ k dλ k (B.3)where we only take the positive one from a pair of eigenvalues, and C ( N, ˜ β, ˜ α ) is definedalso as the corresponding normalization constant. In the Altland-Zirnbauer classification,constants ˜ α and ˜ β are set as (considering the real model of us, we have set the flavor number N f = 0 and the topological index ν = 0 in chiral ensembles)RMT ˜ β ˜ α BDI(chGOE) 1 0AIII(chGUE) 2 1CII(chGSE) 4 3CI(BdG) 1 1D(BdG) 2 0C(BdG) 2 2DIII(BdG) 4 1We will also need the eigenvalue distribution of the hamiltonian which is the square of Q ,so we can take the square distribution of B.3, which will change Gaussian distribution toWishart-Laguerre, which is P ( λ ) dλ = C (cid:48) ( N, ˜ β, ˜ α ) | ∆( λ ) | ˜ β (cid:89) k λ ˜ α − k e − N ˜ β λ k dλ k (B.4)here λ k are nonnegative and C (cid:48) ( N, ˜ β, ˜ α ) is a new normalization constant which is one half of C ( N, ˜ β, ˜ α ). We could also write P ( λ ) dλ ∼ | ∆( λ ) | ˜ β (cid:89) k λ ˜ µk e − N ˜ β λ k dλ k (B.5)– 29 –here ˜ µ = ( ˜ α − / 2. The following table summarize the related index for supersymmetricSYK model N mod 8 Q ˜ α ˜ β ˜ µ − / 22 DIII (BdG) 1 4 04 CII (chGSE) 3 4 16 CI (BdG) 1 1 0In N mod 8 = 0 , 4, the index ˜ µ precisely matches Wishart matrix. Moreover, Although theresult has ˜ µ dependence for N mod 8 = 2 , 6, which does not precisely match Wishart matrixfrom Dyson Gaussian ensemble by index ˜ µ , we could also use the terminology LOE/LSE torefer the universal class from squaring of Gaussian matrix, similar with Altland-Zirnbauerclassification as a subset of Dyson, regardless multiple anti-unitary symmetries. Thus, wecall N mod 8 = 0 , , , References [1] F. J. Dyson, J. Math. Phys. , 140 (1962).[2] M. R. Zirnbauer, J. Math. Phys. , 4986 (1996) [math-ph/9808012].[3] A. Altland and M. R. Zirnbauer, Phys. Rev. B , 1142 (1997).[4] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New Journal of Physics 12, 065010(2010), arXiv:0912.2157 [cond-mat.mes-hall].[5] L. Fidkowski and A. Kitaev, Phys. Rev. B , 075103 (2011), arXiv: 1008.4138[cond-mat.str-el].[6] A. Kitaev, Talks at KITP, April 7, 2015, May 27, 2015, and Feb. 12, 2015; A. Kitaev, Talkgiven at the Fundamental Physics Prize Symposium, Nov. 10, 2014.[7] S. Sachdev and J. Ye, Phys. Rev. Lett. , 3339 (1993) [cond-mat/9212030].[8] A. Almheiri and J. Polchinski, JHEP , 014 (2015) [arXiv:1402.6334 [hep-th]].[9] W. Fu and S. Sachdev, Phys. Rev. B , no. 3, 035135 (2016) [arXiv:1603.05246[cond-mat.str-el]].[10] J. Polchinski and V. Rosenhaus, JHEP , 001 (2016) [arXiv:1601.06768 [hep-th]].[11] A. Jevicki, K. Suzuki and J. Yoon, JHEP , 007 (2016) [arXiv:1603.06246 [hep-th]].[12] J. Maldacena and D. Stanford, Phys. Rev. D , no. 10, 106002 (2016) [arXiv:1604.07818[hep-th]].[13] K. Jensen, Phys. Rev. Lett. , no. 11, 111601 (2016) [arXiv:1605.06098 [hep-th]].[14] M. Cvetic and I. Papadimitriou, JHEP , 008 (2016) Erratum: [JHEP , 120 (2017)][arXiv:1608.07018 [hep-th]].[15] A. Jevicki and K. Suzuki, JHEP , 046 (2016) [arXiv:1608.07567 [hep-th]]. – 30 – 16] D. Bagrets, A. Altland and A. Kamenev, Nucl. Phys. B , 191 (2016) [arXiv:1607.00694[cond-mat.str-el]].[17] J. Maldacena, D. Stanford and Z. Yang, PTEP , no. 12, 12C104 (2016) [arXiv:1606.01857[hep-th]].[18] Y. Gu, X. L. Qi and D. Stanford, arXiv:1609.07832 [hep-th].[19] D. J. Gross and V. Rosenhaus, arXiv:1610.01569 [hep-th].[20] M. Berkooz, P. Narayan, M. Rozali and J. Simon, JHEP , 138 (2017) [arXiv:1610.02422[hep-th]].[21] W. Fu, D. Gaiotto, J. Maldacena and S. Sachdev, Phys. Rev. D , no. 2, 026009 (2017)[arXiv:1610.08917 [hep-th]].[22] O. Parcollet and A. Georges, Phys. Rev. B, , 5341 (1999) [cond-mat/9806119].[23] A. Georges, O. Parcollet, and S. Sachdev, Phys. Rev. Lett. , 840 (2000), cond-mat/9909239.[24] P. Hayden and J. Preskill, JHEP , 120 (2007) [arXiv:0708.4025 [hep-th]].[25] S. Sachdev, Phys. Rev. Lett. , 151602 (2010) [arXiv:1006.3794 [hep-th]].[26] D. Anninos, T. Anous, P. de Lange and G. Konstantinidis, JHEP , 066 (2015)[arXiv:1310.7929 [hep-th]].[27] S. Sachdev, Phys. Rev. X (2015) no.4, 041025 [arXiv:1506.05111 [hep-th]].[28] E. Perlmutter, JHEP , 069 (2016) [arXiv:1602.08272 [hep-th]].[29] D. Anninos, T. Anous and F. Denef, JHEP , 071 (2016) [arXiv:1603.00453 [hep-th]].[30] I. Danshita, M. Hanada and M. Tezuka, arXiv:1606.02454 [cond-mat.quant-gas].[31] P. Betzios, N. Gaddam and O. Papadoulaki, JHEP , 131 (2016) [arXiv:1607.07885[hep-th]].[32] L. Garca-Alvarez, I. L. Egusquiza, L. Lamata, A. del Campo, J. Sonner and E. Solano,arXiv:1607.08560 [quant-ph].[33] D. A. Roberts and B. Yoshida, JHEP , 121 (2017) [arXiv:1610.04903 [quant-ph]].[34] E. Witten, arXiv:1610.09758 [hep-th].[35] A. A. Patel and S. Sachdev, arXiv:1611.00003 [cond-mat.str-el].[36] I. R. Klebanov and G. Tarnopolsky, arXiv:1611.08915 [hep-th].[37] M. Blake and A. Donos, arXiv:1611.09380 [hep-th].[38] T. Nishinaka and S. Terashima, arXiv:1611.10290 [hep-th].[39] R. A. Davison, W. Fu, A. Georges, Y. Gu, K. Jensen and S. Sachdev, arXiv:1612.00849[cond-mat.str-el].[40] N. Sannomiya, H. Katsura and Y. Nakayama, arXiv:1612.02285 [cond-mat.str-el].[41] D. Anninos and G. A. Silva, arXiv:1612.03795 [hep-th].[42] C. Peng, M. Spradlin and A. Volovich, arXiv:1612.03851 [hep-th].[43] Y. Liu, M. A. Nowak and I. Zahed, arXiv:1612.05233 [hep-th]. – 31 – 44] C. Krishnan, S. Sanyal and P. N. B. Subramanian, arXiv:1612.06330 [hep-th].[45] J. M. Magan, arXiv:1612.06765 [hep-th].[46] G. Turiaci and H. Verlinde, arXiv:1701.00528 [hep-th].[47] F. Ferrari, arXiv:1701.01171 [hep-th].[48] A. M. Garca-Garca and J. J. M. Verbaarschot, arXiv:1701.06593 [hep-th].[49] Z. Bi, C. M. Jian, Y. Z. You, K. A. Pawlak and C. Xu, arXiv:1701.07081 [cond-mat.str-el].[50] W. W. Ho and D. Radicevic, arXiv:1701.08777 [quant-ph].[51] Y. Z. You, A. W. W. Ludwig and C. Xu, arXiv:1602.06964 [cond-mat.str-el].[52] A. M. Garcia-Garcia and J. J. M. Verbaarschot, Phys. Rev. D , no. 12, 126010 (2016)[arXiv:1610.03816 [hep-th]].[53] E. Dyer and G. Gur-Ari, arXiv:1611.04592 [hep-th].[54] J. S. Cotler et al. , arXiv:1611.04650 [hep-th].[55] M. F. Sohnius, Phys. Rept. , 39 (1985).[56] K. Efetov, “Supersymmetry in disorder and chaos.”[57] S. S. Lee, arXiv:1009.5127 [hep-th].[58] R J Muirhead, “Aspects of multivariate statistical theory.”[59] E. Witten, Nucl. Phys. B , 513 (1981).[60] F. Cooper and B. Freedman, Annals Phys. , 262 (1983).[61] F. Cooper, A. Khare and U. Sukhatme, Phys. Rept. , 267 (1995) [hep-th/9405029].[62] V. K. Oikonomou and K. Kleidis, Int. J. Theor. Phys. , no. 3, 933 (2015) [arXiv:1408.1552[hep-th]].[63] V. K. Oikonomou, Int. J. Geom. Meth. Mod. Phys. , 1550114 (2015)doi:10.1142/S0219887815501145 [arXiv:1412.0244 [hep-th]].[64] J. Mateos Guilarte and A. Moreno Mosquera, Eur. Phys. J. Plus , no. 2, 93 (2017)doi:10.1140/epjp/i2017-11362-7 [arXiv:1608.06211 [math-ph]].[65] L. Brink, T. H. Hansson, S. Konstein and M. A. Vasiliev, Nucl. Phys. B , 591 (1993)[hep-th/9302023].[66] X. Bekaert, S. Cnockaert, C. Iazeolla and M. A. Vasiliev, hep-th/0503128.[67] M. S. Plyushchay, Mod. Phys. Lett. A , 397 (1996) [hep-th/9601141].[68] M. S. Plyushchay, Annals Phys. , 339 (1996) [hep-th/9601116].[69] L. Van Hove, Nucl. Phys. B , 15 (1982).[70] Junker, Georg. “Supersymmetric methods in quantum and statistical physics”, p. 12-14,Springer Science & Business Media, 2012.[71] F. Calogero, J. Math. Phys. , 2191 (1969).[72] L. Brink, T. H. Hansson and M. A. Vasiliev, Phys. Lett. B , 109 (1992) [hep-th/9206049]. – 32 – 73] Dumitriu, Ioana, and Alan Edelman. ”Matrix models for beta ensembles.” Journal ofMathematical Physics 43.11 (2002): 5830-5847.[74] Edelman, Alan, and N. Raj Rao. ”Random matrix theory.” Acta Numerica 14 (2005): 233-297.[75] Nica A, Speicher R. Lectures on the combinatorics of free probability[M]. Cambridge UniversityPress, 2006.[76] A. D. Jackson, M. K. Sener and J. J. M. Verbaarschot, Nucl. Phys. B , 707 (1996)[hep-ph/9602225].[77] Atas, Y. Y. and Bogomolny, E. and Giraud, O. and Roux, G., Phys. Rev. Lett. , 084101(2013), 084101(2013)