On thermalization in the SYK and supersymmetric SYK models

The eigenstate thermalization hypothesis is a compelling conjecture which strives to explain the apparent thermal behavior of generic observables in closed quantum systems. Although we are far from a complete analytic understanding, quantum chaos is often seen as a strong indication that the ansatz holds true. In this paper, we address the thermalization of energy eigenstates in the Sachdev-Ye-Kitaev model, a maximally chaotic model of strongly-interacting Majorana fermions. We numerically investigate eigenstate thermalization for specific few-body operators in the original SYK model as well as its $\mathcal{N}=1$ supersymmetric extension and find evidence that these models satisfy ETH. We discuss the implications of ETH for a gravitational dual and the quantum information-theoretic properties of SYK it suggests.


Introduction
A deep understanding of the nature of thermalization and an approach to equilibrium in closed quantum systems remains to be seen. For a large class of interacting quantum many-body systems, the eigenstate thermalization hypothesis (ETH) [1,2] is an important conjecture which stands to help explain such a mechanism. The claim is that, in some large but isolated quantum mechanical system, the energy eigenstates look thermal, meaning that operators corresponding to physical observables have diagonal matrix elements given by microcanonical ensemble and off-diagonal elements suppressed by the entropy of the system. This conjecture gives us insight into the evolution of out-of-equilibrium states and is further supported numerically for few-body operators in specific strongly-correlated lattice models (see [3] as well as [4] and references therein). The coincidence and difference between eigenstate expectation values and those in microcanonical ensemble provide a quantitative interpretation of how pure states become thermalized in a quantum chaotic system.
The mechanism of thermalization in a quantum theory is closely related to the statistical mechanics of black holes via the AdS/CFT correspondence, where spacetime geometries are dual to some specific quantum states in the boundary field theories. Understanding thermalization holographically [5] can help illustrate the types of bulk geometries dual to thermal states in the CFT and, more ambitiously, could shed light on black hole formation and evaporation. Moreover, the prototypical examples of eigenstate thermalization are quantum chaotic systems [2,6]. As black holes are both fast scramblers [7,8] and maximally chaotic [9], explaining the emergence of thermalization from chaotic evolution will help us understand how black holes process quantum information [10][11][12].
The Sachdev-Ye-Kitaev (SYK) model [13,14] is a concrete arena to address some of these questions. SYK is a quantum mechanical model of N strongly-interacting Majorana fermions with all-to-all random couplings. The model is solvable at large N and possesses many intriguing features indicative of a gravitational dual, such as an emergent conformal symmetry at low energies and a zero temperature entropy [13,15]. The theory is also maximally chaotic in the sense that an explicit calculation of the out-of-time ordered correlation function [13,15,16] shows a Lyapunov growth which saturates the chaos bound [17], a seemingly unique feature of gravitational theories [18] as well as CFTs with a gravitational dual [19]. The low energy description of the theory is given by a Schwarzian action which also captures dilaton gravity in AdS 2 [20,21]. Finally, SYK exhibits both random matrix universality in the spectrum [22,23] as well as random matrix behavior in the late-time dynamics of the model [24]. Many of these features are present in generalizations of SYK that have been discussed, including a supersymmetric extension of the model.
In this paper, we will discuss ETH in SYK and SYK-like models. We will check the original SYK model and the N = 1 supersymmetric generalization of SYK given by Fu, Gaiotto, Maldacena, and Sachdev [25]. By exact diagonalization, we numerically compute the matrix element in the energy eigenstates and verify that these models satisfy ETH. Our numerics quantitatively illustrate that pure states in these holographic models can thermalize.
There has already been some discussion in the literature of thermalization in SYK, both looking at ETH as well as considering more dynamical questions. It has been shown [26] that ETH is satisfied in the free fermion analog of SYK. [27] discussed a quench setup in the SYK model and analytically showed that the system rapidly thermalizes in a certain limit. Diffusion and spread of entanglement have also been discussed in a generalization of SYK with spatial locality [28,29]. Recently, [30] showed that specific pure states in the SYK model are thermal in the large N limit. Most relevant to this work, [31] investigated ETH in the complex fermion version of SYK model and found evidence for eigenstate thermalization. While the model with complex fermion version of SYK [32][33][34] shares many of the same features, albeit with some subtle differences [35], one is still inclined to believe that thermalization and ETH will hold much in the same way for both Majorana SYK and its supersymmetric extension. We investigate ETH explicitly in these models and, for some few-body operators, check the ansatz for both diagonal and off-diagonal terms, and investigate the Gaussianity of fluctuations. This paper is organized as follows: In Section 2 we review the SYK model and its supersymmetric generalization, then discuss the basics of the eigenstate thermalization hypothesis. In Section 3 we present a numerical investigation, by exact diagonalization, of the matrix elements of specific operators evaluated in pure states. We read off the corresponding behavior of the diagonal and off-diagonal elements. We conclude in Section 4 and discuss the implications for thermal properties of the bulk and quantum information aspects of SYK. In Appendices A and B we mention two technical points, including a constraint from particlehole symmetry in the SYK model, and the issue of unitary averaging of distributions when testing Gaussianities of fluctuations.

SYK and supersymmetric SYK
Before we delve into a numerical investigation, we first review the quantum mechanical models we wish to consider. The Hamiltonian of the (q-point) SYK model [13,15] is given by The system has all-to-all q-body interactions with random couplings J i 1 i 2 ...iq with mean and variance i.e. we disorder average over independent Gaussian random variables J i 1 i 2 ...iq , and J SYK is a positive constant. The model is exactly solvable in the limit 1 βJ SYK N , where an emergent conformal symmetry can be used to compute the correlation functions of the theory. A calculation of the out-of-time-ordered four point function of the theory [13,15] shows that the Lyaponov exponent, governing the time-dependent growth of 1/N corrections, saturates the chaos bound λ = 2π/β at large N . The conformal symmetry is spontaneously and explicitly broken and one can obtain a low-energy description of the model given by the Schwarzian, capturing the dynamics of the pseudo-Goldstone mode. The Schwarzian theory can also be understood as the low-energy dilaton gravity description of the two-dimensional bulk dual of [21,36].
Interesting generalizations to supersymmetric versions of the SYK model were originally given in [25], which have been further studied in [37][38][39][40][41][42]. Here we just discuss the N = 1 generalization. In the (2q − 2)-point N = 1 supersymmetric model, one can construct an N = 1 supercharge Q as a q-body Majorana interaction. The Hamiltonian is then given by the square of the supercharge and again we have random couplings C i 1 i 2 ...iq with mean and variance Just as above, C i 1 i 2 ...iq is an antisymmetric Gaussian random tensor and J N =1 is a constant. The bosonic operators b i of the theory, supersymmetric partners of ψ i , are given by the application of the supercharge as b i = Qψ i . This model shares many of the intriguing large N features observed in the Majorana model, such as chaotic correlation functions, zerotemperature entropy, and an emergent superconformal symmetry that is spontaneously and explicitly broken, giving rise to a Schwarzian-like low-energy action [25]. In this paper we will only consider the simplest non-trivial interactions in the SYK models, meaning we set q = 4 in the SYK model and consider q = 3 supercharges in the N = 1 supersymmetric model. We will also choose J SYK = J N =1 = 1 without loss of generality, as the coupling can simply be regarded as a rescaling of energy.
As our best understood examples of AdS/CFT involve string theory in AdS, with a supergravity description at low-energies, the supersymmetric extension of SYK might be helpful in constructing a more concrete gravity dual in two dimension AdS or understanding how the dilaton gravity might be UV completed. With this in mind, we should note the similarities and differences for thermalization in the original model compared to its supersymmetric extension.

Eigenstate thermalization
ETH is a statement about the thermal nature of the energy eigenstates in a closed quantum system [1,2], where the thermal behavior emerges from pure states for some generic local operators. The claim is that the matrix elements of an operator O, evaluated in the basis of energy eigenstates, should satisfy the following: Moreover, the function f O (Ē, ω) is a smooth function ofĒ and ω, S(Ē) is the entropy of the system, and R mn is a random variable with zero mean and unit variance. ETH states that in a system with many degrees of freedom, the diagonal matrix elements dominate, while the off-diagonal terms are suppressed by the entropy of the system. The distribution of R mn could be real or complex, depending on the symmetry of the system we consider and the basis we choose. In either case complex conjugate constrains the matrix elements of R mn and the function f O . ETH, as formulated above, can help explain the emergence of thermal expectation values of certain operators in the time evolution of quantum many-body systems. Following the discussion in [4,6], we prepare a pure state |ψ in an isolated quantum system with Hamiltonian H and energy eigenstates |n , where the time evolved state is given as with coefficients c n given by c n = n|ψ i . The expectation value of any operator O in the time-evolved state is where we denote matrix elements of the operator O mn = m|O|n . At late times, the offdiagonal terms are weighted with a dephasing exponential factor, which causes them to decay. More precisely, we can take an infinite time average and find assuming that there are no, or at least a nonextensive number of, degeneracies (i.e. not scaling with the system size). This late-time value of the operator is often called its expectation value in the diagonal ensemble. One should note that the leading term in this expectation value at late time is constant and only depends on the construction of the initial state c n . On the other hand, the ensemble average in the microcanonical ensemble is where we sum over the states in some narrow window of energies near E, i.e. sum over n such that |E − E n | < ∆E, and N ∆E denotes the number of states in that window. Now assuming that the initial state is prepared in a narrow energy window around E, according to the assumption Eq. (2.5),Ō(E) is a smooth continuous function and thus near the energy window E the matrix element O nn ≈Ō(E). Thus we have which shows the emergence of thermalization in a system that satisfies ETH. All we needed to assume was that the diagonal and microcanonical distributions were sufficiently narrow. More precisely, Srednicki gave a criteria [43] for the smallness of ∆E relative to the smooth function of diagonal elements, which we identify with microcanonical average of the operator. If the energy fluctuations are subextensive, then we conclude that the microcanonical and long-time averages agree, and thus the system thermalizes. Clearly, ETH cannot hold for all possible operators. In fact, one can construct operators to trivially violate ETH in all possible quantum systems, i.e. projection operators or operators constructed from the Hamiltonian. However, ETH is still far from trivial; these artificial operators are nonlocal and thought not to correspond to real physical observables. It is often claimed that ETH should hold for few-body operators that are sufficiently local. Interestingly [44], it was argued that ETH should hold for operators with spatial support up to half of the system size. ETH was first observed numerically in a bosonic lattice model [3], and further been observed in a variety of strongly-interacting lattice models, including spin chains, fermionic systems, etc. [45][46][47][48][49], including checks in every eigenstate [50].
Another subtlety arises when discussing ETH in a system with disorder. In a disordered system the eigenvalues and eigenstates are not fixed but instead depend on the specific realization of disordered couplings. But ETH still makes sense in the following way. For a given energy E and some number of disorder realizations, the energy eigenvalues of those realizations will vary around E. For a given energy window ∆E, we can look at the energy eigenstates for all possible eigenvalues that are in the range [E − ∆E, E + ∆E] and compute the statistics of those eigenstates (for instance, for fixedĒ and ω, or equivalently, E m and E n , one considers the eigenvalues in the range [E m − ∆E, E m − ∆E] and [E n − ∆E, E n − ∆E] for a given number of random realizations). For ETH in disordered systems, the eigenstates |m and |n are not fixed but rather ensembles of eigenstates in an energy window ∆E. Similar issues are discussed and addressed in [4,51].
There is a long history between ETH and systems which exhibit quantum chaos 1 (as nicely reviewed in [4]). If a chaotic system is highly random it sometimes suffices to consider a Hamiltonian given by a random matrix. In this case, without assuming the initial state is narrow, we can often directly takeŌ(E) outside of the summation. Given the maximal chaos and late-time random matrix behavior in SYK-like models, we expect that they should satisfy ETH. For a discussion of the random matrix behavior of SYK, see [12,24]; for supersymmetric models, see [40,41] and [52]. Recent work also provides evidence which supports this conjecture; for instance, [30] shows that aspecific low-energy pure state correlators are equal to thermal correlators at large N , while [31] gave numerical evidence that the complex SYK model satisfies ETH. In the following section, we will show precisely how ETH is satisfied for the Majorana SYK model as well as its supersymmetric extension.

Setup
In this part we will describe our numerical results. We study the system by exact diagonalization. For SYK model and supersymmetric SYK model, we use the Clifford algebra representation of Majorana fermions, where we employ the standard Pauli matrices as To test ETH, we must restrict to a specific class of operators. One can consider the particle number operator, where k ∈ {1, 2, . . . , N d }, and c k is a complex fermion. We may write Majoranas as The particle number operator n k could also be written as where S k is defined as [30] S k = 2iψ 2k−1 ψ 2k .
where, in a slight abuse of notation, the subscript denotes s ≡ s = {s 1 , . . . , s k }. The 2 N d states |B s form a complete basis of the Hilbert space (and in fact, are unit vectors in the representation Eq. (3.1)). Thus, we may write where half of the λ's in the sum are +1 and other half are −1. As the operator S k has no thermal expectation value in the diagonal terms and we only have random fluctuations in both the diagonal and off-diagonal elements. Note that the diagonal terms for n k are nontrivial and fluctuate, except for when the particle number is fixed by symmetry as we discuss in Appendix A.
To further test ETH, another natural 2-body operator to consider is the hopping operator Note that for k 1 = k 2 , this operator has vanishing expectation value in the microcanonical ensemble Tr(ρ mc h k 1 k 2 ) = 0. So agreement with ETH will mean the diagonal elements are close to zero. In the following sections, we both check ETH for n k and h k 1 ,k 2 . For simplicity, we only consider models with even N d in performing our numerics. In both the SYK model and its supersymmetric extension, fermion number parity is conserved, and thus we can focus on one single parity sector to study the ETH. Without loss of generality, we will simply look at the even parity sector.
We quickly note mention some aspects of the expectation values of particle number operators specific to the Majorana SYK and its supersymmetric generalization. These models have a particle-hole symmetry, where the antiuntiary operator P commutes with the Hamiltonian H [22,24]. P maps a state with fermion number Q to N d − Q, or equivalently sends the filling from ν to 1 − ν when acting on a pure state, where ν = n f /N d . If N mod 8 = 0 where P 2 = 1, then P maps each energy eigenstate to itself, which means we have the constant filling ν = 1/2 for all possible energy eigenstates. In other cases with even N d we are considering, the particle-hole symmetry is non-trivial and enforces a degeneracy in each parity sector. For a more details on particle-hole symmetry and the constraint described above, see Appendix A. Randomly looking at energy eigenstates states from enlarged eigenspace, the particle number expectations should be fluctuating around 1/2. Thus, when calculating diagonal elements we will avoid the N mod 8 = 0 case.

Density of states
Before considering the matrix elements of specific operators projecting on the spectrum and discussing the thermalization of the eigenstate, in Figure 1 we plot the density of states for SYK and N = 1 supersymmetric SYK model for completeness (those results and related numerics are previously obtained in, for instance, [15,40]). The spectral density of SYK has been further studied in [23,24,53]. The density of states are obtained for multiple realization of randomness in the corresponding Hamiltonian. In the (relatively) large N , the distribution of eigenvalues for SYK model will be more closer to a combination between Gaussian distribution (by the central limit theorem) and the Wigner's semicircle, while in the supersymmetric model, the supercharge is Gaussian-like, thus the distribution of energy eigenvalues will approach the square of Gaussian randomness: the Marchenko-Pastur distribution in a large N limit of the one point distribution in the Wishart-Laguerre ensemble (see [40] for more details).

Direct checks of ETH
Now we come to direct checks of ETH in SYK models. First, we make density plots of the matrix elements of number operators and hopping operators projecting onto the energy eigenstates in both models. The evidence of eigenstate thermalization is apparent even in a single random realization of disorder, where we do not need to take a thermal average. For the particle number operators, as shown in Figure 2 and Figure 3 for SYK and supersymmetric model respectively, the density plot clearly shows the diagonal dominance in the matrix elements O mn , while the off-diagonal terms are suppressed by the entropy of the system, with small random fluctuations. However, if we consider hopping operators, as displayed in Figure 4 and Figure 5, there is noŌ(Ē) contribution in the microcanonical ensemble thus confirmation with ETH is simply the appearance of small random fluctuations in the matrix elements, as observed.
Another direct investigation of ETH is to evaluate the difference between the eigenstate projection and the expectation value from the microcanonical ensemble [54], i.e. the difference Figure 3. Eigenstate thermalization density plots for the particle number operator in the supersymmetric SYK model. The setup is the same as described in Figure 2.
between diagonal and microcanonical ensembles. To be concrete, we consider where we sum over all energies (i.e. eigenvalue labels n) in the energy window [E − ∆E, E + ∆E], andŌ(E) is given as the microcanonical ensemble average for the energy E. Specifically, we consider the particle number operator O = n k . As explained previously, the diagonal terms are strictly 1/2 for N mod 8 = 0 for each realization of disorder parameter, thus we specifically consider N = 12 and N = 20 in Figure 6 for SYK and its supersymmetric extension. The results indicate thermalization of eigenstates as the relative difference is very small and becomes further suppressed as we increase the number of fermions. The difference between the microcanonical and diagonal predictions depends on the energy (or the effective temperature in those models). But in the parameter range we are testing, we conclude that the predictions from the microcanonical and diagonal ensembles are very close, although there are still random fluctuations due to the finiteness of the sample size.

Diagonal terms
Some checks of eigenstate thermalization can also made in the diagonal terms O nn . For the single fermion particle number n k , we should observe fluctuations around the microcanonical Figure 4. Eigenstate thermalization density plots for single hopping operators in the original SYK model, with the same setup as described in Figure 2. Note that the microcanonical value of the hopping operator is zero, so agreement with ETH is simply the appearance of the off-diagonal Gaussian fluctuations.
value 1/2 (except for N mod 8 = 0 where there are no fluctuations as the particle number is fixed by symmetry). In the other cases, the system may still have a degeneracy due to the internal symmetry in each parity sector and thus we only choose one eigenstate in each energy level.
One can see that for larger N , and with substantially more eigenvalues, the average deviation from 1/2 should be a function of the eigenvalues and need not be a constant function. Thus we plot the smoothed deviations from the microcanonical value 1/2, for N = 12 and N = 20 with multiple realizations of disorder, in Figure 8. One can see that at lower N (the N = 12 case), the average deviation between diagonal terms and 1/2 is roughly a constant, while for larger N (N = 20 here) there exists stronger fluctuations at larger energies, and more so at even larger N , for instance even in single realization in Figure 7. This fact confirms our observation in Figure 6.

Off-diagonal terms
Another crucial numerical check of ETH is to consider the off-diagonal elements, looking at the statistics of the random fluctuations. We consider the single particle number and hopping operators in the SYK and supersymmetric models and in Figure 9, plot the statistical average of the off-diagonal elements as a function of the energy difference (ω = E m − E n ). One can observe some generic features in the numerics: First, the off-diagonal elements have larger variations in the middle of the curve, where the energy difference is small. Second, the difference between the bulk and the edge in the ω − |O| distribution becomes sharper when N becomes larger. Finally, in the supersymmetric model, the slope of the curve is smaller than the original SYK model. Some of these features are also seen in the complex fermion model [31]. We also note that the decreasing variance of the fluctuations at larger energy differences (large ω) is a universal feature in chaotic systems satisfying ETH, and the averaged decay of fluctuations at larger ω sets the scale of the Thouless energy and the diffusion time.

Gaussianity of fluctuations
The fluctuation term R mn in the off-diagonal contribution is often assumed to be Gaussian (see for instance, [2]), which has been checked numerically in disordered models [51]. Here we will verify that the Gaussian condition is satisfied for both SYK and supersymmetric SYK models in each parity sector. To start, one can directly compute the distribution of off-diagonal terms from multiple realizations of disorder parameters for fixed energy averageĒ and difference ω appearing in the spectrum. Generically, the distribution of R mn should be complex, and an overall complex phases do not change the distribution of R mn (see Appendix B a discussion of this). However, in the case of even N d we are considering, the only exception is N mod 8 = 0 Figure 6. The relative difference between the microcanonical ensemble predictions and the diagonal ensemble predictions in the SYK model (above) and the supersymmetric SYK model (below) as a function of energy E/N . We choose multiple realizations for N = 12 (left, 20000 realizations) and N = 20 (right, 1000 realizations) models respectively. We consider the energy window ∆E = 0.01 in both cases. In this plot, we totally collect 1001 data points and in order to smooth the random fluctuation we use the moving average for 30 nearest energy data points. The remaining random fluctuation may due to the finiteness of the number of random generalizations. case, which corresponds to real representation in random matrix theory classification (see [22,24,40] for reference), for both SYK and supersymmetric SYK models. In this case, one can use the real representation of Clifford algebra to ensure reality of the eigenvectors.
In Figure 10 we check the Gaussian distribution of R mn for fixed energies. Those numerical results clearly show the Gaussianity of these fluctuations, while the numerical non-Gaussianity may due to the finiteness of sample numbers, the finiteness of N , and also the numerical artifacts when smoothing the distribution curves.

Pure state thermalization
There is another related statement of thermalization, which has a clear geometric interpretation in the SYK model, regarding the thermalization of the correlators under the Euclidean evolution of a certain set of pure state [30]. The statement is given as the following: consider an arbitrary eigenstate of S k , namely one of the |B s states we discussed previously, and define a low energy projection with the following Euclidean evolution   with β = 2 . Moreover, a single pure state B is also very useful to simulate the thermal state in the SYK model, namely the correlators are thermalized. Inserting the identity into the above expression, each term in the sum is equal in the large N limit [30], meaning Moreover, we can define where . (3.14) These expressions are valid in the large N limit βJ N . In the limit τ → 0, each term in the sum of the thermal correlator gives ψ i (0)ψ i (0) β /Z = 1/2. The thermalization of these pure state is a consequence of symmetry of the replica symmetric disordered model (a subgroup of O(N ) termed the flip group). Considering that N = 1 supersymmetry only changes the distribution of the Hamiltonian (up to a constant energy shift), the flip group argument is still valid. As a result, we expect that the correlators in the low energy pure are still thermal.
We confirmed these expectations with numerics, comparing the pure state inner product and the partition function (namely, verifying Eq. (3.12)) and comparing the diagonal and off-diagonal correlators with thermal ones (namely, verifying Eq. (3.13)), and find very good numerical agreement. As an example, the comparison with the thermal correlators is shown in Figure 11. This indicates thermalization of low energy projected pure states.

Discussion and Conclusion
In this paper we numerically investigated the thermalization of energy eigenstates in the Majorana SYK and supersymmetric SYK models. We focused primarily on conventional few-body operators, the single particle number and the hopping operators, to verify that these models satisfy ETH. Our results include the direct verification of ETH by checking the diagonal dominance and the relative difference between microcanonical ensembles, the Figure 11. Pure state thermalization: we consider N = 12 models for both SYK and the supersymmetric model and compute the relative difference between the correlation function in the pure state (LHS of Eq. 3.13) and in the thermal state (RHS of Eq. 3.13). We fix β = 2 and consider 500 realizations. The small relative difference supports the pure state thermalization in both models. statistics of the diagonal terms and the off-diagonal terms, and verifying the distribution of off-diagonal fluctuations. The numerical evidence indicates that SYK and its N = 1 supersymmetric extensions satisfy ETH as defined in Eq. (2.5).
Although ETH applies in a broad class of quantum systems, the statement that ETH holds in the SYK model confirms a basic expectation of a system with a gravitational dual: correlation functions in energy eigenstates of look thermal, as anticipated from local bulk effective field theory and in line with results from 2d CFT [55]. There is a direct relation [31] between the thermalization of eigenstates and thermal correlators in low-energy pure states constructed in [30]. In SYK we checked this thermalization for the pure state correlators in the large N limit. As these pure states are suggestive of two-dimensional black hole microstates, operators related to these pure states could be used to study transversable wormholes [56,57].
Beyond making statements about black hole thermalization in a putative gravitational dual, the claim that SYK satisfies ETH is interesting from a quantum information perspective. The frame potential for an ensemble of unitaries E is a measure of how close that ensemble is to forming a unitary k-design [58], i.e. reproducing the first k moments of the Haar ensemble. Building on this, [12] defined the notion of k-invariance as a measure of how scrambled a system is under chaotic time-evolution, giving rise to a random matrix description. If a system achieves k-invariance, the eigenstates may be treated as random vectors. For a system satisfying ETH, the late-time dynamics should be k-invariant and, in this precise sense, random. Given recent information-theoretic approaches to studying chaos and complexity [10][11][12], it would be nice to understand the precise role ETH plays in this framework and whether it is a necessary condition for information scrambling or chaotic decay of correlation function. Moreover, evidence that SYK satisfies ETH hints at some error-correcting properties. As the condition for correctability of an approximate quantum error correcting code is satisfied by finite energy density eigenstates of systems which obey ETH [59], an analytic understanding of ETH in the SYK model might elucidate its underlying code properties.
Generically, one expects ETH to hold for few-body physical operators; but the statement might be more subtle for non-local systems. It would also be interesting to investigate eigenstate thermalization in generalizations of SYK with spatial locality [28] and in tensor models [60][61][62][63], to understand the role all-to-all interactions and quenched disorder play in the thermal structure of eigenstates. There are more fundamental subtleties about which energy windows and for which operators ETH actually applies. Recent progress has been made in addressing these issues, both in the context of chaotic many-body systems [44,64] as well as in 2d CFTs [65,66], alongside work towards generalizing ETH [67,68]. Lastly, as ETH has not been proven to hold in any generic class of interacting models and analytic results of systems that satisfy ETH are few and far between, one might hope that SYK, a rare example of a strongly-coupled but solvable quantum mechanical model, affords the opportunity to prove ETH analytically. We leave these topics to future research. 2644). JL is supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award Number DE-SC0011632. YZ is supported by the graduate student program at the Perimeter Institute.

A Constraints from Particle-Hole Symmetry
Both the SYK model and its supersymmetric extension admit a particle-hole symmetry [22,24] and where K is the antiunitary complex conjugation operator. The particle-hole operator squares to P 2 = (−1) N/4 and acts on the fermions as We can consider the action of P on an energy eigenstate |n . If N/2 is odd the the state is mapped from one charge parity sector to the other, exchanging even and odd. If N/2 is even P maps each sector to itself, but if N mod 8 = 4 then there is a degeneracy in each sector as P 2 = −1 implies P |n must take us to a different eigenstate. If N mod 8 = 0, there is no degeneracy and P maps energy eigenstates to themselves (up to a phase). This implies that when N mod 8 = 0, we have an exact half-filling for energy eigenstates, namely that for any eigenstate |n and any k, the number operator n k = c † k c k has n|c † k c k |n = where we commute P through as P c † k c k = c k c † k P . When N mod 8 = 0, P then takes |n to itself and we find the desired result.
Note that the above argument extends to any operator O 2k of the following form O 2k = ψ i 1 ψ i 2 · · · ψ i 2k , then n|P O 2k P |n = (−1) k n|O 2k |n , (A.5) for operators with no coincident indices i a = i b . Since we can assume that i a = a, O 2k is proportional to a product of 2c † α c α − 1 with α running from 1 to k. In other words, for odd k, O 2k always has zero trace on any energy eigenstates (degenerate or not), in particular, the diagonal terms of O 2k vanish when N mod 8 = 0.

B Unitary Averaging of the Distribution
Considering the ETH ansatz for the SYK model (and its supersymmetric generalizations), m|O|n =Ō(Ē)δ mn + e −S(Ē)/2 f O (Ē, ω)R mn , (B.1) a natural question then arises: if we want to look at the statistics of the off-diagonal random variables R mn , how should we make sense of eigenstates |m and |n when doing numerics by exact diagonalization? More specifically, as we are averaging over random Hamiltonians, each realization comes with a choice of eigenstate up to a nonphysical phase. Here we shall propose a ETH for the SYK model and its supersymmetric generalizations by stating that non-physical phases should be averaged by a uniform distribution to remove all phase fluctuations. To be precise, m, θ m |O|n, θ n θ =Ō(Ē)δ mn + e −S(Ē)/2 f O (Ē, ω)R mn |k, θ k = e iθ k |k , k ∈ {1, 2, · · · , 2 N d } , where θ k are independently chosen from 0 to 2π with uniform distribution, and subscript θ indicates an average over all θ k . We also assume that the choice of |n are not random, i.e. for a set of coupling constants there is only one set of eigenstates with fixed phases. In the numerical realization using Mathematica, the assumption holds because eigenstates are generated with some definite algebraic procedure.
Here we remark that the the distribution of phase of R mn does not depend on the choice of the original (fixed) choice of the phase of |n , because different choice of phase correspond to a rotation in complex plane and averages out over the uniform distribution of θ k . In Figure  12 we check the Gaussian distribution of R mn for fixed energies. Compared to the unmodified distribution in Figure 10, the non-Gaussianity is suppressed.