Complete random matrix classification of SYK models with $\mathcal{N}=0$, $1$ and $2$ supersymmetry

We present a complete symmetry classification of the Sachdev-Ye-Kitaev (SYK) model with $\mathcal{N}=0$, $1$ and $2$ supersymmetry (SUSY) on the basis of the Altland-Zirnbauer scheme in random matrix theory (RMT). For $\mathcal{N}=0$ and $1$ we consider generic $q$-body interactions in the Hamiltonian and find RMT classes that were not present in earlier classifications of the same model with $q=4$. We numerically establish quantitative agreement between the distributions of the smallest energy levels in the $\mathcal{N}=1$ SYK model and RMT. Furthermore, we delineate the distinctive structure of the $\mathcal{N}=2$ SYK model and provide its complete symmetry classification based on RMT for all eigenspaces of the fermion number operator. We corroborate our classification by detailed numerical comparisons with RMT and thus establish the presence of quantum chaotic dynamics in the $\mathcal{N}=2$ SYK model. We also introduce a new SYK-like model without SUSY that exhibits hybrid properties of the $\mathcal{N}=1$ and $\mathcal{N}=2$ SYK models and uncover its rich structure both analytically and numerically.


Introduction
Understanding the mechanism of thermalization and information spreading (scrambling) in nonequilibrium quantum many-body systems is one of the fundamental challenges in theoretical physics. In a classically chaotic system the information on the initial conditions is quickly lost, which can be measured by the Lyapunov exponent that characterizes the sensitivity of the orbit to perturbations of initial conditions. In quantum systems, one clear fingerprint of chaos is the fact that statistical properties of the energy levels are given by random matrix theory (RMT) [1][2][3]. Quantum chaos in this sense has been the subject of research over decades, and its possible role in the relaxation (or thermalization) of a quantum system to equilibrium is still actively debated [4][5][6][7][8][9][10][11]. Further progress was made on the treatment of black holes and holography in terms of quantum information theory [12][13][14][15][16][17][18][19][20][21]. Building on these works, Kitaev suggested to employ the so-called out-of-timeordered correlator (OTOC) [22] to probe information scrambling in black holes and in more general quantum systems [23]. Along this line of thought one can define a quantum analog of the classical Lyapunov exponent, which is argued to have an intrinsic upper bound under certain assumptions [24]. Based on earlier work of Sachdev and Ye [25], Kitaev put forward a (0 +1)-dimensional fermionic model with all-to-all random interactions that can be solved in the large-N limit, with N the number of fermions involved [26]. While it is hard to avoid the spin-glass phase at low temperatures in the original Sachdev-Ye model [27,28], it is ingeniously avoided in Kitaev's model, where fermions are put on a single site. Despite its apparent simplicity, this new Sachdev-Ye-Kitaev (SYK) model has a number of intriguing properties, including the spontaneous breaking of reparametrization invariance, emergent conformality at low energy, and maximal quantum chaos at strong coupling that points to an underlying duality to a black hole [26,[29][30][31][32][33]]. Since the model was announced, a variety of generalizations appeared and computations of the OTOC in various other models were performed [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49], including an SYK-like tensor model without random disorder [50], modified SYK models with a tunable quantum phase transition to a nonchaotic phase [51][52][53], and supersymmetric generalizations of the SYK model [54] (see also [55][56][57][58][59]). An analysis of tractable SYK-type models with SUSY will not only help to better understand theoretical underpinnings of the original AdS/CFT correspondence [60] but also provide insights into condensed matter systems with emergent SUSY at low energy [61][62][63][64][65].
The level statistics of the SYK model was numerically examined in [66][67][68][69] via exact diagonalization and agreement with RMT was found (although sizable discrepancies from RMT were seen in the long-range correlation [67]). An intimate connection between the SYK model and the so-called k-body embedded ensembles of random matrices [70,71] was also pointed out [67]. The algebraic symmetry classification of the SYK model based on RMT in [66][67][68] was recently generalized to the N = 1 supersymmetric SYK model [72]. A random matrix analysis of tensor models has also appeared [73,74].
In this paper, we complete the random matrix analysis of the SYK model. Specifically: 1. We extend the symmetry classification of SYK models with N = 0 and 1 SUSY that were focused on the 4-body interaction Hamiltonian [66][67][68] to generic q-body interactions. The correctness of our classification is then checked by detailed numerical simulations of the SYK model.
2. We provide a detailed numerical examination of the hard-edge universality of energylevel fluctuations near zero in SYK models. 3. We delineate the complex structure of the Hilbert space of the N = 2 SYK model and provide a complete random matrix classification of energy-level statistics in each eigenspace of the fermion number operator.
This paper is organized as follows. In section 2 we review the random matrix classification of generic Hamiltonians to make this paper self-contained. In section 3 we study the non-supersymmetric SYK model. We determine the relevant symmetry classes and report on detailed numerical verifications. In section 4 we study the N = 1 supersymmetric SYK model in a similar fashion. In section 5 we introduce a new SYK-like model that shares some properties (e.g., numerous zero-energy ground states) with the N = 2 SYK model but is theoretically much simpler. In section 6 we investigate the N = 2 supersymmetric SYK model. We explain why the symmetry classification of this model is far more complex than for its N = 1 and 0 cousins. We identify random matrix ensembles for each eigenspace of the fermion number operator and present a quantitative comparison between the level statistics of the model and RMT by exact diagonalization. Section 7 is devoted to a summary and conclusions.
Throughout this paper, we will denote the number of Majorana fermions by N m and the number of complex fermions by N c . The number of fermions in the Hamiltonian is denoted by q and that in the supercharge is denoted byq. Needless to say, q is even andq is odd.
In the early days of RMT, there were just 3 symmetry classes called the Wigner-Dyson ensembles, which can be classified by the presence or absence of the time-reversal invariance and the spin-rotational invariance of the Hamiltonian [87][88][89][90]. It is convenient to distinguish them by the so-called Dyson index β, which counts the number of degrees of freedom per matrix element in the corresponding random matrix ensembles: β = 1, 2, and 4 corresponds, respectively, to the Gaussian Orthogonal Ensemble (GOE), the Gaussian Unitary Ensemble (GUE), and the Gaussian Symplectic Ensemble (GSE). By diagonalizing a random matrix drawn from each ensemble, one finds the joint probability density for all eigenvalues {λ n } to be of the form P (λ) ∝ i<j |λ i − λ j | β n e −V (λn) , where V (x) ∝ x 2 is a Gaussian potential. The spectral density R(λ), also called the one-point function, measures the number of levels in a given interval [λ, λ + dλ]. In RMT, one can show under mild assumptions that for large matrix dimension this function approaches a semicircle R(λ) ∝ √ Λ 2 − λ 2 (Wigner's semicircle law), but in real physical systems R(λ) is typically sensitive to the microscopic details of the Hamiltonian, and one cannot exactly match R(λ) in RMT with the physical spectral density. By contrast, if one looks into level correlations after "unfolding", which locally normalizes the level density to 1, one encounters universal agreement of physical short-range spectral correlations with RMT. 1 Heuristically, larger β implies stronger level repulsion and a more rigid spectrum. A quantum harmonic oscillator exhibits a spectrum with strictly equal spacings, while a completely random point process allows two levels to come arbitrarily close to each other with nonzero probability. RMT predicts a nontrivial behavior that falls in between these two extremes. It is well known that a quantum system whose classical limit is chaotic tends to exhibit energy-level statistics well described by RMT [1,7,8]. Also, Wigner-Dyson statistics emerges in mesoscopic systems with disorder, where the theoretical understanding was achieved by Efetov [91].
An important property of the β = 4 class is the Kramers degeneracy of levels. In general, when there is an antiunitary operator P acting on the Hilbert space that commutes with the Hamiltonian, P −1 HP = H, it follows that for each eigenstate ψ there is another state P ψ that has the same energy as ψ. If P 2 = 1 (GOE), P ψ is not necessarily linearly independent of ψ, hence levels are not degenerate in general, whereas if P 2 = −1 (GSE) their linear independence can be readily shown, so that all levels must be twofold degenerate. We note that the existence of such an operator is a sufficient, but not necessary, condition for the degeneracy of eigenvalues.
Long after the early work by Wigner and Dyson, 7 new symmetry classes were identified in physics. Hence there are now 10 classes in total. (Some authors count them as 12 by distinguishing subclasses more carefully, as we will describe later.) The salient feature pertinent to those post-Dyson classes is a spectral mirror symmetry: the energy levels are symmetric about the origin (also called "hard edge"). This means that, while they show the standard GUE/GOE/GSE level correlations in the bulk of the spectrum (i.e., sufficiently far away from the edges of the energy band), their level density exhibits a universal shape near the origin, different for each symmetry class. (Such a property is absent in the Wigner-Dyson classes since the spectrum is translationally invariant after unfolding and there is no special point in the spectrum.) The physical significance of such near-zero eigenvalues depends on the specific context in which RMT is used. In Quantum Chromodynamics (QCD), small eigenvalues of the Dirac operator in Euclidean spacetime are intimately connected to the spontaneous breaking of chiral symmetry and the origin of mass [79,92,93]. In mesoscopic systems that are in proximity to superconductors, small energy levels describe low-energy quasiparticles and hence affect transport properties of the system at low temperatures. In supersymmetric theories the minimal energy is nonnegative, and it takes a positive value when SUSY is spontaneously broken [94][95][96].
The three chiral ensembles [79,[97][98][99][100] relevant to systems with Dirac fermions such as QCD and graphene are denoted by chGUE/chGOE/chGSE (also known as the Wishart-Laguerre ensembles) and have the block structure 0 * * 0 , which anticommutes with the chirality operator γ 5 = 1 0 0 −1 . This accounts for the spectral mirror symmetry in these 3 1 A cautionary remark is in order. When unitary symmetries are present, the Hamiltonian can be transformed to a block-diagonal form, where each block is statistically independent. The spectral correlations must then be measured in each independent block. If one sloppily mixes up all eigenvalues before measuring the spectral correlations, the outcome is just Poisson statistics (see section III.B.5 of [2] for a detailed discussion).
classes. We remark that chiral symmetry (i.e., a unitary operation that anticommutes with the Hamiltonian) is often called a sublattice symmetry in the condensed matter literature. A unique characteristic of the chiral classes in contrast to the other 4 mirror-symmetric classes is that there can be an arbitrary number of exact zero modes. This is easily seen by making the matrix block * rectangular, say, of size m × n. When |m − n| is large, the nonzero levels are pushed away from the origin due to level repulsion. In the limit m, n → ∞ with m/n → 1 the macroscopic spectral density fails to approach Wigner's semicircle and rather converges to what is called the Marčenko-Pastur distribution [101]. In the thermodynamic limit of QCD with nonzero fermion mass, the number of zero modes |m − n| ∝ V 1/2 4 [93] while m, n ∝ V 4 , where V 4 is the Euclidean spacetime volume, and hence the physical limit is m/n → 1.
The other 4 post-Dyson classes are referred to as the Bogoliubov-de Gennes (BdG) ensembles. They were identified by Altland and Zirnbauer [75,102]. It is the particlehole symmetry that accounts for the mirror symmetry of the spectra in these classes. This completes the ten-fold classification of RMT as summarized in table 1. There is a one-to-one correspondence between each ensemble and symmetric spaces in Cartan's classification, so the RMT ensembles are often called by abstract names such as A, AI, and AII due to Cartan [76]. In recent years this classification scheme was found to be useful in the classification of topological quantum materials [86,[103][104][105][106][107].
We refer the reader to [75-77, 108, 109] for the detailed mathematics of the Altland-Zirnbauer theory and only recall the essential ingredients here. Let T + (T − ) denote an antiunitary operator that commutes (anticommutes) with the Hamiltonian. 2 (Note that any antiunitary operator can be expressed as the product of a unitary operator and the complex conjugation operator K.) The chirality operator (a unitary operator that anticommutes with the Hamiltonian and squares to 1) is denoted by Λ from here on. The first step is to check whether T + , T − , and Λ exist for a given Hamiltonian. If both T + and T − exist, one always has chiral symmetry, Λ = T + T − . The second step is to check if the antiunitary symmetry squares to +1 or −1. This allows one to figure out which class the Hamiltonian belongs to. However, there is an additional subtlety in the symmetry classes BD and DIII. There one has to distinguish two cases according to the parity of the dimension of the Hilbert space (see table 1), which results in the presence/absence of exact zero modes. The classes B and DIII-odd have physical applications to superconductors with p-wave pairing [110][111][112][113]. The functional forms of the universal level density near zero for all the 7 post-Dyson classes are explicitly tabulated in, e.g., [113,114]. Note that, because class B and class C share the same set of indices α and β, their level density near zero coincides, except for a delta function at the origin in class B. 2 Here we conform to the notation of [72]. Rather than calling T± time-reversal symmetry or spinrotational symmetry, we prefer to denote them by abstract symbols, since the proper physical interpretation of each operator depends on the specific system. Table 1. Classification of RMT symmetry classes. In the first three rows we list the Wigner-Dyson classes. β is the Dyson index defined in the main text. In the remaining rows we list the chiral and BdG classes. The joint probability density for energy levels in these ensembles assumes the form P (λ) ∝ i<j |λ 2 i − λ 2 j | β n |λ n | α , and the indices β and α are presented in the third and fourth column, respectively. α is related to the number of exact zero modes. The index ν defined in the last column is related to the topological charge of the gauge field in non-Abelian gauge theories.
Here we define ν to be nonnegative. The symbol "-" implies that there is no symmetry in that class. The classes B and DIII-odd are sometimes omitted in other references, but we include them here for completeness. T + (T − ) denotes an antiunitary operator that commutes (anticommutes) with the Hamiltonian, and Λ is the chirality operator. If both T + and T − are present, there is chiral symmetry, but the converse is not true in general. Our notation in this table is such that A is the complex conjugate of A and A † is the conjugate transpose of A, i.e.,

Definitions of relevant operators
To begin with, recall that when we speak of a non-SUSY SYK model, there are actually two models, one involving Majorana fermions [26,30,31] and another involving complex fermions [29,35,44,66,115]. In either case it is useful to start with the creation and annihilation operators of complex fermions, denoted by c a and c a , respectively, obeying These operators can be represented as real matrices by adopting the Jordan-Wigner construction [35,66] c a = ( 1≤b<a σ z b )(σ x a + iσ y a )/2 and c a = (c a ) † . 3 We also define the fermion number operator (3. 2) The total Hilbert space V of dimension 2 Nc splits into two sectors with even/odd eigenvalue of F , i.e., (−1) F = ±1.
One can construct N m = 2N c Majorana fermions χ i from complex fermions as The antiunitary operator of special importance in the SYK model is the particle-hole operator [66][67][68]116] where K is complex conjugation. One can show [66][67][68] P c a P = ηc a , P c a P = ηc a , P χ i P = ηχ i , (3.5) Here x denotes the greatest integer that does not exceed x. We stress that all of the above formulas hold irrespective of the form of the Hamiltonian.

Classification
Let us begin with the non-supersymmetric SYK model with N m Majorana fermions for N m even. 4 For a positive even integer 2 ≤ q ≤ N m , the Hamiltonian [26,30,31] is given by where J i 1 ···iq are independent real Gaussian random variables with the dimension of energy, J i 1 ···iq = 0 and J 2 The prefactor i q/2 is necessary to make H Hermitian. This model is conjectured to be dual to a black hole in the large-N limit [26,30,31] and for βJ 1 saturates the bound on quantum chaos proposed in [24]. While the q = 4 version has attracted most of the attention in the literature, it is useful to consider general q because the theory simplifies in the large-q limit [26,31]. Now, due to the Majorana nature of the fermions, the fermion number is only conserved modulo 2. The Hilbert space naturally admits a decomposition into two sectors of equal dimensions, with a definitive parity of the fermion number. Since H does not mix sectors with (−1) F = +1 and −1, H acquires a block-diagonal form Hermitian square matrices of equal dimensions. By examining the commutation relation of H and P , one finds that q = 0 (mod 4) and q = 2 (mod 4) have to be treated separately because HP = (−1) q/2 P H. The spectral statistics for q = 2 (mod 4) did not receive attention in [66][67][68]72], 5 and we shall work it out below. This is a new result. q = 0 (mod 4) In this case [H, P ] = 0. Thus P corresponds to T + in table 1. For N m = 0 and 4 (mod 8), P is a bosonic operator and maps each parity sector onto itself. For N m = 0 (mod 8), P 2 = +1 so that H = GOE⊕GOE. For N m = 4 (mod 8), P 2 = −1 so that H = GSE⊕GSE. In both cases the two blocks of H are independent in general. Finally, for N m = 2 and 6 (mod 8) P is a fermionic operator and exchanges the two sectors. Hence H = belongs to GUE. It follows that the eigenvalues are twofold degenerate for N m = 2, 4 and 6 (mod 8), and unpaired only for N m = 0 (mod 8). This is summarized in table 2, which is consistent with [66][67][68]72]. q = 2 (mod 4) Now {H, P } = 0. Thus P corresponds to T − in table 1 and the spectrum enjoys a mirror symmetry λ ↔ −λ. 6 For N m = 0 and 4 (mod 8), P is a bosonic operator and maps each parity sector onto itself. For N m = 0 (mod 8), P 2 = +1 so that H = BdG(D) ⊕ BdG(D). (It is not class B because the dimension 2 Nm/2−1 of each sector is even.) For N m = 4 (mod 8), P 2 = −1 so that H = BdG(C) ⊕ BdG(C). In both cases the two blocks of H are independent in general. For N m = 2 and 6 (mod to GUE, for the same reason as above. This is summarized in table 3. As a generalization one can also consider a Hamiltonian that includes both a q = 0 (mod 4) term and a q = 2 (mod 4) term. Then H has no antiunitary symmetry and the result is just GUE ⊕ GUE, i.e., H = A 0 0 B with A and B independent Hermitian matrices. 5 An exception is the simplest case q = 2, which was analytically solved at finite Nm [56] and in the limit Nm → ∞ [31,55] (see also [68,117]). Note that H in this theory is just a random mass with no interactions, so one cannot extrapolate features of q = 2 to the more nontrivial q ≥ 4 cases. 6 What is meant here is that the mirror symmetry is present for every single realization {Ji 1 ,··· ,iq } of the disorder.  [66][67][68]72].
Even when the symmetry class of H is known, it is highly nontrivial whether the level correlations of H quantitatively coincide with those of RMT. In the SYK model (3.7) there are only O(N q m ) independent random couplings, while a dense random matrix has O(2 Nm ) independent random elements. The level statistics of H for q = 4 has been studied numerically via exact diagonalization [66][67][68][69] and agreement with the RMT classes in table 2 was found for not too small N m . This is consistent with the quantum chaotic behavior of the model [26,31].

Numerical simulations Level correlations in the bulk
Here we report on the first numerical analysis of the bulk statistics of energy levels for the N = 0 SYK model with q = 6 via exact diagonalization to test table 3. To identify the symmetry class we employ the probability distribution P (r) of the ratio r = (λ n+2 − λ n+1 )/(λ n+1 − λ n ) of two consecutive level spacings in a sorted spectrum, as it does not require an unfolding procedure [66,118,119]. We used accurate Wigner-like surmises for the Wigner-Dyson classes derived in [119], For Poisson statistics we have P P (r) = 1/(1 + r) 2 [119]. Our numerical results are displayed in figure 1. Without any fitting parameter, they all agree excellently with the GUE (β = 2) as predicted by table 3. This indicates that quantum chaotic dynamics emerges in this model even for such small values of N m .

Universality at the hard edge
In class C and D the origin is a special point due to the spectral mirror symmetry, and the level statistics near zero shows universal fluctuations different from those in the bulk of the spectrum [75]. Their form is solely determined by the global symmetries of the Hamiltonian and is insensitive to microscopic details of interactions. In figure 2 we compare the distributions of the near-zero energy levels of the N = 0 SYK model with q = 6 and those of RMT, finding nearly perfect agreement. 7 The nonzero (zero) intercept at λ = 0 in class D (class C) directly reflects the fact that α = 0 for class D (α = 2 for class C), where α is the index listed in table 1.

Overview of the N = 0 SYK model with complex fermions
We finally comment on the non-supersymmetric SYK model with complex fermions [29,35,44,66,115]. The Hamiltonian reads where µ is the chemical potential for the fermion number operator F in (3.2) and the coupling is a complex Gaussian random variable obeying Since H preserves the fermion number, H as a matrix has a block-diagonal structure representing each eigenspace of F = 0, 1, . . . , N c . There is no antiunitary symmetry for H and consequently the levels collected in each block of H would obey GUE. Intriguingly, one can amend H by adding correction terms so that it commutes with P [35,66]. In this case, the half-filled sector 7 To obtain these plots we determined the RMT curves numerically for matrix size 10 3 using the mapping to tridiagonal matrices invented in [120]. We then rescaled the RMT curves as p(x) → cp(cx) and tuned the parameter c to achieve the best fit to the data, where c is common to the three curves in each plot. 4 N = 1 SYK model

Classification
The supersymmetric generalization of the SYK model was introduced in [54] (see also [55][56][57][58][59]). The model with N = 1 SUSY has the Hamiltonian H = Q 2 with supercharge where 1 ≤q ≤ N m is an odd integer. (Note that Q † = Q.) In this case H involves terms with up to 2q − 2 fermions. The couplings C i 1 i 2 ···iq are independent real Gaussian variables with mean C i 1 i 2 ···iq = 0 and variance C 2 The groundstate energy of this model is evidently nonnegative. In [54] a strictly positive ground-state energy that decreases exponentially with N was obtained numerically, indicating that SUSY is dynamically broken at finite N and restored only in the large-N limit.
It is easy to verify the simple relation between the spectral densities of H and Q, where ρ H (λ) ≡ Tr δ(λ − H) and ρ Q (X) ≡ Tr δ(X − Q) . Equation (4.2) reveals that the level density of H would blow up as λ −1/2 near zero if Q had a nonzero density of states at the origin for large N m . This blow-up was indeed seen in the exact diagonalization analysis [72] as well as in analytical studies of the low-energy Schwarzian theory [54,121,122]. Since Q is more fundamental than H we will Table 4. Symmetry classification of Q in the N = 1 SYK model forq = 1 (mod 4). For the block structure of each class we refer to table 1.
focus on the level structure of Q below, viewing it as a matrix acting on the many-body Fock space.
The random matrix classification forq = 3 has recently been put forward in [72]. Here we will generalize this to all oddq, with emphasis on the difference of symmetry classes betweenq = 1 (mod 4) andq = 3 (mod 4). The main theoretical novelty in the N = 1 SYK model is the fact that Q anticommutes with the fermion parity operator (−1) F . Thus (−1) F plays the role of γ 5 for the Dirac operator in QCD and naturally induces a block structure 0 * * † 0 for Q. The spectrum of Q is therefore symmetric under λ ↔ −λ. Since the block * is a square matrix, there are no topological zero modes, i.e., all eigenvalues of Q are nonzero unless fine-tuning of the matrix elements is performed. From the relation H = Q 2 we conclude that all eigenvalues of H should be at least twofold degenerate.
Following [72] we introduce a new antiunitary operator R ≡ P (−1) F . We have where N c = N m /2 as before and η is given in (3.6). These relations, combined with table 1, lead to the classification of Q shown in table 4 forq = 1 (mod 4) and table 5 forq = 3 (mod 4). By comparing the (anti-)commutators in each table, we see that the roles of P and R are exchanged forq = 1 and 3. Consequently the positions of BdG(CI) and BdG(DIIIeven) are exchanged. In these tables we made it clear that we are considering chGOE and chGSE in the topologically trivial sector ν = 0.
One can also consider a superposition of multiple fermionic operators in the supercharge, e.g, real Gaussian couplings. Then Q fails to commute or anti-commute with P and R and the symmetry class is changed: Q now belongs to the β = 2 chGUE (AIII) class with ν = 0. There is no degeneracy of eigenvalues for Q while all eigenvalues of H = Q 2 are two-fold degenerate since {(−1) F , Q} = 0.
In all cases considered above for N = 1, the symmetry classes differ from the Wigner-Dyson classes because of the presence of chiral symmetry (−1) F . This difference manifests itself in distinctive level correlations near the origin (universality at the hard edge). In Table 5. Symmetry classification of Q in the N = 1 SYK model forq = 3 (mod 4). This table is consistent with [72]. For the block structure of each class we refer to table 1.
order to expose this in the thermal N = 1 SYK model, the temperature must be lowered to the scale of the smallest eigenvalue of H. This is exponentially small in N m .

Numerical simulations
Level correlations in the bulk Previously, the level statistics in the bulk of the energy spectrum for the N = 1 SYK model withq = 3 was studied in [72] and results consistent with table 5 were reported. Here we report the first numerical analysis of the bulk statistics for the N = 1 SYK model witĥ q = 5 via exact diagonalization, to test table 4. To identify the symmetry class, we again used the ratio of two consecutive level spacings. Our numerical results are displayed in figure 3. Excellent agreement with the RMT curves of the symmetry classes predicted by table 4 is observed. This evidences the existence of quantum chaotic dynamics in this model and corroborates our classification scheme.  The smallest eigenvalue approaches zero from above for larger N m , indicating restoration of SUSY in the large-N m limit as already reported in [54]. We note that the RMT classes chGOE (BDI) and chGSE (CII) were originally invented and exploited in attempts to theoretically understand fluctuations of small eigenvalues of the Euclidean QCD Dirac operator with special antiunitary symmetries in a finite volume [99,100,[123][124][125], 8 related to spontaneous breaking of chiral symmetry through the Banks-Casher relation [92]. The RMT predictions agree well with the Dirac spectra taken from lattice QCD simulations [131]. It is a nontrivial observation that the smallest energy levels of the N = 1 SYK model, which set the scale for the spontaneous breaking of SUSY, obey the same statistics as the eigenvalues of the Dirac operator in QCD, which has totally different microscopic interactions compared to the SYK model. This is yet another example for random matrix universality.
5 Interlude: a simple model bridging the gap between N = 1 and 2

Motivation and definition
The SYK model with N = 2 SUSY [54] has the Hamiltonian H = {Q, Q} with two supercharges Q and Q, each comprising an odd number of complex fermions. This model preserves the U(1) fermion number exactly, so that the Hamiltonian is block-diagonal in the fermion-number eigenbasis. As shown by the Witten-index computation in [54], the Hamiltonian has an extensive number of exact zero modes 9 and SUSY is unbroken at finite N c . These features are in marked contrast to the N = 1 SYK model, where the fermion number is only conserved modulo 2, the Hamiltonian is positive definite with no exact zero modes, and SUSY is spontaneously broken at finite N c .
While there is no logical obstacle to moving from N = 1 to 2, it is helpful to have a simple model that serves as a bridge between these two theories. The model we designed for this purpose is defined by the Hamiltonian H = M 2 with the Hermitian operator where 1 ≤ p ≤ N c is an even integer and Z j 1 ···jp are independent complex Gaussian random variables with mean zero and Z ab Z ab = 2J/N 2 c for some J > 0. The creation and annihilation operators c a and c a were introduced in section 3.1. Because of M = M † we have H ≥ 0, similarly to the supersymmetric SYK models. If we forcefully substitute p = 3 and let i p/2 → i, then M = Q + Q and H = M 2 = {Q, Q}, i.e., the N = 2 SYK model is recovered (see section 6). What difference emerges if we retain an even number of 8 See also [126][127][128][129][130] for related works in mathematics. 9 The existence of a macroscopic number of ground states is a familiar phenomenon in lattice models with exact SUSY [132][133][134][135][136][137]. fermions in M ? Of course it makes M a bosonic operator and destroys SUSY. At this cost, however, we gain three new features that were missing in the N = 1 SYK model: (i) the fermion number is conserved modulo 2p (rather than modulo 2), (ii) H has a large number of exact zero modes, and (iii) an interplay between N c and F emerges in the symmetry classification of energy-level statistics. The last point is especially intriguing since this property is shared by the N = 2 SYK model (section 6). This is why we regard this model as "intermediate" between the N = 1 and N = 2 SYK models. Studying the level structure of this exotic model provides a useful digression before tackling the N = 2 case.
By exact diagonalization we have numerically computed the spectral density of M for p = 2 and 4, see figure 6. In all plots there is a delta function at zero due to the macroscopic number of zero-energy states. Interestingly, the global shape never resembles Wigner's semicircle but rather depends sensitively on both p and N c . For p = 2 we observe oscillations in the middle of the spectrum, for which we currently do not have a simple explanation. The case p = 2 could be more the exception than the rule, 10 much like the q = 2 SYK model that is solvable and nonchaotic [31,56,68] unlike its q > 2 counterparts.
For both p = 2 and 4, a close inspection of the plots near the origin reveals that for odd N c there is a dip of the density around the origin, indicating that small nonzero levels are 10 We speculate that the spectral density for this case may even be computed exactly since M is just a fermion bilinear, but this is beyond the scope of this paper.
repelled from the origin, while there is no such repulsion for even N c . The same tendency of the spectral density (albeit with the parity of N c reversed) has been observed for the N = 2 SYK model, too [121]. We will give a simple explanation of this phenomenon later.

Classification for p = 2
To make the presentation as simple as possible, we shall begin with p = 2, in which case the fermion number F is conserved modulo 4. The Hilbert space V of N c complex fermions can be arranged into a direct sum of four spaces V 0,1,2,3 , where V f is the eigenspace of F corresponding to F = f (mod 4), i.e., The numbers D 0,1,2,3 are listed for 3 ≤ N c ≤ 10 in table 6. Since there is no nonzero matrix element of M between states with different parity of F we have M = . The chiral structure in each term is due to the chiral symmetry {i F , M } = 0, which ensures the spectral mirror symmetry of M . It should be stressed that A 0 and A 1 are in general rectangular. When they become a square matrix can be read off from table 6. These cases are colored in red and green. They only occur for even N c (which is also true for p = 4, see table 7 below). On the other hand, for odd N c , both A 0 and A 1 are rectangular. As is well known from studies in chiral RMT [79,100], in that case the nonzero eigenvalues of M (i.e., the nonzero singular values of A 0 and A 1 ) are pushed away from the origin by the large number of exact zero modes. Indeed, α in table 1 is proportional to the number of zero modes, and large α suppresses the joint probability density of eigenvalues near zero. This leads to the dip around the origin in the left plots of figure 6. However, for even N c , in the subspaces without exact zero modes there is no repulsion of the nonzero modes from the origin, and thus no dip of the density (which is summed over all subspaces) shows up near zero. In order to understand the level degeneracy in each sector correctly, we must figure out the antiunitary symmetries of the matrix M . We use the particle-hole operator P in (3.4) again. In addition, we define another antiunitary operator S ≡ P · i F . One can show Both P 2 and S 2 are tabulated in table 6, but extra care is needed for S because S 2 is not just ±1 but a nontrivial operator that depends on F . For even N c , each chiral block belongs to one of chGSE (CII) β=4 , BdG (DIII-even) β=4 , chGOE (BDI) β=1 , and BdG (CI) β=1 according to the values of P 2 and S 2 (cf. table 1).
For odd N c , P maps a state in V 0 ⊕ V 2 to V 1 ⊕ V 3 and vice versa. Therefore the nonzero levels of M in V 0 ⊕ V 2 must be degenerate with those in V 1 ⊕ V 3 . Since there is no antiunitary symmetry acting within each chiral block, all uncolored sectors in table 6 belong to chGUE (AIII).
This completes the algebraic classification of the model (5.1) for p = 2 based on RMT. This classification is periodic in N c with period 4 as can be seen from table 6. We have numerically checked the level degeneracy of M in each sector for various N c and confirmed consistency with our classification. In this process we found, surprisingly, that levels often show a large (e.g., 16-fold) degeneracy that cannot be accounted for by our antiunitary symmetries P and S. Such a large degeneracy, which presumably is responsible for the wavy shape in the upper plots of figure 6 and makes the level spacing distribution for p = 2 deviate from RMT, was not observed for p = 4. We interpret this as an indication that the model with p = 2 is just too simple to show quantum chaos and therefore do not investigate it further.

Classification for p = 4
As a more nontrivial case we now study the p = 4 model, which preserves F (mod 8). This time the Hilbert space decomposes as Table 7. Model (5.1) for p = 4. We list the dimensions (5.5) of the eigenspaces of F (mod 8). Uncolored blocks belong to chGUE (AIII) β=2 , while : chGSE (CII) β=4 with ν = |D i − D i+4 |/2, : BdG (DIII-even) β=4 , : chGOE (BDI) β=1 with ν = |D i − D i+4 |, and : BdG (CI) β=1 . Details of each class can be found in table 1. The mark (2) after the number of positive levels of M indicates that those levels are twofold degenerate, e.g., 20 (2) means 10 pairs. In each block of given N c there is an equal number of positive and negative levels because of chiral symmetry, {κ F , M } = 0. Also shown are the squares of the antiunitary operators P and S. The symmetry pattern is periodic in N c with period 8.
, where the terms correspond to V 0 ⊕V 4 , V 1 ⊕V 5 , V 2 ⊕V 6 , and V 3 ⊕V 7 , respectively. As a consequence, the spectrum of M enjoys a mirror symmetry as in the model with p = 2. Let us define an antiunitary operator S ≡ P · κ F , where κ ≡ e iπ/4 is the 8-th root of unity and P was defined in (3.4). One can show The dimension of each subspace of V is listed for 7 ≤ N c ≤ 14 in table 7. As for p = 2, the particle-hole operator P generates degeneracies between distinct chiral blocks. For instance, at N c = 11, the 166 distinct positive levels in V 0 ⊕V 4 are degenerate with those in V 3 ⊕V 7 . The symmetry classification is just a rerun of our arguments for p = 2 and therefore omitted here. We have numerically confirmed that table 7 gives the correct degeneracy of levels. (Unlike for p = 2, we did not observe any unexpected further degeneracies.) Table 7 not only provides a symmetry classification but also enables us to derive a fairly simple analytic approximation to the global spectral density. Let us recall the so-called Marčenko-Pastur law [101]: suppose X is a complex L × N matrix with L ≤ N whose elements are independently and identically distributed with X ij = 0 and |X ij | 2 = σ 2 < ∞. Let us denote the L eigenvalues of √ XX † by {ξ i } ≥ 0. Then for L, N → ∞ with L/N ∈ (0, 1] fixed, the probability distribution of {ξ i } takes on the limit

Global spectral density
This function satisfies the normalization ∞ 0 dx F (α, x) = 1 for all α ∈ (0, 1]. We now exploit this law to describe the global density of our p = 4 model, shown previously in figure 6. Whether (5.7) works quantitatively or not is not obvious a priori because the matrix elements of (5.1) are far from statistically independent, but rather strongly correlated. Putting this worry aside, let us consider the N c = 9 case first. According to table 7, there are four chiral blocks, and two of them are copies of the other two, so we should sum just two Marčenko-Pastur distributions. For N c = 10, we have to sum three. Taking into account that the global density in figure 6 counts both positive modes and exact zero modes, we obtain formulas with the correct normalization, P (p=4,Nc=9) (σ; ξ) = 2 10 · P 10,126 (σ; ξ) + 36 · P 36,84 (σ; ξ) 2 9 − 2(10 + 36) , (5.9a) P (p=4,Nc=10) (σ; ξ) = 2 · 46 · P 46,210 (σ; ξ) + 20 · P 20,252 (σ; ξ) + 120 · P 120,120 (σ; ξ) 2 10 − (2 · 46 + 20 + 120) . (5.9b) The parameter σ has to be tuned to achieve the best fit to the data because RMT does not know the typical energy scale of the model. The results of the fits displayed in the bottom plots of figure 6 show impressive quantitative agreement. We also notice a shortage of levels near the peak density, as well as a leakage of levels toward larger values. Even though the agreement is not perfect it is intriguing that a naïve ansatz such as (5.9) is sufficient to account for the shape of the global density. We tried a similar fit for p = 2 as well but did not find any agreement even at a qualitative level, probably due to the nonchaotic character of the p = 2 model as described before.

Numerical simulations
Level correlations in the bulk We numerically checked the bulk statistics (GOE/GUE/GSE). As there are quite a few chiral blocks in table 7 we did not check all of them but concentrated on three cases: (i) the V 3 ⊕ V 7 sector for N c = 10, (ii) the V 3 ⊕ V 7 sector for N c = 11, and (iii) the V 0 ⊕ V 4 sector for N c = 12. To identify the symmetry classes we again used the probability distribution of the ratio of two consecutive level spacings. Our numerical results are displayed in figure 7, where excellent agreement with the respective symmetry classes predicted by table 7 is found. This corroborates our symmetry classification scheme. 6 N = 2 SYK model

Preliminaries
The N = 2 SYK model [54,58,59] has significantly different properties from its N = 1 cousin. The Hamiltonian is defined by H = {Q, Q} with two supercharges that are nilpotent, Q 2 = Q 2 = 0, where the couplings X ijk are independent complex Gaussian random variables obeying X ijk X ijk = 2J/N 2 c . Apart from the random disorder, this model is somewhat similar to lattice models with exact SUSY [132][133][134][135][136][137]. The model can be generalized so that Q and Q involveq fermions withq odd [54]. We postpone this generic case to section 6.5 and for the moment focus onq = 3, i.e., (6.1). As for the operator P in (3.4), we have As shown in [54,121,122], H possesses a number of exactly zero eigenvalues, so SUSY is not spontaneously broken in contrast to the N = 1 model. Moreover, the N = 2 model has U(1) R-symmetry.
[H, F ] = 0 ensures that H and F can be diagonalized simultaneously. The total Hilbert space V has the structure where V f is the eigenspace of F with eigenvalue f . The level density of H in the low-energy limit has been derived analytically from the large-N c Schwarzian theory [121,122], whereas analysis of the level statistics and symmetry classification of H based on RMT has not yet been done for the N = 2 SYK model. In the remainder of this section we fill this gap.

Naïve approach with partial success
In this subsection we briefly review a simple approach to the N = 2 model that is a natural extrapolation of our treatment for the N = 0 and 1 SYK models but is beset with fatal problems and eventually fails. This subsection is included for pedagogical reasons and can be skipped by a reader interested only in final results. In section 3.4 we have reviewed the symmetry properties of the N = 0 SYK model with complex fermions, which had the virtue of the exactly conserved fermion number, just like the N = 2 SYK model. If one were to boldly extrapolate the statements in section 3.4 to the N = 2 case, one would conclude that the levels of H in all V f except for V Nc/2 belong to GUE while those in V Nc/2 belong to GOE or GSE depending on P 2 = ±1. However, numerical analysis of the level correlations clearly reveals disagreement with the expected statistics. This failure can be traced back to the fact that in this approach all the fine structure of H imposed by N = 2 SUSY is neglected.
So let us change the strategy and try to move along the path we have followed in sections 4 and 5. First of all, note that in the N = 2 SYK model one can write H = M 2 with a Hermitian operator M ≡ Q + Q . Since M preserves F (mod 3) and anticommutes with (−1) F , it is useful to divide V into subspaces V f on which F = f (mod 6), i.e., Closed analytic expressions for D f are given in appendix A. Then M assumes a block- , where the terms correspond to V 0 ⊕ V 3 , V 1 ⊕ V 4 , and V 2 ⊕ V 5 , respectively. The spectrum of M has a mirror symmetry for every single realization of {X ijk }. As a consequence, every nonzero eigenvalue of H is at least twofold degenerate. From the above structure, a lower bound on the number N z of exact zero modes of M and hence of H can readily be obtained (cf. appendix A) as The same bound was obtained via the Witten index in [57]. 11 In numerical simulations we found that this bound is saturated for N c ∈ {0, 2, 3} (mod 4), while a strict inequality holds for N c = 1 (mod 4) due to the presence of O(1) "exceptional" zero modes [54,57] (see also appendix B). We will explain their origin later. We note in passing that the present argument based on M does not tell us how many zero modes exist in each V f .

Global spectral density
Utilizing the decomposition of M into three chiral blocks, we can derive an approximate analytic formula for the global level density based on the Marčenko-Pastur law (5.7), repeating the steps that led to (5.9). (We note that the level densities of M and H are linked by formula (4.2), where Q should be replaced by M here.) Figure 9 displays the numerically obtained global spectral density of M for N c = 9 and 10 together with the analytic approximations obtained by tuning the parameter σ for optimal fits. The quality of the  Figure 9. Spectral density of M = Q + Q in the N = 2 SYK model from exact diagonalization for N c = 9 and 10, averaged over random samples. Since the spectra are symmetric about 0, only the nonnegative part is shown. The delta peaks at the origin represent exact zero modes, as in figure 6. In both plots J = 1 and the total density is normalized to 1. The blue dashed lines are the best fits of analytic approximations based on the Marčenko-Pastur law. 11 We emphasize that the extensive number of zero-energy states in this model owes their existence to the mismatch of D f and D f +3 (f = 0, 1, 2). If one adds an arbitrarily small perturbation that breaks the U(1) R-symmetry down to Z2, the Hamiltonian would lose its triple chiral-block structure and is left with just the two eigenspaces of (−1) F , which have equal dimension. Then nothing protects zero modes from being lifted and SUSY gets broken, as reported in [57,138]. agreement is worse than for the previous model (figure 6). In particular, the pronounced sharp peak of the density cannot be reproduced with the Marčenko-Pastur law. This could be an indication that the N = 2 SYK model indeed has a more complex structure than the model in section 5.
In figure 9 there is a spectral gap for N c = 9 but not for N c = 10. The peculiar dependence of the level density of H on the parity of N c was also noted in [121]. Intriguingly, this can easily be accounted for by the fact that a chiral block with D f = D f +3 is present only for odd N c (cf. appendix A). This can be shown by elementary combinatorics.

Symmetry of M
To classify M based on RMT we can again make use of P and R ≡ P (−1) F in the same way as for the N = 1 SYK model (section 4). For N c = 1 (mod 4), it can easily be shown for the three sets of f and N c specified above. The point is that d odd is an odd integer. This means that the spectrum of M on V f ⊕ V f +3 cannot consist of d odd positive levels and d odd negative levels, since this would contradict the Kramers degeneracy. We conclude that M (and H) must have at least 2 zero modes in V f ⊕ V f +3 . This explains why we encounter "exceptional" zero modes for N c = 1 (mod 4), and is corroborated by our exact diagonalization analysis of H (see appendix B). 12 It turns out, however, that the current approach is incapable of describing the actual level structure of M in full detail. For instance, although M in the sector V 0 ⊕ V 3 with N c = 12 is classified as class chGOE (BDI) β=1 , exact diagonalization shows that all nonzero eigenvalues of M in this sector are in fact twofold degenerate. The reason that the symmetry classification based on M is doomed to be incomplete is that M does not manifestly reflect the fermion-number conservation of H. We have no access to the level statistics in the individual eigenspaces V f of F as long as we see H through the lens of M . The upshot is that since the structure of the N = 2 SYK model is qualitatively different from its cousins with N = 0 and 1 SUSY, we need an entirely new approach to carry out its symmetry classification. This is the subject of the next subsection.
(In figure 11 below we will show a graphical representation of the interrelations of the V ±,z f .) Next we introduce notation for the dimensions of the subspaces,  Up to now we have not mentioned how to compute N ± f explicitly for given f and N c . Actually this proves to be a straightforward (albeit tedious) task if we posit the following premise: For any N c ≥ 3, all exact zero modes of H reside in V f with |f −N c /2| ≤ 3/2, where the equality holds only for exceptional zero modes that occur when N c = 1 (mod 4). 14 (6.15) This rather strong condition on the ground states of H is not only corroborated by detailed numerical simulations (see appendix B and [54]) but also derived from the Schwarzian effective theory valid in the large-N c and low-energy limit [121,122]. If (6.15) is accepted, one can fully clarify the relation of Hilbert spaces linked by Q as in table 8. The sequences tabulated there are exact sequences in the terminology of mathematics, in the sense that the kernel of Q acting on V f coincides exactly with the image of V f −3 by Q. Two examples of these sequences, extended up to V Nc , are graphically illustrated in figure 11 for N c = 15. Table 8. Exact sequences of the Hilbert spaces generated by the linear map Q. Complementary exact sequences descending from V Nc , V Nc−1 , and V Nc−2 by way of Q can be obtained by applying the particle-hole operator P to the sequences in the table. The spaces V contain an exponentially large number of "typical" zero modes, see (6.5). The spaces V * contain no zero modes for N c = 3 (mod 4), or 1 or 3 "exceptional" zero modes for N c = 1 (mod 4).
Although we do not provide a rigorous proof of (6.15), there is a heuristic argument to convince oneself that (6.15) is correct. Let us consider a sequence · · · If Q in the middle were a completely random linear map, it is a matrix of size dim(V f ) × dim(V f +3 ) whose rank is almost surely dim(V f ) (in the absence of fine-tuning or a special symmetry). This is of course an oversimplification for Q, because it is not a generic linear map but a nilpotent map. Taking this into account, let us next view Q as a random matrix of size dim has been left out. Then the rank of Q is almost surely dim V f \ Q(V f −3 ) , i.e., there is no "nontrivial" zero mode of Q in V f . This argument may be repeated along the sequence as long as the condition dim(V f ) < dim(V f +3 ) is fulfilled. A completely parallel argument can also be given for a "descending" sequence · · · . By pinching the sequence from both ends like this, we find at the end of the day that all zero modes (Qψ = Qψ = 0) must be concentrated in the subspace V f with the largest dimension in the sequence. This is equivalent to the condition (6.15). Now it is straightforward to work out N ± f . Let us begin with the case of even N c . First, for 3 ≤ f < N c /2 − 1, V f does not contain zero modes under the assumption (6.15). Hence, with the help of (6.12), we find This recursion relation for N + f is to be solved with the initial conditions N + 0 = 1, N + 1 = N c , and N + 2 = N c (N c − 1)/2. The result reads where f 0 ≡ f − 3 f /3 ∈ {0, 1, 2}. Equation (6.11) was used in the first equalities of (6.17) and (6.18). These formulas hold in the range 0 ≤ f < N c /2 − 1. We verified (6.17) numerically for N c up to 17. Finally, to derive N ± f for f close to N c /2 , we need to know N z f . Recalling the premise (6.15) and the fact that the inequality (6.5) is saturated except when N c = 1 (mod 4) (see appendix A for D f and section 6.2 for the origin of the 1 or 3 "exceptional" zero modes in this case), we readily arrive at the following summary: which fully agrees with numerical results in [54]. This input should be plugged into 20) where N + f −3 has been obtained by (6.17). This completes our discussion of even N c . For odd N c , (6.17) and (6.18) still hold in the range 0 ≤ f < (N c − 3)/2 (see table 8). For f near N c /2 we only have to substitute (6.19) into The numerical results in appendix B agree with the formulas derived in this subsection.

Generalization toq > 3
We now generalize the preceding classification scheme to the N = 2 SYK model with H = {Q, Q} andq complex fermions in the supercharge, whereq is odd, i.e., This is a counterpart of (4.1) with N = 1. Forq = 3 it reverts to (6.1). The tables (6.2) and (6.13) forq = 3 are now generalized to We also analyzed the dimensions N ±,z f of the subspaces, for which formulas similar to those in section 6.4 can be derived. Forq = 5, we have numerically confirmed up to N c = 17 that all exact zero modes of H reside in V f with |f − N c /2| ≤ 5/2. The last inequality is saturated only for N c = 7 and 11 by just 2 zero modes in each case. This is not only consistent with our heuristic argument in section 6.4 but also conforms to the claim at large N c [121,122] that all zero modes should satisfy |f − N c /2| <q/2. In the regime N c 1 one can ignore O(1) exceptional zero modes and the strict inequality may be justified.

Conclusions
In this paper we have completed the symmetry classification of SYK models with N = 0, 1, and 2 SUSY on the basis of the Altland-Zirnbauer theory of random matrices (table 1). The symmetry classes of RMT not only tell us the level degeneracies of the Hamiltonian but also offer a diagnostic tool of quantum chaos through level correlations in the bulk of the spectrum. Furthermore, when the spectral mirror symmetry is present, RMT precisely predicts universal level correlation functions in the vicinity of the origin (also known as hard edge or microscopic domain [79]). The present work can be viewed as a generalization of preceding works [66][67][68][69]72] that analyzed the level statistics of the N = 0 and 1 SYK models solely with a 4-body interaction. 15 Our new results include the following: 1. The symmetry classification of the N = 0 SYK model was given for a Hamiltonian with the most generic q-body interaction. The result, summarized in tables 2 and 3, includes the RMT classes C and D that did not show up in the preceding classification of [66][67][68][69]72]. Our results were corroborated by detailed numerics (figure 1).
2. We numerically compared the smallest eigenvalue distributions in the N = 0 SYK model with q = 6 with the RMT predictions of class C and D, finding excellent agreement ( figure 2).
3. The symmetry classification of the N = 1 SYK model was given for a supercharge with the most generic interaction ofq Majorana fermions (tables 4 and 5). This extends [72] which investigated onlyq = 3. Our results were corroborated by detailed numerics ( figure 3). 4. We numerically compared the smallest eigenvalue distributions in the N = 1 SYK model withq = 3 and 5 with the RMT predictions, finding excellent agreement (figures 4 and 5). This confirms the hard-edge universality of the N = 1 SYK model for the first time and is relevant for the thermodynamics of this model at low temperatures comparable to the energy scale of the SUSY breaking.
5. We proposed an intriguing new SYK-type model which lacks SUSY but whose Hamiltonian is semi-positive definite and has an extensive number of zero-energy states (section 5). The symmetry classification based on RMT was provided, and a detailed numerical analysis of the spectra both in the bulk and near the origin was performed, resulting in agreement with the RMT predictions.
6. We completed the RMT classification of the N = 2 SYK model for the first time. This model is qualitatively different from its N = 0 and 1 cousins in various aspects. It is a model of complex fermions rather than Majorana fermions, and it has a U(1) Rsymmetry. The symmetry classification of this model is nontrivial because the structure of its Hilbert space is far more complex (see figure 11 for an example) than that of the N = 0 SYK model with complex fermions considered previously in [29,35,44,66,115].
7. In section 6.2 we succeeded in giving a logical explanation for the curious fact [54,57] that, in the N = 2 SYK model, the number of zero-energy ground states exactly agrees with the lower bound from the Witten index in some cases but not in other cases. In short, this is due to the dichotomy between the odd dimensionality of the Hilbert space and Kramers degeneracy.
This work can be extended in several directions. First, our analysis of spectral properties of the Hamiltonian could be further deepened by using probes that are sensitive to long-range correlations of energy levels, like the level number variance Σ 2 (L) and the spectral rigidity ∆ 3 (L) [2,80]. Investigating the spectral form factor of the N = 2 SYK model and making a quantitative comparison with RMT along the lines of [68] is another future direction, although physical interpretation of the ramp, dip, etc., of the spectral form factor as a signature of quantum chaos is rather subtle [45]. Finally, we note that there is no analytical result for the global spectral density of the N = 1 and 2 SYK models, although an accurate formula is already known for the N = 0 model [67][68][69]. We wish to address some of these problems in the future.