Supersymmetric SYK model and random matrix theory

In this paper, we investigate the effect of supersymmetry on the symmetry classification of random matrix theory ensembles. We mainly consider the random matrix behaviors in the $\mathcal{N}=1$ supersymmetric generalization of Sachdev-Ye-Kitaev (SYK) model, a toy model for two-dimensional quantum black hole with supersymmetric constraint. Some analytical arguments and numerical results are given to show that the statistics of the supersymmetric SYK model could be interpreted as random matrix theory ensembles, with a different eight-fold classification from the original SYK model and some new features. The time-dependent evolution of the spectral form factor is also investigated, where predictions from random matrix theory are governing the late time behavior of the chaotic hamiltonian with supersymmetry.


Introduction
Physical systems with some stochastic or chaotic properties have some randomness in the setup of the fundamental hamiltonian, which could be effectively simulated in the context of random matrix theory. When choosing an ensemble from random matrix theory for a chaotic hamiltonian, we often need to consider the symmetries in the dynamics of the related physical system. The choice of standard matrix ensembles from symmetries, historically comes from the invention of Dyson [1], which is called three-fold way when classifying Gaussian Unitary Ensemble (GUE), the Gaussian Orthogonal Ensemble (GOE), and the Gaussian Symplectic Ensemble (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 the practical usage, one of the most celebrated works would be the classification of interaction inside topological insulators and topological phases in a ten-fold way [4,5].
In the recent study, the rising interests of studies on Sachdev-Ye-Kitaev (SYK) model gives another profound application in the random matrix theory classification. SYK model [6,7] is a microscopic quantum hamiltonian with random Gaussian non-local couplings among majonara fermions. As is maximally chaotic and nearly conformal, this model could be treated as a holographic dual of quantum black hole with AdS 2 horizon through the (near) AdS/CFT correspondence [8][9][10][11][12][13][14][15][16][17]. In the recent research people have also discussed several generalizations of the SYK model [18][19][20][21], such as higher dimensional generalizations and supersymmetric constraints. Some other related issues and similar models are discussed in . In the recent discussions, people have discovered that the SYK hamiltonian has a clear correspondence with the categories of the three fold standard Dyson ensembles, unitary, orthogonal and sympletic ensembles, in the random matrix theory [51][52][53][54]. In the recent work, [53,54], it is understood that the time-dependent quantum dynamics of the temperature-dependent spectral form factor, namely, the combinations of partition functions with a special analytic continuation in SYK model, is computable in the late time by form factors in the random matrix theory with the same analytic continuation, as a probe of the discrete nature of the energy spectrum in a quantum black hole, and also a solid confirmation on the three-fold classification [54].
In the route towards Dyson's classification, one only considers the set of simple unitary or anti-unitary operators as symmetries when commuting or anticomuting with the hamiltonian. An interesting question would be, what is the influence of supersymmetry, the symmetry between 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 influences in the study of disorder system and statistical physics [56], and could be emergent from condensed matter theory models [57]. Originating from particle physics, supersymmetry will enlarge the global symmetry group in the theory, has fruitful algebras and strong mathematical power used in several models in quantum mechanics and quantum field theory, and is extremely useful to simplify and clarify classical or quantized theories. In the recent study of SYK model, the supersymmetric generalization for the original SYK has been discussed in 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 in two dimension in a supersymmetric theory, and also a related conjecture in [54] for spectral form factor and correlation functions in super Yang-Mills theory.
In order to explore the supersymmetric constraints on the random matrix theory classification, in this paper we will study the symmetry classification and random matrix behavior of the N = 1 supersymmetric extension of SYK model by Fu-Gaiotto-Maldacena-Sachdev's prescription [21]. The effect of supersymmetry in the symmetry classification could be summarized in the following aspects, • Supersymmetry will cause the hamiltonian to show a quadratic expression. Namely, we could write H as the square of Q. This condition will greatly change the distribution of the eigenvalues. From random matrix language [58], if Q is a Gaussian random matrix, then H should be in a Wishart-Laguerre random matrix, with the eigenvalue distribution 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 a positivity condition for all eigenvalues. Namely, if Q has the eigenvalue distribution that eigenvalues come in pair with positive and negative signs, the squaring Q will cause 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 is not even. This will change the previous classification in the hamiltonian based on the representation of Clifford algebra from mathematical point of view.
• We find the Witten index or Witten parity operator (−1) F , which is well-known as a criterion for supersymmetry breaking [55,[59][60][61], is crucial in classifying the symmetry class for supercharge Q. Some evidence of this point also could be found in some other models or setups. For instance, Witten parity is the Klein operator which separates the bosonic and fermionic sectors in the N = 2 supersymmetric systems [62,63]. [64] provides a more nontrivial example, where the odd parity operators are used to move states along a chain of different fermion sectors. Reversely, in some systems where one can define a graded algebra, Klein operator serves as a key factor in realizing supersymmetry, which is helpful in models of bosonization and higher spin theories, etc. [65][66][67][68]. For example, [67] constructs the bosonized Witten supersymmetric quantum mechanics by realizing the Klein operator as a parity operator. [68] realize a Bose-Fermi transformation 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 that incorporating the Witten operator is crucial in some computation in supersymmetric systems with finite temperature. In the supersymmetric SYK model we are considering, Witten parity and the anti-unitary operator together become a new anti-unitary operator, which will significantly enlarge the set of symmetries in the hamiltonian, and change 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 construction and thermodynamics of SYK model and its supersymmetric extensions. In Section 3 we will discuss the random matrix classification for models, especially supersymmetric extensions of the SYK model. In Section 4 we will present our numerical confirmation for symmetry classifications from the exact diagonalization, including the computation of the density of states and spectral form factors. In Section 5, we will arrive at a conclusion and discuss the directions for future works. In the appendix, we will review some knowledge to make this paper self-contained, including basics on Altland-Zirnbauer theory and a calculation on the random matrix theory measure.

Introduction on models
In this paper, we will mostly focus on SYK models and their extensions. Thus before the main investigation, we will provide a simple introduction on the necessary knowledge of related models to be self-contained.

The SYK model
In this part, we will simply review the SYK model mainly following [12]. The SYK model is a microscopic model with some properties of quantum black hole. The hamiltonian 1 is given by where ψ i are Majorana fermions and they are coupled by the four point random coupling with Gaussian distribution where J SYK and J SYK are positive constants, and J SYK = √ 2J SYK . The large N partition function is given by where E 0 is the total ground state energy proportional to N and it is roughly E 0 = −0.04N [54]. s 0 is the ground state entropy contributed from one fermion, and one can estimate it theoretically [12], where G is the Catalan number. c is the specific heat, which could be computed by and α S = 0.0071 is a positive constant. This contribution c/β is from the Schwarzian, the quantum fluctuation near the saddle point of the effective action in the SYK model. The Schwarzian partition function is where the path integral is taken for all possible reparametrizations τ (u) of the thermal circle in different equivalent classes of the SL(2, R) symmetry. The Schwarzian corresponds to the broken reparametrization symmetry of the SYK model. One can compute the one-loop correction from the soft mode of the broken symmetry, As a result, one can consider the correction from the soft mode if we consider an external one-loop factor (βJ SYK ) −3/2 . The density of states could be also predicted by the contour integral of the partition function as Following [21], in the supersymmetric extension of SYK model, firstly we define the super- for Majonara fermions ψ i . C ijk is a random tensor with the Gaussian distribution as the coupling, where J N =1 is also a constant with mass dimension one. The square of the supercharge will give the hamiltonian of the model 1 One could also generalize the SYK model to general q point non-local interactions where q are even numbers larger than four. The hamiltonian should be where Sometimes we will discuss the general q in this paper but we will mainly focus on the q = 4 case.
where [· · · ] is the summation of all possible antisymmetric permutations. Besides the shifted constant E c , the distribution of J ijkl is different from the original SYK model because it is not a free variable of Gaussian distribution, which changes the large N behavior of this model. In the large N limit, the model has an unbroken supersymmetry with a bosonic superpartner b i . The Lagrangian of this model is given by In this model, the Schwarzian is different from the original SYK model. We also have the expansion for the large N partition function But the results of E 0 and s 0 are different (while the specific heat is the same for these two models). In the large N limit, the supersymmetry is preserved, thus we have the ground state energy E 0 = 0. The zero temperature entropy is given by Moreover, the one-loop correction from Schwarzian action is different. As a result of supersymmetry constraint, the one-loop factor is (βJ which predicts a different behavior for the density of states For the generic positive integerq we can also define the N = 1 supersymmetric extension with non-local interaction of 2q − 2 fermions. The supercharge should be where Andq = 3 will recover the case in the main text.

Random matrix classification
It is established that SYK model is classified by random matrix theory in that the random interacting SYK hamiltonian fall into one of the three standard Dyson ensembles in the eight-fold way [51][52][53][54]. It is natural to believe that the supersymmetric extension can also be described by random matrix theory. To sharpen the argument, we derive the exact correspondence between each SYK hamiltonian and some random matrix ensembles, in other words, the eight-fold rule for supersymmetric case. A priori, the supersymmetric SYK hamiltonian 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 two major differences, which have been also mentioned in the previous discussions.
• The degeneracy of the two hamiltonian matrices are different. The degeneracy of supersymmetric SYK model is also investigated by [21], which we derive again using some different discussion in Section 3.2.2. The degeneracy space is enlarged by supersymmetry. Generally, the energy level distribution of random matrices is sensitive to the degeneracy and is thus sensitive to the supersymmetric extension.
• Another difference is the apparent positive semidefiniteness of the hamiltonian being the square of the supercharge. We will see later that the positive constraint leads to a new eigenvalue distribution different from those of Gaussian ensembles.
Symmetry analysis is crucial in classifying the random matrix statistics of hamiltonian matrices. [51,54] argue that the particle-hole symmetry operator determines the class of random matrix theory statistics. The random matrix classification dictionary is determined by the degeneracy and the special relations required by having the symmetry. The systematic method of random matrix classification is established as the Atland-Zirnbauer theory [2,3], reviewed in appendix A. The anti-unitary operators play a central role in the classifications. The Atland-Zirnbauer also applies to extended ensembles different from the three standard Dyson ensembles, which we find useful in classifying the supersymmetric SYK theory. In Section 3.1 we derive again the eight-fold way classification of original SYK hamiltonian using Atland-Zirnbauer theory and find unambiguously the matrix representations of hamiltonian in each mod eight sectors. We notice that the matrix representation of hamiltonian takes block diagonal form with each block being a random matrix from a certain ensemble. This block diagonal form is also found by [51] in a different version.
Naively one would apply the same argument to the supersymmetric hamiltonian, since it also enjoys the particle-hole symmetry. But this is not the full picture. First, one need to take into account of hamiltonian being the square of the supercharge and is thus not totally random. In Section 3.2.1 we argue that the supercharge Q has a random matrix description which falls into one of the extended ensembles. Using the Atland-Zirnbauer theory on Q we obtain its matrix representation in block diagonal form and use it to determine the matrix representation of the hamiltonian in Section 3.2.2. Second, in order to obtain the correct classification one needs to consider the full set of symmetry operators. Apparently particlehole is not enough since supersymmetry enlarges the SYK degeneracy space. We argue that the Witten index operator, (−1) F , is crucial in the symmetry analysis of any system with supersymmetry. Incorporating (−1) F we obtain the full set of symmetry operators. Finally, the squaring operation, will change the properties of the random matrix theory distribution of supercharge Q, from Gaussian to Wishart-Laguerre. The quantum mechanics and statistics in supersymmetric SYK models, based on the main investigation in this paper, might be a non-trivial and compelling example of supersymmetric symmetry class.

SYK
Now we apply the Altland-Zirnbauer classification theory (see appendix A for some necessary knowledges) to the original SYK model [51][52][53][54]. This is accomplished by finding the symmetry of the theory (and has been already discussed in other works, see [51,54]). First, one can change the majonara fermion operators to creation annihilation operators c α andc α by Hilbert space with two different charge parities. One can define the particle-hole operator where K is the complex conjugate operator (c α andc α are real). The operation of P on fermionic operators is given by From these commutation relation we can show that To compare with the Altland-Zirnbauer classification, we need to know the square of P and this is done by direct calculation Now we discover that P can be treated as a T + operator and it completely determines the class of the hamiltonian. Before we list the result, it should be mentioned that the degeneracy of 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 as a 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 2 = −1, so there is a two-fold degeneracy. There is no further independent symmetries. From Altland-Zirnbauer theory we know that in each parity block there is a GSE matrix. Also, where two copies of GSEs are independent.
• N mod 8 = 0: The symmetry P is a parity-invariant mapping and P 2 = 1. There is no further symmetries so the degeneracy is one. From Altland-Zirnbauer theory we know that in each parity 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, where the level statistics means some typical numerical evidence of random matrix, for instance, Wigner surmise, number variance, or ∆ 3 statistics, etc. Although the SYK hamiltonian can be decomposed as two different parity sectors, we can treat them as standard Dyson random 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 independent sectors, where two mixed sectors will show a many-body localized phase statistics instead of a chaotic phase statistics, which has been discussed originally in [51].) In the following we will also numerically test the random matrix behavior, and based on the numerical testing range of N we can summarize the following table for practical usage.

N = 1 supersymmetric classification
Supersymmetry algebra is a Z 2 -graded algebra, where states and operators are subdivided into 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 operators with even parity. The Klein operator of supersymmetry algebra is naturally the Witten index operator.
Witten index might plays a role in the symmetric structure and block decomposition in the supersymmetric quantum mechanics. A simple example is [70], in N = 2 supersymmetry algebra, Define W be the Witten operator. The Witten operator has eigenvalue ±1 and separates the Hilbert space into two parity sectors We can also define projection operators P ± = 1/2(1 ± W ). In the parity representation the operators take 2 × 2 block diagonal form Because of Q 2 = 0 and {Q, W } = 0 the complex supercharges are necessarily of the form which imply In the above equation, A takes H − → H + and its adjoint A † takes H + → H − . The supersymmetric hamiltonian becomes diagonal in this representation In this construction, the Hilbert space is divided by Witten parity operator. The hamiltonian is shown to take the block diagonal positive semidefinite form without even referring to the explicit construction of the hamiltonian. It is remarkable that the above computation is very similar to our work from Section 3.2.1 to 3.2.2. Applications of this property can be found in [62,63]. They describe a supersymmetric Quantum Mechanics system where fermions scatter 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 2 sector is simply the kernel of the differential operator and the Witten index is computed. A more non-trivial example is provided by [64]. In this work, the Hilbert space is divided into an N fermions Fock space. Thus the Hamiltonian 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 a bosonic system such as a Calogero-like model [71], a system of one dimensional Harmonic oscillators with inverse square interactions and extensions. What makes supersymmetry manifest is the Klein operator. The model and its various extensions are studied in [65][66][67][68]72]. A trivial simple Harmonic operator has algebra [a − , a + ] = 1. The algebra describes a bosonic system. Z 2 grading is realized by introducing an operator K = cos(πa + a − ). The new operator anti-commutes with a − and a + thus is a Klein operator. Based on the Klein operator one can construct the projection operators on both sectors and also the supercharge. In this way the simple harmonic oscillator is "promoted" to have supersymmetry. A generalization to simple hamornic oscillator is the deformed Heisenberg algebra, [a − , a + ] = 1 + νK. The corresponding system is an N = 2 supersymmetric extension of the 2-body Calogero. The model is also used in considerably simplifying Calogero model.
These evidences strongly support the argument that supersymmetry will change the classification of symmetry class in quantum mechanical models. In the following work, we will show that supersymmetric SYK model symmetry class can be explicitly constructed and change the classification of random matrix theory ensembles.

Supercharge in N = 1 SYK
In the N = 1 supersymmetric model, it should be more convenient to consider the spectrum of Q instead of H, because H is the square of Q. Although Q is not a hamiltonian, since we only 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 about the symmtry of the theory. We notice that the Witten index (−1) F is which is the fermionic parity operator up to a sign (−1) N d . Witten index and particle-hole symmetry have the following commutation relation: Now we define a new operator, R = P (−1) F . It has a compact form R and P are both anti-unitary symmetries of Q, with commutation relations: and squares Thus, in different values of N , the two operators P and R behave different and replace the role in T + and T − in the Altland-Zirnbauer theory. Now we can list the classification for the matrix ensemble of N = 1 supersymmetric SYK model One can also write down the block representation of Q. Notice that the basis of block decomposition is based on the ±1 eigenspaces of anti-unitary operators, namely, it is decomposed based on the parity.

Hamiltonian in N = 1 theory
Now we already obtain the random matrix type of the supercharge. Thus the structure of the square of Q could be considered case by case. Before that, we can notice one general property, that unlike the GOE or GSE group in SYK, in the supersymmetric model there is a supercharge Q contains odd number of Dirac fermions as a symmetry of H, thus it always changes 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 Hermitian and H = Q 2 > 0. Thus the random matrix class of N = 1 will be some classes up to positivity constraint.
• N = 0 mod 8: In this case Q is a BDI (chGOE) matrix. Thus we can write down the block decomposition as where A is a real matrix. Thus the hamiltonian is obtained by Since AA T and A T A share the same eigenvalues ({R, Q} = 0 thus R flips the sign of eigenvalues of Q, but after squaring these two eigenvalues with opposite signatures become the same), and there is no internal structure in A (in this case P is a symmetry of Q, [P, Q] = 0, but P 2 = 1, thus P cannot provide any further degeneracy), we obtain that H has a two-fold degeneracy. Moreover, because AA T and A T A are both real positive-definite symmetric matrix without any further structure, it is nothing but the subset of GOE symmetry class with positivity condition. These two sectors will be exactly degenerated.
• N = 4 mod 8: In this case Q is a CII (chGSE) matrix. Thus we can write down the block decomposition as where B is a quaternion Hermitian matrix. Thus after squaring we obtain Since BB † and B † B share the same eigenvalues, and each block has a natural two-fold degeneracy 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 opposite signatures become the same. Also, in this case P is a symmetry of Q, [P, Q] = 0, and P 2 = −1), we get a four-fold degeneracy in the spectrum of H. Because BB † and B † B are quaternion Hermitian matrices when B is quaternion Hermitian 3 , BB † = B † B are both quaternion Hermitian positive-definite matrix without any further structure. As a 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 the block decomposition as where Y is a complex, skew-symmetric matrix. Thus after squaring we obtain Firstly let us take a look at the degeneracy. Since −YȲ and −Ȳ Y share the same eigenvalues and each block has a natural two-fold degeneracy because in skew-symmetric matrix the eigenvalues come in pair and after squaring pairs coincide (Physically it is because {P, Q} = 0 thus P flips the sign of eigenvalues of Q, but after squaring these two eigenvalues with opposite signatures become the same. Also, in this case R is a symmetry of Q, [R, Q] = 0, and R 2 = −1), we obtain a four-fold spectrum of H.
Now take the operator Q as a whole, from the previous discussion, we may note that it is quaternion Hermitian because it could be easily verified that QΩ = ΩQ and Q † = Q. Thus Q 2 = 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 parity decomposition thus in each part it is also a quarternion Hermitian matrix. Because in the total matrix it is a subset of GSE symmetry class (with positivity constraint), in each degenerated parity sector it is also in a subset of positive definite GSE symmetry class (one can see this by applying the total measure in the two different, degenerated part).
• N = 6 mod 8: In this case Q is a CI (BdG) matrix. Thus we can write down the block decomposition as where Z is a complex symmetric matrix. Thus after squaring we obtain Since ZZ andZZ share the same eigenvalues ({P, Q} = 0 thus P flips the sign of eigenvalues of Q, but after squaring these two eigenvalues with opposite signatures become the same), and there is no internal structure in Z (in this case R is a symmetry of Q, [R, Q] = 0, but R 2 = 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 as the whole matrices instead of blocks. For H we notice that the transposing operations make the exchange of these two sectors. However, the symmetric matrix statement is basis-dependent. Formally, similar with the quarternion Hermitian case, we can extend the definition of symmetric matrix by the following. Define . We can check easily that Q satisfies this condition, thus Q 2 = H must satisfy. Thus we conclude that the total matrix H in a subset of GOE symmetry class (with positivity constraint).
Although from symmetric point of view, the hamiltonian of N = 1 model should be classified in the subsets of standard Dyson ensembles. But what the subset exactly is? In fact, the special structure of the squaring from Q to H will change the distribution of the eigenvalues from Gaussian to Wishart-Laguerre [58,73,74] (Although there are some differences in the powers of terms in the eigenvalue distributions.) We will roughly called them as LOE/LUE/LSE, as has been used in the random matrix theory research. Some more details will be summarized in the appendix B.
However, the difference in the details of the distribution, beyond numerical tests of the distribution 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 be given as follows. From the supercharge point of view, because Q is in the Altland-Zirnbauer distribution with non-trivialα (see appendix B), the squaring operation will not change the level statistics such as Wigner surmise and spectral form factors (which could also be verified by 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 degeneracies have been already calculated in [21]), N mod 8 Deg. RMT Block Type Level stat.

Exact Diagonalization
In this part, we will present the main results from numerics to test the random matrix theory classification in the previous investigations. One can diagonalize the hamiltonian exactly with the representation of the Clifford algebra by the following. For operators acting on N d = N/2 qubits, one can define where σ p means standard Pauli matrices acting on the p-th qubit, tensor producting the identity matrix on the other parts, and ζ = 1, 2, ......, N d . This construction is a representation of the Clifford algebra And one can exactly diagonalize the hamiltonian by replacing the majonara fermions with gamma matrices to find the energy eigenvalues. Thus, all quantities are computable by brute force in the energy eigenstate basis.
The main results of the following investigation would be the following. In the density of supercharge eigenstates and energy eigenstates in the supersymmetric SYK model, the behavior is quite different, but coincides with our estimations from the random matrix theory classification: the spectral density of supercharge Q shows clearly the information about extended ensembles from Altland-Zirnbauer theory, and the spectral density of energy H shows a 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, GUE and GSE, the numerics from Wigner surmise and spectral form factor will show directly these eight-fold features. 3 We say a matrix M is a quaternion Hermitian matrix if and only if

Density of states
The plots for density of states in SYK model and its supersymmetric extension are shown in Figure 1 for comparison. For each realization of random hamiltonian, we compute all eigenvalues. After collecting large number of samples one can plot the histograms for all samples as the function ρ(E). For density of states in SYK model, in small N tiny vibrations are contained, while in the large N the distribution will converge to a Gaussian distribution besides the small tails. However, in the supersymmetric SYK model the energy eigenvalue structure is totally different. All energy eigenvalues are larger than zero because H = Q 2 > 0. Because of supersymmetry the lowest energy eigenvalues will approach zero for large N , and the figure will come to a convergent distribution. The shape of this distribution matches the eigenvalue 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 acquires a dip at zero, which is a clear feature for extended ensembles and could match the averaged density of eigenvalues of random matrices in CI, DIII [3] and chiral [76] ensembles at large N .

Wigner surmise
There exists a practical way to test if random matrices from a theory are from some specific ensembles. For a random realization of the hamiltonian, we have a collection of energy eigenvalues E n . If we arrange them in ascending order E n < E n+1 , we define, ∆E n = E n − E n−1 to be the level spacing, and we compute the ratio for the nearest neighbourhood spacing as r n = ∆E n /∆E n+1 . For matrices from the standard Dyson ensemble, the distribution of level spacing ratio satisfies the Wigner-Dyson statistics [77]) (which is called the Wigner surmise for GOE universal class,β = 1, Z = 8/27; for GUE universal class,β = 2, Z = 4π/(81 √ 3); for GSE universal class,β = 4, Z = 4π/(729 √ 3) (In fact, these are labels for the field of representation. See appendices for more details). Practically we often change r to log r, and the new distribution after the transformation is P (log r) = rp(r). Standard Wigner surmises are shown in the Figure 2. [51] has computed the nearest-neighbor level spacing distribution of 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 investigation shows 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 supersymmetric SYK models, as has been predicted in the previous discussions.
Some comments should be given in this prediction. Firstly, one have some subtleties in obtaining correct rs. 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, 4 in SYK, the particle hole operator P maps each sector to itself, thus if we take all r n the distribution will be ruined, serving as a many-body-localized distribution (the Poisson distribution). For N mod 8 = 2, 6 in SYK, the particle hole operator P maps even and odd parities to each other, and one can take all possible rs in the distribution because all fermionic parity sectors are degenerated. Similar things are observed for all even N in the supersymmetric SYK model. As we mentioned before, the reason is that the supercharge Q is a symmetry of 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, where the correspondence works only for N ≥ 5, because there is no thermalization if N is too small [51]. However, the threshold for obtaining a standard random matrix from N = 1 supersymmetric extension is much larger.
In Section 3.2.1, we argued that the supercharge operator Q in N = 1 supersymmetric SYK theory are also random matrices in some extended ensembles [2,3]. We compute the level statistics of Q and compare it with the Wigner surmises of three standard Dyson ensembles in 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 confirms the relationship between Q's random matrix ensemble and that of the corresponding H.  see all the information in the ensembles.

Spectral form factors
Before presenting the numeric results of spectral form factors, we will review the discreteness of spectrum and the spectral form factor following [54]. For a quantum mechanical system, the partition function could be continued as The analytically continued partition function Z(β, t) is an important quantity to understand a discrete energy spectrum. Typically, people will compute the time average to understand the late time behavior, but for Z(β, t), it vibrates near zero at late time and the time average should be zero. Thus, we often compute Z(β,t) . For a discrete energy eigenvalue spectrum, we have It's hard to say anything general directly for a general spectrum, but one can use the long-term average for large enough T (n E means the degeneracy). For a non-degenerated spectrum, it should have a simple formula However, for a continuous spectrum, the quantity has vanishing long-term average. Thus, the quantity should be an important criterion to detect the discreteness. In this paper, we will use a similar quantity, which is called the spectral form factor In the SYK model, these quantities will have similar predictions with the hamiltonian replaced by random matrix from some specific given Dyson ensembles. For example, for a given realization M from a random matrix ensemble with large L, we have the analytically continued partition function The properties of spectral form factors given by random matrix theory, g rmt (t), have been studied in [54]. There are three specific periods in g rmt (t). In the first period, the spectral form 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 pattern is extremely similar with SYK model. Theoretically, in the early time (before t d ), g(t) should not obtained by g rmt (t) because of different initial dependence on energy, while in the late time these two systems are conjectured to be coincide [54].
With the data of energy eigenvalues one could compute the spectral form factors, which have been shown in Figure 5 for supersymmetric SYK model. We perform the calculation for three different functions g(t), g d (t) and g c (t) with β = 0, 5, 10 and several N s. Clear patterns similar with random matrix theory predictions are shown in these numerical simulations. One could directly see the dip, ramp and plateau periods. For small βs there exist some small vibrations in the early time, while for large β this effect disappears. The function g d is strongly vibrating because we have only finite number of samples. One could believe that the infinite number of samples will cancel the noisy randomness of the curves.
A clear eight-fold correspondence has been shown in the spectral form factor. Near the plateau time of g(t) one should expect roughly a smooth corner for GOE-type, a kink for GUE-type, and a sharp peak for GSE-type. Thus, we observe roughly the smooth corners for N = 14, 16, 22, 24, while the sharp peaks for N = 18, 20, 26, 28 (although the peaks look not very clear because of finite sample size). For N = 10, 12, as shown in Figure 3 there is  no clear random matrix correspondence because N is too small, thus we only observe some vibrations 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 the statistics of the corresponding H. The spectral form factors of Q are slightly different from those of H, yet they show exactly the same eight fold behavior.

Dip time, plateau time and plateau height
More quantitative data could be read off from the spectral form factors. In Figure 7, Figure  8 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 linear fitting 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 the error bar could be computed as the standard deviation.
It is claimed in [54] that polynomial and exponential fitting could be used to interpret the dip time as a function of N with fixed temperature. We apply the same method to the supersymmetric extension. However, we find that in the supersymmetric extension, the fitting 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 polynomial fitting, the fitting effect of exponential function is relatively better. On the exponential fittings with respect to different degeneracy groups, the coefficients before N are roughly the same (0.24N for β = 5) while the overall constants are different. That indicates that the eight-fold degeneracy class or random matrix class might influence the overall factors in the dip time exponential expressions.
One could also read off the plateau time and exponentially fit the data. Similar with dip time, we could also separately fit the plateau time with respect to two different random matrix classes, and one could find a difference in the overall factors of these two groups, while the coefficients before N are the same. There is a non-trivial check here. Theoretically from random matrix theory one can predict that the plateau time scales like t p ∼ e S(2β) [54]. In the large β limit, the entropy should be roughly the ground state entropy. Analytically, the entropy is predicted by S(β = ∞) = N s 0 = 0.275N . Now check the largest β we take (β = 10), we can read off the entropy by 0.277N (GSE-type), 0.275N (GOE-type), or 0.277N (two groups together), which perfectly matches our expectation.  For the plateau height, one can clearly see an eight-fold structure. From the previous discussion we obtain that the plateau height should equals to Z(2β)/Z(β) 2 times a contribution from the degeneracy, which is clearly shown in the figure. For N = 14, 16, 22, 24 (GOE-type), the degeneracy is two thus points should be on the lower line, while for N = 18, 20, 26, 28 (GSE-type), the degeneracy is four thus points should be on the upper line. These observations match the prediction from random matrix theories.

Conclusion and outlook
In this paper, we use analytic arguments and numerical evidence to explore the supersymmetric constraints on the random matrix theory symmetry class. We focus on the N = 1 supersymmetric SYK model, a supersymmetric generalization of nonlocal-coupled majonara fermions 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 supersymmetric SYK model has a different behavior for N mod 8 structure. These arguments might be made to be more general: supersymmetry could directly change the universal class of Hamiltonian (GOE/GUE/GSE) by classifying the symmetry class of supercharge, where combinations 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 type of distribution from Gaussian to Wishart-Laguerre. These points may happen for generic supersymmetric statistical physics models.
We also use numerical method, exact diagonalization to confirm the random matrix theory classification 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 feature from one-point function of extended random matrix theory ensembles, while the Hamiltonian shows 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 capture important physical features, and the new eight-fold rule could be verified.
Several future directions could be investigated. Firstly, one may consider higher supersymmetry constraints on the SYK model, such as N = 2 generalization. Many thermodynamical and field theory properties of higher supersymmetric SYK theory are non-trivial, and it might be interesting to connect these properties to random matrix theory. Moreover, to understand the spectral form factor with supersymmetric constraints, one could also try to study superconformal field theory partition functions at late time. Finally, introducing supersymmetries in the symmetry classification of phases in the condensed matter theory will bring more understanding at the boundary of condensed matter and high energy physics. We leave these interesting possibilities to future works. metric SYK models. TL, JL, YX and YZ are supported by graduate student programs of University 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]) that brings hamiltonians to ten different random matrix classes. In a physical system, symmtries can appear and they consist a group G, then the space of physical states is a projective representation of the symmetry group. A fundamental question we can ask is, what is the most general 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 homomorphism from group G to Z 2 which labels unitarity of operators. Let G 0 be the subgroup of unitary operators, then V splitts into irreps of G 0 : where V i are irreps and m i are their multiplicities in V . If there is no antiunitary operators then followed by Schur's lemma, the most general Hamitonians are those belong to the set plus Hermicity. This is called Type A in the Altland-Zirnbauer theory, without any antiunitary 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 → T U T −1 , is an automorphism of G 0 , thus T maps a component V i ⊗ C m i to another V j ⊗ C m j . A simple case is when i = j, which is easy to see that the most general hamiltonian is of form [2,3] (H, T HT −1 ) (A. 3) where H is an Hermitian operator in component i and T HT −1 acts on component j. Thus it's also of Type A.
The Type A is the simplest structure without any further symmetries. However, if we consider i = j, and consider more anti-unitary operators, the situation is much more technical. It turns out that possible hamiltonians with specific symmetric structures can be classified into ten classes. Here we skip the detailed analysis and directly present the final results. These classes are classified by the following three different operators, • T + , antiunitary, commutes with hamiltonian, and T 2 + = ±1 • T − , antiunitary, anticommutes with hamiltonian, and T 2 − = ±1 • Λ, unitary, anticommutes with hamiltonian, and Λ 2 = 1 If two of these three operators exist, the third will be determined by the following identity, The properties of these three operators can classify the hamiltonian in the following ten classes, where there are no values in some corresponding operators we mean that there is no such a symmetry in the system. We also present the block representation in this table, where blocks are classified by the ±1 eigenspace of anti-unitary operators. The first three ensembles in this table are original Dyson ensembles, while other extended ensembles are their subsets with higher symmetries. These classifications are widely used in theoretical physics, for example, the symmetry classifications of topological insulators and topological phases [4,5].

B Eigenvalue distribution
This appendix is a simple introduction on the random matrix theory eigenvalue distribution (for instance, see [73,74] We 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 to Wishart-Laguerre, which is P (λ)dλ = C (N,β,α)|∆(λ)|β