Eigenstate thermalization in the Sachdev-Ye-Kitaev model

The eigenstate thermalization hypothesis (ETH) explains how closed unitary quantum systems can exhibit thermal behavior in pure states. In this work we examine a recently proposed microscopic model of a black hole in AdS$_2$, the so-called Sachdev-Ye-Kitaev (SYK) model. We show that this model satisfies the eigenstate thermalization hypothesis by solving the system in exact diagonalization. Using these results we also study the behavior, in eigenstates, of various measures of thermalization and scrambling of information. We establish that two-point functions in finite-energy eigenstates approximate closely their thermal counterparts and that information is scrambled in individual eigenstates. We study both the eigenstates of a single random realization of the model, as well as the model obtained after averaging of the random disordered couplings. We use our results to comment on the implications for thermal states of the dual theory, i.e. the AdS$_2$ black hole.


Introduction
Via holographic duality (or 'AdS/CFT') black holes are described by thermal states of a dual quantum field theory. The process of black hole formation and evaporation corresponds to the process of thermalization of certain (unitary) quantum field theories, evolving from nonequilibrium initial states towards thermal equilibrium. However, many important questions pertaining to this attractively simple picture remain mysterious, as expressed notably in the information loss paradox [1], and its more recent ramifications [2]. An important aspect of a satisfactory resolution of these puzzles within holography is a precise understanding of the relation between information loss, unitary evolution and thermalization of the boundary theory. Put succinctly the question of how thermal [3,4] a generic high-energy eigenstate looks in a theory with a holographic dual is closely tied to the question of how legitimately one may consider the dual geometry of this pure state to be a black hole. It is therefore important to understand the detailed mechanism of thermalization in any given dual field theory model, as well as to extract lessons for the general class of theories with bulk duals, wherever these are available. As is often the case, the most promising starting points are the simplest instances which still contain enough structure to allow one to address the subtleties also present in the higher-dimensional, more complicated theories of primary concern. For this reason, lower-dimensional holographic models, such as AdS 3 /CFT 2 (i.e. three-dimensional gravity) [5][6][7][8], as well as the even lower-dimensional case of NAdS 2 /NCFT 2 ('Near AdS 2 ', 'Near CFT 2 ') [9][10][11][12] have once again risen to the fore.
The subject of the present study is a recently proposed many-body quantum mechanical system of interacting fermions [9,10], closely related to an older one with a controllable spinglass ground state [13], the so-called Sachev-Ye-Kitaev (SYK) model. This is a quenched disorder quantum system with random all-to-all couplings, which shows emergent approximate conformal invariance in the infrared as well as an extensive entropy at low temperature [9][10][11]14]. Furthermore, it has also been demonstrated [9] that the model exhibits maximal quantum chaos (in the sense of [15]), diagnosed by a certain Lyapunov exponent (λ L ) extracted from thermal four-point correlation functions [9,11,14] (see also [16,17] for interesting aspects of these correlation functions). In the SYK model then, the question of thermalization is one about the behavior of a simple 1 many-body Hamiltonian under unitary evolution, i.e. the study of thermalization in closed many-body quantum systems. Much has been learned about such situations (and we refer the reader to the review [22] for more details and references). The simplicity of the SYK model means that it is a perfect candidate for a detailed study of the mechanism of thermalization within the setting of holographic duality with access to the deeply quantum regime.
It should be kept in mind that owing to the ensemble average implicit in the model, there are potential subtleties in connecting results in the SYK model to black-hole physics. To 1 Simple enough that various experimental approaches have been put forward [18][19][20][21]. start addressing this issue, [23][24][25] have shown that certain tensor theories defined without averaging over an ensemble, give rise to the same large-N limit as the SYK models. Here we also study the question of thermalization in a given random realization of the SYK model, for which the previous cautionary remark does not apply. To the extent that one limits oneself to self-averaging quantities the same then holds for the disorder averaged model.
In this paper we demonstrate in detail that the SYK model satisfies the properties of the eigenstate thermalization hypothesis and therefore establish this scenario as the appropriate mechanism for thermalization in this prototypical example of holographic duality. We expect that other versions of holographic duality also exhibit eigenstate thermalization, and we offer some further comments in the discussion section. A discussion on ETH and holography may be found in [26]. We are here studying a version of the SYK model in which four fermions interact, but one can extend this to a q−fermion interaction [11]. This was used by [27] in order to analytically study thermalization following a coupling quench. Our study elucidates the thermal structure inherent in pure states more generally, and it would be enlightening to derive some of our numerical results analytically, perhaps starting at large q. The free model q = 2 was already shown to satisfy eigenstate thermalization in [28].

Summary of results
Let us briefly summarize the main results of this paper. First and foremost we establish, by exactly diagonalizing the complex SYK Hamiltonian for up to 17 sites, that expectation values of non-extensive -that is those involving a few sites only -operators are to a very good approximation thermal. In particular their matrix elements take on the expected form encapsulated in the eigenstate thermalization hypothesis. This has interesting ramifications for the holographic dual, and we address some of these in the discussion section.
By studying off-diagonal matrix elements of non-extensive operators we establish that the SYK model behaves like a random-matrix theory (RMT) for a certain range of energies, but more generally deviates from such RMT behavior. By analogy with the theory of disordered conductors we refer to the energy scale at which deviations from RMT are observed as the Thouless energy E T . We find that the Thouless energy is controlled by the coupling as E T ∼ J 2 . The more strongly coupled the system, the larger the energy range for which it exhibits RMT behavior.
Having established that eigenstates behave thermally we compare correlation functions of non-extensive operators in eigenstates with their corresponding thermal expectation values.
We find that two-point and four-point correlation functions in eigenstates are qualitatively similar to thermal averages, but do differ in their detailed structure. We also study correlations in large superpositions of eigenstates which approximate thermal averages with respect to a canonical density matrix. To the extent that we can extrapolate these results to large values of N this suggests that individual eigenstates can in some sense be considered to be dual to the black-hole geometry in the putative dual, although we make no statement about its interior (see discussion section).
This motivates us to consider measures of scrambling in eigenstates, which we find to behave in accordance with expectations from combining known results in the canonical ensemble with our results on eigenstate thermalization. We therefore expect that there is a large−N eigenstate equivalent of the maximal scrambling exponent satisfied by the SYK model, as detailed in (2.11) below.

Background
In this section we introduce the model to be examined, discuss its properties, and introduce some pertinent notions of many-body thermalization necessary to follow the remainder of the work. We begin by preparing the ground for our investigation of eigenstate thermalization, followed by a definition of the model along the lines of [9,10]. Much of the material in this section is well known, and we refer the reader to the review [22] for further details.
We nevertheless include some background material in an effort to make the paper more self contained.
The manifestation of thermal behavior, defined as the applicability of equilibrium thermodynamics to isolated quantum systems at first seems puzzling from a microscopic perspective: how is a seemingly mixed thermal state reached, starting from a pure initial quantum state?
A common test case is the so-called quench scenario 2 , where a system is prepared in a farfrom-equilibrium state and its thermalization is studied as the system evolves in time. In classical mechanics the emergence of thermal behavior as the end point of such an evolution is a consequence of dynamical chaos, or more formally of ergodicity: the delocalization of a general initial state over phase space is understood as a consequence of dynamical chaos.
The end-point is then given as the state of maximal entropy, consistent with the values of its (few) integrals of motion at the initial time. In quantum systems, where the dynamics is linear and therefore not chaotic in the classical sense, the detailed mechanism is different, even though the notion of (quantum) chaos again plays a central role. A powerful framework for quantum thermalization is given by the eigenstate thermalization hypothesis [3,4,32] (ETH). In this scenario a quantum system thermalizes as a consequence of the structure of its eigenstates, rather than chaotic dynamics. In effect thermal behavior is already apparent at the level of individual eigenstates of the many-body Hamiltonian.

Eigenstate thermalization
The Eigenstate Thermalization Hypothesis [3,4] is an assumption on the properties of the spectrum of many-body Hamiltonians for systems showing thermal behavior. We start by stating the ETH ansatz for matrix elements of nonextensive operators 3 and then proceed to a discussion of its motivation and some consequences.

Complex SYK model
The model involves fermionic microscopic degrees of freedom on N sites. The spatial configuration of these sites is immaterial, due to the non-local nature of the interaction. In this paper we work with complex spinless fermions subject to the Hamiltonian [10,33] with complex coupling parameters J ij;kl . We define a set of fermion creation and annilihation operators, satisfying Antisymmetry of the fermions as well as Hermiticity of the Hamiltonian impose the constraints J ij;kl = −J ji;kl , J ij;kl = −J ij;lk , J ij;kl = J * kl;ij . (2.6) The remaining independent coupling constants in J ij;kl are drawn from a Gaussian distribution with zero mean. The variance sets the average strength of the coupling, which has to scale with the number of sites N as shown to ensure a well-behaved large-N limit. We will keep this normalization despite the fact that we always work at strictly finite N . The Hamiltonian, as written above (2.4) does not respect particle-hole symmetry. However, the particle-hole violating effects come from terms where two indices of J ij;kl take the same value and so they are suppressed by powers of 1/N . One could add a chemical potential that restores particle hole symmetry at any N , as in [33]. Since both versions connect to the same large-N limit we have chosen here instead to work with the simpler Hamiltonian (2.4).

Properties
The Hamiltionian (2.4) is a variant of a set of models considered previously in the context of quantum spin glasses [13], but crucially does not itself have a spin-glass phase. A closely related model has recently been proposed by Kitaev in terms of Majorana fermions at N sites [9].
The model is solvable 4 at large-N and can be shown to exhibit a number of striking properties [9,11,14], chiefly among them the fact that it exhibits maximal chaos, in the sense that an appropriately defined (see section 4.3) out of time order four-point function (OTOC) decays exponentially for times up to the scrambling time, where α is a coefficient, e.g. α ∼ βJ/N for the large-N limit of SYK [11]. This defines a quantum version of a Lyapunov exponent. The average is taken in the thermal state, as indicated. This Lyapunov exponent takes the maximal [15] value 2π/β in the SYK model, which is the same as that of a Schwarzschild black hole in Einstein gravity [9,39]. Here we provide numerical evidence that this Lyapunov exponent can also be extracted by considering instead energy eigenstates, that is to say correlation functions of the form The value of λ L in eigenstates can be meaningfully compared with the thermal value by appealing to the eigenstate thermalization hypothesis in order to associate a temperature T = β −1 to the individual eigenstate |E . To this end one may either follow [40] and associate an effective temperature β −1 with energy E via the canonical average Alternatively, for many of our computations, we determine all thermodynamic quantities in a microcanonical ensemble centered at the average energy E = (E n + E m )/2 for ease of comparison with the ETH ansatz.
Here we study (2.9) using exact diagonalization. It would clearly be interesting to establish some of our results analytically at large N , or perhaps via the large-D limit of [41]. Our work can be viewed as motivating the conjecture that λ L , defined in terms of eigenstates as in (2.9) will satisfy when evaluated at large N 1 and large coupling βJ 1, where β is defined as the inverse of the map (2.10) above. Let us now move on to a discussion of the methods and results of this work.

One point functions and eigenstate thermalization
The main analysis technique in this work is numerical diagonalization of the many-body Hamiltonian. For most applications we numerically diagonalize the Hamiltonian (2.4) for up to 17 sites 5 . This allows us to explicitly calculate the matrix elements of non-extensive operators made up of creation and annihilation operators involving a small number of sites. 5 We emphasize that this corresponds to a Hilbert space dimension 2 N = 2 17 which is the same as that of the Majorana SYK model with M = 34 sites, corresponding to a Hilbert space dimension of 2 M/2 . It is easily seen that total fermion number commutes with the Hamiltonian allowing us to work in sectors of fixed fermion number n F , denoting the filling fraction ν = n F N . This is useful numerically as it allows us to cut down the effective matrix sizes in the actual diagonalization process. For the most part we will work in the half-filling sector ν = 1 2 . If the number of sites N is odd, we mean ν = N +1 2N when we refer to 'half filling'.

On-diagonal terms are thermal
According to the ETH ansatz (2.1), diagonal matrix elements O nn = n|O|n are smooth functions of the average energyĒ, while off-diagonal elements are suppressed by the entropic factor e −S/2 .
We start by illustrating this exponential suppression in Figure (3.1), where it is easily seen that only the diagonal entries of the matrix are appreciable. We illustrate this behavior for N = 10 in Figure (3.1), but have checked it extensively for other accessible values of N finding excellent agreement with ETH expectations. For N = 10, one can make out the fluctuating nature of the off-diagonal matrix elements, which we will characterise precisely in section 3.2 below. Of course so far this is at the level of qualitative observation and we will turn to a more quantitive analysis of the off-diagonal matrix elements below.
Before we do so, let us analyse the on-diagonal matrix elements, O nn Ē in detail. The expectation for finite values of N is that the diagonal matrix elements condense more and more tightly onto a limiting smooth curve O(Ē) which is defined by extrapolation to the thermodynamic limit N → ∞. This is illustrated in Figure 3

Off-diagonal terms
Moving on to the off-diagonal terms of the matrix elements we demonstrate that the remainder terms are indeed well described by a Gaussian random distribution with zero mean and unit variance, in other words In    Figure 3.3, with a cross-over between RMT and non-RMT behavior at a characteristic energy E T , set by the coupling strength J. Despite the lack of spatial diffusion in the model we refer to this energy as the Thouless energy E T . Due to this lack of spatial structure (SYK is effectively a zero-dimensional model), the Thouless energy cannot be set by a diffusion time, and it is thus natural that it be set instead by the coupling strength J (at fixed energyĒ). The uppermost panel of to non-constant, i.e. non-RMT behavior. In higher dimensional local models this energy is often associated with the diffusive Thouless energy, but, as indicated before, such a physics interpretation appears not to be available in the zero-dimensional SYK case. We note that this corresponds well with the observations of [43], who previously presented evidence for This accords well with intuition as well as previous results in the thermal ensemble indicating chaotic properties to be most pronounced in the strong-coupling regime. Let us now discuss further chaotic aspects of the model.

Correlation functions and chaotic behavior
In section 3 we established the applicability of the ETH ansatz by studying one point functions of nonextensive operators in the complex SYK model. We will now study higher-point correlations in order to elucidate dynamical aspects of the model related to quantum chaos. Again, it is interesting to compare the behavior in a single random realization of the model versus the behavior of the same quantity after averaging over a large number of realizations.
We shall start, however, with the spectral form factor where the distinction between pure states and the thermal ensemble is meaningless, as can be seen from its definition in terms of an analytically continued partition function. One may also construct the spectral form factor as the fidelity of a certain pure state [46].

The spectral form factor
The spectral form factor is a well studied quantity in random matrix theory (see e.g. [47]) as it gives a clean probe of the eigenspectrum of a model. In particular its late-time behavior is sensitive to the discreteness of the spectrum as well as level statistics. The spectral form factor is most conveniently defined in terms of the analytically continued partition function,  In both cases we see the characteristic decay followed by linear ramp and plateau behavior. Left panel: β = 1. One sees several partial revivals in the decaying region with a power-law envelope. Right panel: β = 10. The red markers show points we used for a fit in the slope region, where we find a decay of ∝ t −4.53... , which is consistent with the value reported in [45]. also see the characteristic ramp and plateau behavior [44,46] at late times which, once more, looks qualitatively like RMT. The difference between RMT and the SYK model is to be found in the precise timescales of dip and plateau times [44]. For completeness and ease of comparison we discuss the RMT spectral form factor in some detail in appendix B.

The two-point function
Let us now turn to the study of correlation functions of non-extensive operators. Here we work with the two-site hopping operator h ij whose one-point functions are studied in Appendix A. The operator is defined as follows: pick two arbitrary sites i and j and write We have also verified that analogous results hold for the connected and full correlation functions of the on-site number operator n i .

Eigenstates
We now study correlation functions of the hopping operator in energy eigenstates. For definiteness we will take sites i = N − 1 and j = N , but any two sites will give essentially the same answer. We consider the two point function of the hopping operator, for some excited eigenstate |n with energy E n , as well as its connected cousin We find that the connected correlation function in eigenstates quickly decays and subsequently oscillates around zero as shown in Figure 4.2. This latter fact is easily established by averaging over time. One finds where we temporarily denoted the Hermitian operator h ij by O to avoid too much index clutter. It follows that the connected eigenstate correlation function averages to zero, G n c (t) = 0, at late times. The typical size of the late-time fluctuations follows from eigenstate thermalization, (2.1), to be ∼ e −S/2 . We show explicit computations of G n c (t) for different values of N in Figure 4.2 (top left).
As implied by (2.1), the behavior of the two-point function in eigenstates approximates very closely the corresponding microcanonical quantities. This agreement is expected to be perfect in the thermodynamic limit N → ∞. By a microcanonical two point function of an operator O we mean the quantity where N E is the number of states in a window of energies of given width ∆E around the average energyĒ and the sum over n runs over exactly those states. This agreement is illustrated in Figure 4.2 (top right). As a basic check one can convince oneself that (2.1) implies that its long time average gives exactly (4.5), which is also borne out in Figure  4.2 (top right). We conclude therefore that two-point functions in the disorder-averaged theory become arbitrarily close to their thermal (microcanonical) averages. This agreement is perhaps not surprising if one realizes that averaging the correlation function over couplings is operationally similar to a microcanonical average in the first place.
However, we now want to compare the behavior of G n and G n c , to the corresponding thermal correlation functions with respect to the canonical density matrix ρ = e −βH , at energy E(β), determined by the map (2.10), that is and its connected cousin G β (t) can differ.
The connected correlation function G n c (t) in eigenstates oscillates around zero at late times, unlike its thermal equivalent G β c (t), which oscillates around a non-zero average value as seen in Figure 4.2 (bottom right). These differences are subleading in the size of the Hilbert space and are expected to become negligible in the N → ∞ limit. The latter limit is of special interest for a putative gravitational dual as it corresponds to the semi-classical regime where a geometric description should become possible.
As a side comment, in the limit β → ∞, the thermal expectation value starts approximating the eigenstate one arbitrarily closely. This, of course, has nothing to do with thermalization, as it corresponds to the zero-temperature limit, where the 'thermal' average projects on the ground state, and is thus manifestly equal to the eigenstate correlation function. In conclusion then we find that the two-point function in individual eigenstates of

Superposition states
Up to this point we have considered mostly eigenstates. From what we have found one can conclude also that arbitrary pure states with narrow support in energy thermalize to microcanonical averages at late times, consistent with eigenstate thermalization. However, thermalization of states with broad support in energy do not thermalize in this way. We will next consider pure states, closely related to the ones considered in [29], with very broad spread over the energy spectrum 7 and demonstrate that they nevertheless thermalize, but more precisely to canonical averages. Note, again, that at finite N microcanonical and canonical averages do not have to exactly agree, and consequently one or the other may be a better approximation to a thermalizing pure state correlation function.
We consider pure states which are superpositions of eigenstates, as one would obtain, e.g. as a result of a sudden quench. An interesting class of such pure states can be constructed as follows. Let |C be a canonical state at half filling. Select a set of N C creation operators where |Ω is the state with no fermions. We assumed that N is even and we work at half filling ν = N C N = 1 2 , but such states can be constructed also for odd N and other filling fractions. What is important is that we think of these states at large N as being finitely excited, i.e. ν is held fixed as N → ∞. Now we define a one-parameter family of pure states via Euclidean evolution of the canonical state, where, according to ETH, the coefficients c α are Gaussian distributed complex random variables [4]. We expect these states to behave approximately thermally at a temperature β −1 that we can determine by the map (2.10), together with the fact that the eigenstates satisfy (2.1). To this end we compute The random nature of the coefficients c α ensures that this expectation value behaves like a thermal average. We can make this more precise by averaging the expansion coefficients over the eigenstate ensemble [4], using With the help of this expression, we find where we have averaged numerator and denominator independently. This quantity is equivalent to the thermal expectation value E β=2 since the state |C typically has support over the whole spectrum. It is thus natural to compare expectation values in |C with thermal expectation values in the canonical ensemble at β = 2 . The behavior of these states in a sense is similar to the thermofield double state, with the role of the trace over the second copy taken over by the random distribution of expansion coefficients. In Fig (4.2) we show the correlation function and its connected version G c (t), defined in the obvious way, in comparison with the analogous canonical averages. We see that in fact they are rather close to their thermal counterparts at inverse temperature β = 2 as expected. Let us next move on to four-point functions, an in particular the issue of chaos in eigenstates.

The four-point function
We now finally turn to studying the four-point function 8 which serves to characterize early time chaos via the Lyapunov exponent λ L defined in equation (2.9). Again we emphasize that this has been studied extensively in the thermal ensemble [9,11,15,33], while here our main focus is to study it in eigenstates.
We again focus on the two-site hopping operator h ij . For an SYK model on a grid of size N , let us define two operators together with their out-of-time-order four-point function The expectation value is taken, as before, in a finite-energy eigenstate, denoted F n (t), or 8 The results in this section were obtained in collaboration with Jérémie Francfort.
for comparison in the (micro-)canonical ensemble, denoted FĒ(t), F β (t), respectively. The expected form (see (2.11)) of this function at times up to the scrambling time is where the coefficient α = βJ/N in the canonical ensemble, at strong coupling and large N [15], the same circumstances under which the Lyapunov exponent takes its maximal value λ = 2π/β [9,15]. For our OTO correlation functions FĒ(t), which well approximates that in the corresponding eigenstate, we show a fit to the form (4.18) in Figure 4.3. Our results here are consistent with those of [33] who considered F β (t) in exact diagonalization and concluded that at finite N , the Lyapunov exponent is not maximal and does not vary inversely with β, but rather that the parameter governing λ is the coupling J. The behavior of the OTO correlator in eigenstates thus accords very well with expectations from eigenstate thermalization. In particular the early-time physics, including the scrambling time, of this correlator in eigenstates coincide to numerical accuracy with the thermal results. This suggests that the large-N OTO correlator in eigenstates F n (t) will also take the form (4.18) with Lyapunov exponent given by (2.11).

Discussion
In this work we have endeavoured to establish the mechanism of thermalization in the complex spinless SYK model, as a toy model of a strongly correlated quantum system with a holographic dual 9 . We focused on the complex model, but we expect the Majorana model to exhibit qualitatively similar behavior, that is to say that we believe it also satisfies eigenstate thermalization.
We were careful to establish eigenstate thermalization both for an individual random realization -with increasing accuracy as the Hilbert space dimension 2 N is increased -as well as for a fixed Hilbert space dimension at moderate N -with increasing accuracy as one averages over more and more realizations.
We have also studied the extent to which two-point and four point correlations in finiteenergy pure states approximate those in the thermal ensemble at the corresponding tem-perature. Our results support the conclusion that individual eigenstates in the SYK model behave thermally. We established that the agreement between pure and thermal expectation values becomes better for a single realization of the model as N is increased and for a fixed finite N , as we average over a larger and larger ensemble of Hamiltonians. However, consistent with earlier findings [44] we observe that correlation functions are self-averaging at early times, but lose this property at late times. This property is shared by the spectral form factor. Herein lies an important subtlety: a bulk dual 10 of SYK has been proposed for the disorder averaged theory, which means that any bulk solution is strictly dual to an ensemble of boundary Hamiltonians. One should therefore not associate a single eigenstate of an individual realization of the boundary theory with the late-time behavior of a bulk solution. This point does not apply to the tensor models of [23][24][25]. It should, however, be kept in mind during the rest of the discussion section.

Comments on putative bulk dual
Let us now address the question of the interpretation of our results in terms of a putative holographic dual, keeping in mind our previous comments on the status of such a dual.
A crucial issue in the holographic description of black holes is the representation of their interior from the boundary theory point of view [57]. One application of our results in this respect concerns the relationship between entanglement and geometry [58,59] (see also [60,61] for a world-sheet analog). One may appeal to the approach of [62] to argue that a typical highly entangled eigenstate of the SYK model does not have a dual with a smooth geometrical connection. The argument of [62] uses eigenstate thermalization as a hypothesis to roughly reason as follows. We note that a typical two-sided correlation function in the eternal black hole geometry will be of order e S at early times, coming from the wormhole connecting the two boundaries [58,59]. One then appeals to the eigenstate expectation values of the form (2.1) to argue that the same operator in a generic highly entangled finite energy state does not have the required e S correlations at early time, in fact it is exponentially suppressed. This way one arrives at a contradiction with the assumption that the correlator can be computed in a smooth geometry with a wormhole connection between the two boundaries. Closely related ideas have been advanced in [63]. By establishing eigenstate thermalization in the SYK/NAdS 2 context, one important implication of our 10 See previous footnote.
work is that a generic highly entangled state of the SYK model either does not have a smooth geometric dual, or that entanglement does not generically correspond to having a geometrical wormhole in the putative bulk dual of SYK. However, without entering into the details, if one allows the state-dependent construction of interior operators by [64], smooth black hole interiors may be constructed.
Of course more directly eigenstate thermalization tells us that one-sided correlation functions look thermal even in eigenstates. This means that correlations in individual eigenstates are well approximated by dual computations in the black hole geometry. Similar results apply in AdS 3 /CFT 2 see e.g. [5,65], where two-point functions of light operators in states created by heavy primary operators were shown to be well approximated by the corresponding results in the BTZ black hole and non-equilibrium initial states thermalize exactly to this state [7,8].
In higher dimensions a direct approach, such as in this paper, seems out of reach. It would therefore be interesting to gain a more analytical understanding of our results, and we hope to address this in the future. It will be interesting to try and carry our results over into a more widely applicable picture of thermalization in theories with holographic duals. In this respect it may be interesting to map out the applicability of eigenstate thermalization in more SYK-like models, such as [68][69][70][71]. Conversely, if one instead accepts that eigenstate thermalization should hold in theories with holographic duals, one might hope to use the constraints on matrix elements due to eigenstate thermalization, in order to further refine the requirements on CFTs with a well-behaved holographic dual.
It would be interesting to repeat the study in the present paper using the tensor models of [23,24,72] and in particular to study how much of what we uncover here survives away from the large-N semi-classical limit in such theories, which have the advantage of being defined without the quenched disorder average. Certain spectral and chaotic properties have been studied via exact diagonalization by [73].
In conclusion we believe that the detailed study of thermalization via eigenstates in SYK, both numerically and analytically, gives us a concrete opportunity to better understand the physics of quantum black holes, at least at the level of toy models.

Acknowledgements
We would like to thank Bela Bauer for his early invaluable help in writing an efficient numerical algorithm for the complex SYK model. We thank Dmitry Abanin, Shouvik Datta, Thierry Giamarchi, Wen Wei Ho, Pranjal Nayak, Kyriakos Papadodimas, Eliezer Rabinovici, and Benjamin Withers for discussions. Our special thanks go to Jérémie Francfort for his collaboration at an early stage. The simulations in this paper were performed on the baobab HPC cluster at the University of Geneva and we thank Yann Sagon for his assistance.

A Full spectrum & ETH for other operators
This appendix is concerned with filling in some details on the spectrum of the model, as well as to supply more details on eigenstates thermlization for a different non-extensive operator, namely the hopping operator. A careful analysis of the spectral properties of the SYK model  was carried out in [11,14,43,44,74]. Here we present some essential features on the complexspinless case, to set the context, and also to benchmark our numerics. More details are presented in the aforementioned references.

A.1 Density of states
The full spectrum of the model is most efficiently computed by considering each allowed filling fraction ν separately. It is both of interest to consider the spectrum of a single randomly chosen realization as well as averaged over a large number of realizations (see. Figure A.1).
As has been observed before, for the Majorana model, the spectrum self averages very well even for moderate values of N as can be surmised from comparing left and right panels of

A.2 ETH for the hopping operator
Everything we said about thermalization in eigenstates should hold for any non-extensive operator in SYK model. In order to illustrate that this is indeed the case, we collect here some illustrative results for an operator which differs considerably from the on-site number operator, namely the two-site hopping operator.
A further Hermitian operator of interest is the two-site hopping operator h ij . It is defined by fixing two arbitrary sites i and j and then writing One might think that the simpler operator c † i + c i would also have been a possible choice, but it is easy to see that does not conserve fermion number and so its matrix elements vanish in a fixed fermion sector as considered in this work.
We have extensively studied the matrix elements of the hopping operator, finding that they also satisfy the ETH ansatz. Fig A.2 illustrates the exponential suppression of offdiagonal matrix elements by an entropy factor. The subtlety for the hopping operator is that its thermal value, i.e. its on-diagonal matrix elements, is actually zero. This is shown in Figure A

B Random matrix theory
In this appendix we collect some results on random matrix theory, which have been referred to occasionally in the main text. These serve as a reference for the behavior of the complex spinless SYK model, which, as we explained, shows RMT-like behavior for some parameter and energy ranges. The RMT results also served as helpful test cases to verify our algorithm with the aid of known results. practice by taking the average over a large number of individual draws from the ensemble.

B.2 Spectral form factor
The RMT spectral form factor has recently been studied by various groups [44,46,75] as a reference and illustration for certain aspects of the SYK case, which are qualitatively well captured in RMT. For convenience we also present this quantity here, focusing on the GUE.
In RMT, more specifically the GUE, it is actually possible to analytically calculate the spectral form factor [46,47] using the method of orthogonal polynomials. Let L be the Hilbert space dimension, that is to say we consider the ensemble of L × L Hermitian matrices. For SYK the Hilbert space dimension was 2 N so, when comparing the two one should obviously set 2 N = L. If we define ν = β + it, the answer takes the form [46] S(β, t) = g(β, t) Z(β) 2 (B.2) with g(β, t) = e 1 4 (ν 2 +ν 2 ) (g c (β, t) + g d (β, t)) + Z(2β). One can clearly appreciate the qualitative ressemblance to the SYK case. A detailed discussion of time scales and various power laws can be found in [46].