Eigenstate Thermalization Scaling in Majorana Clusters: from Chaotic to Integrable Sachdev-Ye-Kitaev Models

# Eigenstate Thermalization Scaling in Majorana Clusters: from Chaotic to Integrable Sachdev-Ye-Kitaev Models

## Abstract

The eigenstate thermalization hypothesis (ETH) is a conjecture on the nature of isolated quantum systems that guarantees the thermal behavior of subsystems when it is satisfied. ETH has been tested in various forms on a number of local many-body interacting systems. Here we examine the validity of ETH in a class of nonlocal disordered many-body interacting systems — the Sachdev-Ye-Kitaev Majorana (SYK) models — that may be tuned from chaotic behavior to integrability. Our analysis shows that SYK (with quartic couplings), which is maximally chaotic in the large system size limit, satisfies the standard ETH scaling while SYK (with quadratic couplings), which is integrable, does not. We show that the low-energy and high-energy properties are affected drastically differently when the two Hamiltonians are mixed.

## I Introduction

For well over a century it has been understood that physical observables of a system coupled to a heat bath in the long time limit are such that they may be computed from a statistical average over all states of the system with a weight that depends on only a few parameters including the temperature of the bath.

Since heat baths are idealized objects coupled to but distinct from the system of interest, a natural question is: under what circumstances can a subsystem of a closed (isolated) quantum system be thermalized by the remainder, in the sense of physical observables being expressed in terms of a trace over some density matrix with a handful of constraints including particle number and temperature? This led to the eigenstate thermalization hypothesis (ETH) which is a conjecture about the nature of matrix elements of physical observables that, if true, reconciles the predictions of thermodynamics (with respect to an ensemble whether microcanonical, canonical or grand canonical) with those of quantum mechanics in the long-time limit Srednicki (1994); Deutsch (1991); D’Alessio et al. (2016).

The essence of ETH is that diagonal matrix elements of physical observables in the energy eigenstate basis, that is distinguished by time evolution, are expected to be quasi-randomly distributed with a width that falls off rapidly with the Hilbert space dimension leading to a smooth variation of the observable with energy. The off-diagonal matrix elements are also expected to fall off rapidly with increasing system size which can be viewed as a de-phasing of the subsystem by the remainder of the system. Taken together, these conditions on the matrix elements ensure that subsystems of closed quantum systems reach steady states that correspond to thermal equilibrium.

Advances in the study of optically trapped cold atoms in recent years have made the study of reasonably well isolated quantum systems experimentally feasible to the extent that ETH has transformed from a foundational issue in statistical mechanics into a matter of practical significance. With such an experimental impetus, there is now a sizable literature exploring ETH mainly using numerical techniques on finite clusters in a variety of many-body local interacting quantum systems Rigol et al. (2008); Rigol (2009); Biroli et al. (2010); Rigol and Santos (2010); Santos and Rigol (2010); Neuenhahn and Marquardt (2012); Beugeling et al. (2014); Brandino et al. (2012); Ikeda et al. (2013); Steinigeweg et al. (2013); Beugeling et al. (2015); Sorg et al. (2014); Kim et al. (2014); Mondaini et al. (2016); Alba (2015); Mondaini and Rigol (2015); Lashkari et al. (2016); Nandy et al. (2016); Luitz and Bar Lev (2016); D’Alessio et al. (2016); Chandran et al. (2016); Mondaini and Rigol (2017). There is by now overwhelming evidence that ETH is obeyed over most of the many-body spectrum in systems with local interactions, with important exceptions close to integrability and in many-body localized states.

Despite the widespread success of ETH in such numerical investigations, a number of questions remain open. In this paper, we extend studies of ETH to a class of many-body interacting models that are both disordered and nonlocal. One member of this class is the Sachdev-Ye-Kitaev model, SYK, which consists of Majorana fermions with infinite range random four-point couplings Kitaev (2015); Maldacena and Stanford (2016). This model recently shot to prominence having been shown to exhibit a number of remarkable properties Kitaev (2015); Maldacena and Stanford (2016) including, in the large and strong coupling limit (i) an extensive zero temperature entropy (ii) an approximate emergent reparametrization invariance at large and (iii) an out-of-time-ordered correlator with a Lyapunov exponent that saturates a conjectured bound Maldacena et al. (2016) implying that it is “maximally chaotic”.

A possible intuition for the dynamics of such zero-dimensional disordered models is that the non-locality allows information to propagate rapidly across all the sites while the randomness in the couplings translates to dynamics that approximates to the action of random unitaries that should scramble information rapidly. This intuition is borne out by a numerical study of the quench dynamics Eberlein et al. (2017) which shows that SYK approaches equilibrium at some temperature when the initial state is the equilibrium state for some other model. However, this intuition fails for the close relative SYK which is a random quadratically coupled model for which the Lyapunov exponent vanishes García-García et al. (2017). The presumption that scrambling and thermalization should be connected suggests that SYK should obey ETH while SYK should not. This is the question that we explore in this paper. Previous work has showed that the Majorana version Hunter-Jones et al. (2017) and the complex fermion analog of the SYK model Sonner and Vielma (2017) satisfy ETH. Here, we explore the Majorana version, which splits into different sectors depending on . We tune parametrically between the SYK and SYK models.

Like Refs. Sonner and Vielma, 2017; Hunter-Jones et al., 2017 and the usual ETH literature, we first present a study of the middle of the many-body spectrum, i.e., excluding the spectral edges. We carry out a finite size scaling analysis that tests a strong (scaling) version of ETH Beugeling et al. (2014). Previous work has argued that the complex fermion analog of SYK obeys ETH Magán (2016). Our numerical and analytical work demonstrates scaling behavior of the SYK matrix elements that contrasts sharply with the scaling behavior of ETH but which is very similar to the behavior observed in other integrable models. The ETH scaling is valid also at parameters intermediate between the SYK and SYK points, which is expected as the model becomes non-integrable as soon as quartic coupling is included, i.e., everywhere except at the SYK point.

In addition to these results pertaining to the middle of the spectrum, we show a peculiarity of the low-energy part of the spectrum, which is more relevant to the holography literature. The low-energy eigenstates of SYK show the ETH behavior, which is not expected in local models but may occur in nonlocal models because the low energy eigenstates exhibit volume law entanglement Liu et al. (2017). The low-energy eigenstates of SYK show the scaling typical of integrable models. However, at intermediate points, the ETH scaling is no longer seen — an admixture of SYK destroys ETH scaling at low energies. We believe this can be traced to the fact that the SYK is an attractive fixed point in the RG sense. In other words, ETH scaling of low-energy eigenstates is not governed by non-integrability but by the RG flow.

The paper is organized as follows. In the next section, we outline the statement of ETH, highlighting the scaling form (size-dependence). Section III has a discussion of the symmetries of SYK and SYK and the constraints they place on the matrix elements of different classes of “local” operators. Sections II and III are mostly pedagogical reviews, meant for readers unfamiliar with ETH scaling or with the finite-size properties of SYK clusters. Our main results are in Sections IV, V and VI. In Sections IV and V we study the scaling behaviour of the diagonal matrix elements and, in Section VI, the off-diagonal matrix elements. We then conclude with a discussion of some of the significance of our findings.

## Ii The Eigenstate Thermalization Hypothesis

ETH is a generic mechanism for thermalization of observables in a closed, finite, non-integrable system out of equilibrium, i.e., a mechanism for why the long-time average value of an observable should be describable by a Gibbs ensemble. The class of operators for which ETH is expected to be valid includes local operators. ETH is a statement about the diagonal matrix elements of the operator in the energy eigenstates of the Hamiltonian, i.e., the eigenstate expectation values . ETH states that the diagonal matrix elements of are smooth functions of the energy in the case when the system size is large, with the fluctuations being exponentially small in the system size. This follows from the idea that high-energy eigenstates of non-integrable systems are complicated enough that the eigenstate coefficients are effectively random. For systems with a finite Hilbert space dimension , the central limit theorem yields the scaling to be Neuenhahn and Marquardt (2012); Beugeling et al. (2014):

 Onn=f(1)O(En)+1D1/2f(2)O(En)Rnn (1)

where is a pseudo-random number with unit width and are smooth functions of the eigenenergy .

A further important ingredient of ETH concerns the off-diagonal matrix elements: for . These matrix elements should be suppressed exponentially in the system size. Use of the central limit theorem again shows the scaling to be :

 Omn=1D1/2f(3)O(ϵ,ω)Rmnform≠n (2)

where is a pseudo-random number with unit width, , and . Together, these are usually written together in the following form:

 Omn=f(1)O(En)δmn + e−S(E)/2fO(ϵ,ω)Rmn. (3)

This common form is more general than Eqs. (1), (2), because it applies also to systems with unbounded Hilbert spaces. For systems with finite Hilbert spaces, the entropy in the middle of the spectrum scales as , so that the two forms are equivalent.

Although this is not commonly stressed, ETH is a scaling statement. The idea of eigenstate coefficients being effectively random leads unavoidably to the scaling. Thus, a system in which the state-to-state fluctuation of diagonal matrix elements decreases polynomially with system size would not be considered as satisfying the ETH, as formulated commonly through Eq. (3). Polynomial scaling has been observed in integrable systems Ziraldo and Santoro (2013); Alba (2015); Nandy et al. (2016), which are considered not to obey the ETH. Also, diagonal matrix elements being equal to the canonical value in the infinite-size limit (used as a definition of ETH in Ref. Magán (2016)) is not equivalent to the usual statement of ETH, outlined above.

For local Hamiltonians, the ETH is expected to be valid only in the middle of the spectrum (high-energy states). Since we are studying non-local Hamiltonians in this work, we will also examine whether it is valid for low-energy states.

The physical content of the ETH Ansatz becomes clear when we consider computing thermodynamic observables as long time averages. For example, consider on state :

 O(t)=∑m|cm|2Omm+∑m,nm≠nc⋆mcnei(Em−En)tOmn.

The long-time average of this operator expectation value is

 ⟨O⟩≡1T∫T0dtO(t)=∑m|cm|2Omm

which is the trace of the operator computed over the so-called diagonal ensemble Rigol et al. (2008). One may show that this is equivalent to the microcanonical average if ETH is satisfied, assuming that the state has mean energy with subextensive fluctuations. The condition of small off-diagonal elements ensures that the time average of converges to its thermal value in the long time limit.

ETH in the form stated above has been tested and found to be valid on a range of local many-body interacting models using, primarily using full spectrum numerical diagonalization on finite systems. The effect of tuning parameters toward an integrable point has also been widely explored Rigol (2009, 2009); Rigol and Santos (2010); Beugeling et al. (2014); Sorg et al. (2014). As the couplings are tuned towards an integrable point, there is a continuous departure from the conditions of ETH, e.g., a broadening of the distributions of matrix elements. The significance of the breakdown of ETH close to integrability is that integrable points have an extensive number of constants of the motion that prevent thermalization to the usual thermodynamic ensembles. Instead, such systems are conjectured to thermalize to a generalized Gibbs ensemble Rigol et al. (2007) that enforces the extensive number of integrals of the motion in such systems. Another class of systems that do not respect ETH are Anderson localized or many-body localized systems in which disorder inhibits thermalization Nandkishore and Huse (2015).

In the following, we address the question of whether ETH is respected in a class of models that are simultaneously fully connected, so that the concept of locality is absent, and intrinsically disordered. We test ETH scaling, as described above, for SYK and SYK.

## Iii SYK Models at Finite N

In this section, we introduce the SYK and SYK models and their symmetries. We introduce the few-body operators whose matrix elelments we will study, and show how the model symmetries place constraints on these matrix elements.

### iii.1 Hamiltonian and Symmetries

 H=(cosθ)HSYK4+(sinθ)HSYK2 (4)

where

 HSYK4=∑1≤i

and

 HSYK2=∑1≤i

Here ’s are Majorana operators. The couplings and are gaussian distributed random variables with mean zero and respective variances and . The parameter runs from to , with and being the SYK and SYK points, respectively.

The anticommutation relation satisfied by Majoranas, , is identical to the Euclidean Clifford algebra up to normalization. Therefore, we may represent the site Majoranas with Hilbert space dimension by gamma matrices of size . The particular representation we use for numerical analysis is given below in Section III.3.6. We introduce the parity operator and note (i) that (ii) that anticommutes with the for a given . Property (ii) implies that the SYK Hamiltonian commutes with so the parity is a good quantum number.

We now consider the time reversal operator : an anti-unitary operator which may be written as the product of a unitary operator and the complex conjugation operator . The gamma matrices obey . Our expectation is that complex fermions map to under time reversal, which implies that . This is ensured by choosing

 UT=PM+1∏m=1Γ2m−1. (7)

### iii.2 Observables

Our investigation into the validity of ETH in SYK models rests on calculations of matrix elements for particular few-body operators. The operators we consider are

 Aij≡iχiχj

and

 Bijkl≡χiχjχkχl

where the Majorana labels should be distinct but are otherwise arbitrary.

### iii.3 Periodicity of SYK4

The SYK Hamiltonian is a matrix composed of independent random nonzero elements but with elements. The matrices are thus sparse, as is typical for many-body Hamiltonians. As is the case with usual condensed matter Hamiltonians, it can still be useful to compare the symmetries with dense random matrix ensembles (Wigner-Dyson classes).

As pointed out first in Ref. You et al. (2017), the nature of the spectrum and matrix elements of the SYK model have a periodic dependence on , because the underlying symmetries correspond to different Wigner-Dyson classes (GOE, GUE, GSE), depending periodically on . The periodic correspondence is shown in the table on page 1 for up to 30. This correspondence to the random matrix classes has been verified numerically, e.g., through the level spacing statistics You et al. (2017); García-García and Verbaarschot (2016); Cotler et al. (2017) and from the spectral form factor , where is the partition function Cotler et al. (2017). Both measures exhibit departures from random matrix predictions above some Thouless energy scale as is common in many-body models Beenakker (1997). In the spectral form factor, ramp and plateau features occur in both random matrix theory and SYK with the Thouless energy corresponding to short time scales than the ramp onset Cotler et al. (2017). Below, we review the association between the SYK model of different sizes to the Wigner-Dyson classes. We also present numerical data illustrating features complementary to those in the earlier literature.

Because the Majoranas are invariant under time reversal, so too is the SYK Hamiltonian, . By direct calculation one may show that for mod and for mod . It also follows that for mod and for mod . So there are four distinct cases to consider based on the and . We always work in the basis with

 P=(100−1).

To see how all the symmetries constrain the spectrum, we now show that the algebraic relations on and may be satisfied by a particular matrix within each class that then fixes the form of the Hamiltonian.

#### N=0 mod 8

This class is specified by the conditions and . One can find a representation for the gamma matrices for which

 UT=(100−1)

which fulfils the algebraic relations on and . It follows from time reversal invariance, , that the Hamiltonian takes the form

 H=(H(1)R00H(2)R).

Each block is a distinct real random matrix. Therefore, we expect gross properties of the eigenstates of each block to correspond to those of the GOE random matrix ensemble.

#### N=2 mod 8

In this case,

 UT=(0110)

fulfils the conditions and and then implies that the Hamiltonian takes the form

 H=(HC00H∗C).

The two chirality sectors thus have identical GUE spectra. In other words, the eigenvalues are doubly degenerate with distinct chirality eigenvalues.

#### N=4 mod 8

For this class,

 UT=(Ω00−Ω)

ensures that and where which may be, for example,

 Ω=(01−10).

Then

 H=(H(1)H00H(2)H).

where and are each a matrix of real quaternions represented by blocks of complex numbers of the form

 (z1z2−z⋆2z⋆1).

Random Hamiltonians built from real quaternions belong to GSE by definition. The condition enforces a double degeneracy on the eigenstates in a similar manner to the Kramers degeneracy of time reversal invariant half odd integer spin systems. The blocks have distinct spectra, but the eigenvalues within each sector are doubly degenerate.

#### N=6 mod 8

The fourth and final case has and so that

 UT=(01−10)

from which it follows that

 H=(HC00H∗C).

So, the properties are the same as for the case mod .

#### Symmetry constraints on diagonal matrix elements

Time reversal places constraints on certain observables. In particular, for the GOE ensemble, since and , the diagonal matrix elements vanish, . For the classes with double degeneracy , ensured by time reversal symmetry, one finds

 ⟨+|iχiχj|+⟩=−⟨−|iχiχj|−⟩

Some of these properties of the SYK spectrum and matrix elements are illustrated in Fig. 1 and summarized in the table.

Referring to Fig. 1, we observe the following features. For , the GOE case, energies are singly degenerate and the two-point operator matrix element vanishes as shown above. For the other three system sizes, the spectra are doubly degenerate and the two-point matrix elements for the pairs are equal and opposite. For the GSE class, the double degeneracy arises for states with identical parities as illustrated by the up-down reflection symmetry in the upper panel. For the GUE classes, the doubly degenerate states have opposite parities so the upper panels have up-down reflection symmetry up to a swap in the symbols (open and closed symbols swap). For the four-point function matrix elements for doubly degenerate pairs are the same so for GUE the different parity symbols are overlaid and in the GSE case identical parity symbols are overlaid.

#### Numerics on SYK Models

We may build a particular gamma matrix representation starting from the Pauli matrices , and , the first two of which represent two Majoranas. Then, given the gamma matrices for Majoranas together with the operator we may obtain a representation for Majoranas

 Γ(N+2)i =σ1⊗Γ(N)ii=1,…,N Γ(N+2)N+1 =σ1⊗P(N) Γ(N+2)N+2 =σ2⊗I.

In this representation,

 P=(100−1)

and the Hamiltonian is block diagonal with each block giving states of distinct .

### iii.4 Syk2: Operators and periodicity

We now consider the SYK model. Once again, the Hamiltonian is block diagonal in the gamma matrix representation given above because it commutes with . SYK is not time reversal invariant but instead has the feature that implying that the spectrum has paired eigenvalues with . Despite the lack of time reversal symmetry, the periodicities in mod are present but less distinctive than in SYK. Similar arguments to those for SYK reveal that all the Hamiltonians have complex entries but mod have pairs arising with the same parity while mod have pairs coming from different blocks. In SYK, time reversal antisymmetry places the following constraints on matrix elements

 ⟨+|iχiχj|+⟩ =−⟨−|iχiχj|−⟩ ⟨+|χiχjχkχl|+⟩ =⟨−|χiχjχkχl|−⟩

where, here, have energies with opposite sign and the site indices on each of the operators are distinct. These features of the matrix elements are illustrated in Fig. 2 where denotes some two point Majorana operator and a four-point Majorana.

The SYK model has a Poissonian level spacing distribution as is typical of integrable models. This is most easily understood in the language of complex fermions. In the single particle sector, the hopping model is simply a random matrix in GOE with a semi-circle density of states. The multi-particle sectors are straightforwardly obtained from the single particle sector and the resulting multiparticle density of states should be gaussian as an application of the central limit theorem. For a given window of energies and at sufficiently high energies, the occupation number distribution is effectively random because there is no mutual level repulsion in this integrable case so the level spacings are Poisson distributed.

### iii.5 Combined Hamiltonian (SYK4+Syk2) — spectrum and level statistics

Fig. 3 shows the evolution of the spectrum for different small system sizes as coupling is tuned from corresponding to SYK to which is SYK. The parity symmetry is common to models at all . For , the three system sizes considered correspond to GSE (), GUE () and GOE () which respectively have double degeneracies within each parity sector, double degeneracies within opposite parity sectors and only accidental degeneracies. As increases, the energies flow to the symmetry at which is present in identical parity sectors for and opposite parity sectors for . For intermediate there are direct level crossings between different parity sector eigenstates and avoided level crossings between identical parity sector eigenstates as one would expect.

Further connections between random matrix classes and SYK models are drawn out through a calculation of the mean ratio of consecutive level spacings Oganesyan and Huse (2007). This measure is based on , the set of level spacings in an ordered list of eigenenergies. The ratio

 rn=min(sn,sn−1)max(sn,sn−1)

is defined for each pair of consecutive level spacings. For Poisson statistics, the probability distribution of is with mean . For the Wigner-Dyson ensembles, we have and Atas et al. (2013).

Fig. 4 shows as a function of the coupling . The spectrum of the odd parity sector is used in each case, and values are combined from a large number of disorder realizations. The value of corresponds to Poisson statistics at regardless of Wigner-Dyson class. At , the value for matches our expectation that it belongs to GOE. A similar result holds for the GUE cases . For the GSE cases at , there is a double degeneracy (the Kramers degeneracy) within each parity sector. Thus every second level spacing is zero, leading to .

As increases from zero, the statistics in all cases match those of GUE because the matrix elements are irreducibly complex and the degeneracies are lifted. For larger there is a smooth crossover to Poissonian statistics. The Poissonian to GUE crossover near the SYK point has previously also been presented in Ref. García-García et al. (2017). The -statistics has also been described for the SYK point in Ref. You et al. (2017) and for a coupled chain of alternating SYK and free-Majorana ‘sites’ in Ref. Jian and Yao (2017).

Comparing different system sizes (e.g., and for the GSE class) shows that all the crossovers become steeper for larger system sizes. At larger , the deviation from GUE behavior gets increasingly confined to the singular points and .

## Iv Diagonal Matrix Elements: High-energy states

In this section we report a study of the eigenstate expectation values (diagonal matrix elements) in the eigenstates in the middle of the many-body spectrum. (We refer to the middle of the spectrum as high-energy states.) ETH scaling in the low-energy part of the spectrum is discussed in a later section.

As usual in studies of the ETH, it is important to restrict to a single symmetry sector. We present results for the odd-parity sector.

### iv.1 Non-integrable behavior of SYK4

We begin with SYK. As we have seen, various properties of SYK have a periodic dependence on the system size that can be seen to follow from a corresponding random matrix ensemble. We will present results for within a common ensemble. The GOE systems have the disadvantage that the operator is identically zero for all eigenstates, while the GSE systems have additional symmetry sectors within the odd-parity half of the Hilbert space. In addition, there are more numerically accessible finite-size instances of GUE than either of the other two. Hence we present results for the systems with GUE symmetry, i.e., for the sequence .

We first consider the diagonal matrix elements of the two-point operator for some arbitrary choice of and . Fig. 5(a) shows the diagonal matrix elements for the full spectrum from exact diagonalization plotted against energy for four different system sizes falling into the GUE ensemble. The matrix elements are observed to be distributed around zero for a given disorder realization independent of the location of the eigenstate within the spectrum. The distribution has a clear energy dependence, at least for the and , being broader at the extremes of the spectrum. The distribution of the matrix elements irrespective of energy is gaussian and the width of the distribution narrows as increases, as shown in Fig. 5(b).

As explained in Section II, ETH is a scaling statement — the width of the diagonal matrix elements should scale as . In the SYK models, . Numerical results for six system sizes, shown in Fig. 5(d), shows that the width does indeed scale as for SYK. (For a single parity sector, the relevant Hilbert space is half of , i.e., . However, this is a change of overall factor and does not affect the scaling form, .)

There is some shot-to-shot variation of the width as calculated from individual realizations. Fig. 5(c) shows the distribution of values obtained from different disorder realizations. The width of these distributions falls off rapidly with system size.

### iv.2 Integrable behavior of SYK2

We now consider SYK, which is an integrable (quadratic) model and hence is not expected to follow ETH scaling. We consider diagonal matrix elements of the two-point Majorana operator (operator ) and four-point operator (operator ). Numerical diagonalization shows a scatter in the matrix elements centred on zero (top panels of Fig. 6) for both operators and the distribution does not apparently narrow with increasing system size. This is qualitatively consistent with the behaviour observed in other integrable models. A more quantitative picture is seen by explicitly plotting the width against system size, as in the lower panel of Fig. 6. This clearly shows that scaling (i.e., scaling) does not hold. Different disorder realizations reveal a large scatter in the width for each system size as shown by the different green points at each . The average of all of these appears to fall off with a power law , in strong contrast to standard ETH scaling. For the operators we find scaling.

We may derive the scaling of the width of the diagonal matrix element distribution: for the operators and for the operators. We may pair the Majoranas of SYK arbitrarily and here we choose complex fermion to be composed of Majoranas and so that and . The SYK Hamiltonian can be seen to contain number conserving and non-conserving terms. As the Hamiltonian is quadratic we may diagonalize via

 cI=∑α(uαIdα+uα+NId†α)

so that

 HSYK2=∑αϵαd†αdα

and the many-body eigenstates are

 |ψ{m}⟩=(d†1)m1(d†2)m2…(d†N/2)mN/2|0⟩ ≡|m1,m2,…,mN/2⟩

corresponding to eigenenergies . The diagonal matrix element may be chosen without loss of generality to be so that the complex fermions belong to the same site

 Missing or unrecognized delimiter for \right +(uα+(N/2)⋆Iuα+(N/2)I−uαIuα⋆I)dαd†α}

In the spirit of exploring fluctuations between different eigenstates within the same disorder realization as is typical in the ETH literature, we make the assumption that the single-particle eigenstate coefficients are Haar-random variables. Then

 ⟨⟨uαIuβ⋆J⟩⟩=2NδIJδαβ ⟨⟨uαIuβJuγKuδL⟩⟩=4N2−4(δIKδJLδαγδβδ+δILδJKδαδδβγ) −8N(N2−4)(δIKδJLδαδδβγ+δILδJKδαγδβδ)

where the double brackets indicate an average over the Haar measure. It follows that that is, the diagonal two point matrix elements for SYK fluctuate around zero mean. We may compute the variance of the distribution under the same randomness assumption with the result that this varies as for large . This is consistent with our numerically observed scaling for the standard deviation. A similar calculation for the four-point operator shows that and . Note that the above results may also be regarded as coming from a disorder average over the SYK couplings for a fixed energy window.

Taken together, the numerical results we and others have obtained for SYK — the Poissonian statistics, vanishing of large number of matrix elements and power law fall-off of the spread in the diagonal matrix elements — are all consistent with results obtained for integrable many-body models D’Alessio et al. (2016); Vidmar and Rigol (2016); Ziraldo and Santoro (2013); Alba (2015); Nandy et al. (2016). At integrability, the number of conserved quantities is extensive and these are thought to bring about thermalization to a generalized Gibbs ensemble, for a wide range of observables and initial conditions, coinciding with a failure of ETH scaling. We note that Ref. Magán, 2016, which claims ETH in complex SYK, is based on a different and non-standard statement of ETH as discussed in Section II. We discuss thermalization at integrable points further in the concluding section VII.

### iv.3 Between the SYK4 and SYK2 points

Our general Hamiltonian, , is not integrable, and therefore the middle of the spectrum should show standard ETH scaling and generic thermalization behavior. The behavior is demonstrated in Figure 7, top center. For large enough systems one expects this to be true for all values of , i.e., excepting the integrable point . At points closer to , one expects the behavior to set in at larger sizes. This is the typical behavior close to integrability, as explored, e.g., in Ref. Beugeling et al. (2014).

## V Mid-spectrum versus low-energy matrix elements

Fig. 7 shows the scaling of the standard deviation of the distribution formed from the operator matrix elements at and points. The upper and lower panels show, respectively, drawn from the middle quarter and lower th of the full spectrum.

Unlike Hamiltonians with local interactions, the SYK models (with ) have the property that low-energy states in the large- limit are ‘thermal’ and even maximally chaotic Kitaev (2015); Maldacena and Stanford (2016); Kourkoulou and Maldacena (2017). In finite-size SYK systems, this is manifested, for example, in the ground-state entanglement entropy between two halves of the system growing linearly with system size, i.e., ‘volume law’-like behavior instead of ‘area law’-like behavior Liu et al. (2017). It is therefore natural to conjecture that the low-energy eigenstates might display ETH scaling, i.e., that the low-energy part of the spectrum might show behavior typical of effectively random states. Fig. 7(b) shows ETH scaling for the low-energy eigenstates of SYK. While the asymptotic behavior is apparently , the behavior sets in at larger sizes than for the middle of the spectrum. For local Hamiltonians, the low energy states show quite different behavior from typical states in the middle of the spectrum - the low energy states tend to have area law entanglement entropy and ETH ceases to hold. The SYK model is exceptional in this sense since the ETH scaling is present throughout the spectrum.

For (SYK), shows the previously reported scaling, not only at the middle of the spectrum as reported in Fig. 6, but also at low energies, Fig. 7(e,f).

The most interesting behavior is seen for the combined Hamiltonian, as exemplified by the data in Fig. 7(c,d). At the point, in the middle of the spectrum both operators exhibit the scaling we have come to expect from the SYK model; indeed is the expected behavior for any non-integrable system. However, the low-energy behavior is determined not by the non-integrability but rather by the RG properties. The scaling dimension of the fermion operator in SYK is , so the quadratic SYK perturbation has positive mass dimension and is RG relevant. This implies that low energy observables will behave differently to expectations coming from SYK when the SYK coupling is present. Indeed, Fig. 7(d) shows that, when we consider only states from the lowest th of the spectrum, the scaling of the widths departs considerably from the SYK prediction , consistent with the expectation from the RG argument above. Thus, for , the high-energy physics is governed by non-integrability, while the low-energy physics is governed by flow to the integrable fixed point.

## Vi Off-Diagonal Matrix Elements

In this section, we discuss the off-diagonal matrix elements, Fig. 8, for the operators and . The off-diagonal matrix elements have arbitrary phase, since the eigenstates each are defined only up to a phase. On average, this guarantees that the real and imaginary parts of the off-diagonal elements will be distributed symmetrically around zero.

For SYK, the distribution of off-diagonal elements is gaussian, as in other non-integrable models Mondaini and Rigol (2017); Beugeling et al. (2015). As for the case of diagonal matrix elements, we study the finite-size scaling of the standard deviation of this distribution. Once again, we find or scaling, Fig. 8(a), consistent with ETH scaling. The behaviour of the two operators differ only in the coefficient.

For SYK, one finds that most off-diagonal matrix elements vanish, Fig. 8(b). Since SYK is a free fermion model, the eigenstates may be organized by quasiparticle number. The operator written in the quasiparticle basis has number conserving terms and terms that change the quasiparticle number by two. For the operator, there are terms that change the quasiparticle number by zero, two and four. Therefore off-diagonal matrix elements in may be non-vanishing only when the eigenstates differ by two in the quasiparticle number which is an exponentially small number of possible off-diagonal matrix elements. Similar statements hold for any -Majorana operator for .

Because of this anomalously large weight at zero, the distribution of off-diagonal matrix elements has a delta function at zero, and so cannot be meaningfully plotted from finite-size data. However, for close to (close to the SYK point), the distribution retains signatures of this feature, as seen in the inset to Fig. 8(c). This is quantified by plotting the kurtosis of the distribution, which is 3 for a Gaussian distribution, as a function of , Fig. 8(c). (With the number of disorder realizations used, the kurtosis curves are not completely smooth, but the overall behavior is clear.) The highly peaked structure due to the many zeros at/near the SYK point leads to a sharp rise of the kurtosis. Comparison of the three sizes also shows that the region of proximity shrinks with increasing system size — for larger system sizes, a given departure from gaussian occurs for larger .

A similar feature was observed near integrability in Ref. Beugeling et al. (2015) for multiple non-random local models. This suggests that our explanation above for zero off-diagonal elements can be adapted to a wide class of integrable models.

## Vii Discussion

In this paper, we have extended tests of ETH to a class of random zero dimensional interacting models that may be tuned from chaotic to integrable. Our detailed numerical analysis shows that the chaotic model SYK obeys the standard scaling form of the eigenstate thermalization hypothesis. The key results are that the two and four point matrix elements in the eigenstate basis scale as leaving a sharp distribution around zero mean for large finite both at low energies and around the middle of the spectrum. In contrast, the scaling behaviour of the integrable SYK model matrix elements is algebraic in over the entire spectrum.

The fact that SYK obeys ETH indicates that the model has thermal correlators in the long time limit even within single eigenstates. SYK, in contrast, might be expected to thermalize to a generalized Gibbs ensemble averages determined by the extensive number of conserved quantities. It turns out that the canonical and GGE predictions coincide in free fermion models for single particle operators Ziraldo and Santoro (2013), such as our operator and the operator considered in Ref. Magán, 2016. It may be interesting to explore further the quench dynamics of non-quadratic operators in complex and Majorana SYK.

The off-diagonal matrix elements in the energy basis determine the way in which subsystems approach equilibrium. Since ETH is obeyed by SYK one may make statements about the long-time behavior of different correlation functions for example, as observed in Ref. Cotler et al., 2017 that the time-dependent two-point correlators tend to zero in the long time limit for (mod while a nonvanishing limit is possible for the (mod sequence.

The off-diagonal matrix elements also control the behaviour of thermal out-of-time-ordered correlation functions (OTOCs) such as at inverse temperature . Such correlation functions have been much studied recently in the context of quantum chaos and scrambling. At early times, the local operators mutually commute while the unitary evolution causes a spreading of information in sufficiently chaotic systems leading to an exponential increase in the OTOC - the exponent having the interpretation of a Lyapunov exponent in the semiclassical limit. One might expect such scrambling of information to be an integral part of thermalization. The fact that the Lyapunov exponent is finite in SYK and zero in SYK García-García et al. (2017) is intuitively consistent with these expectations.

SYK has the remarkable feature in the large and limits that the Lyapunov exponent corresponds to maximal chaos . In contrast, it has been observed in numerical computations of the OTOC that the finite behaviour of SYK has a Lyapunov exponent that is only weakly temperature dependent. The gulf between the latter result and the large limit cannot be bridged through the finite size scaling admitted by exact diagonalization.

Similarly, there is only a weak sense in which SYK for small system sizes is a particularly efficient thermalizer. In local many-body interacting systems, insofar as ETH has been tested, one finds similar behavior to SYK in non-integrable models while SYK behaves much like local integrable models. However, significant departures are expected at very low energies between the random D models considered here and local Hamiltonians. This is because the low energy eigenstates of SYK exhibit a volume law entanglement entropy Fu and Sachdev (2016); Liu et al. (2017) whereas local models have, instead, area law entanglement at low energies. By showing that ETH scaling holds down to the bottom of the spectrum we have uncovered a new aspect of this unusual low energy behavior. One might speculate that these features are present in generic fully connected and nonintegrable interacting systems.

We have explored ETH for parameters interpolating between SYK and SYK. Remarkably, the low energy behavior departs from the usual ETH scaling because of the RG relevance of SYK whereas the middle of the spectrum of the mixed Hamiltonian is governed by nonintegrability and ETH scaling persists.

###### Acknowledgements.
We acknowledge conversations with A. Bäcker, J. Bardarson, W. Beugeling, S. Bhattacharjee, R. Dantas, M. Heyl, I. Khaymovich, A. Lazarides, I. Mandal, R. Moessner, P. Surowka, and D. Trapin.

### References

1. M. Srednicki, Phys. Rev. E 50, 888 (1994).
2. J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
3. L. D’Alessio, Y. Kafri, A. Polkovnikov,  and M. Rigol, Advances in Physics 65, 239 (2016).
4. M. Rigol, V. Dunjko,  and M. Olshanii, Nature 452, 854 (2008).
5. G. Biroli, C. Kollath,  and A. M. Läuchli, Phys. Rev. Lett. 105, 250401 (2010).
6. M. Rigol and L. F. Santos, Phys. Rev. A 82, 011604 (2010).
7. L. F. Santos and M. Rigol, Phys. Rev. E 82, 031130 (2010).
8. C. Neuenhahn and F. Marquardt, Phys. Rev. E 85, 060101 (2012).
9. W. Beugeling, R. Moessner,  and M. Haque, Phys. Rev. E 89, 042112 (2014).
10. G. P. Brandino, A. De Luca, R. M. Konik,  and G. Mussardo, Phys. Rev. B 85, 214435 (2012).
11. T. N. Ikeda, Y. Watanabe,  and M. Ueda, Phys. Rev. E 87, 012125 (2013).
12. R. Steinigeweg, J. Herbrych,  and P. Prelovšek, Phys. Rev. E 87, 012118 (2013).
13. W. Beugeling, R. Moessner,  and M. Haque, Phys. Rev. E 91, 012144 (2015).
14. S. Sorg, L. Vidmar, L. Pollet,  and F. Heidrich-Meisner, Phys. Rev. A 90, 033606 (2014).
15. H. Kim, T. N. Ikeda,  and D. A. Huse, Phys. Rev. E 90, 052105 (2014).
16. R. Mondaini, K. R. Fratus, M. Srednicki,  and M. Rigol, Phys. Rev. E 93, 032104 (2016).
17. V. Alba, Phys. Rev. B 91, 155123 (2015).
18. R. Mondaini and M. Rigol, Phys. Rev. A 92, 041601 (2015).
19. N. Lashkari, A. Dymarsky,  and H. Liu, arXiv preprint arXiv:1610.00302  (2016).
20. S. Nandy, A. Sen, A. Das,  and A. Dhar, Physical Review B 94, 245131 (2016).
21. D. J. Luitz and Y. Bar Lev, Phys. Rev. Lett. 117, 170404 (2016).
22. A. Chandran, M. D. Schulz,  and F. J. Burnell, Phys. Rev. B 94, 235122 (2016).
23. R. Mondaini and M. Rigol, Phys. Rev. E 96, 012157 (2017).
24. A. Y. Kitaev, “A simple model of quantum holography,”  (2015), KITP Program: Entanglement in Strongly-Correlated Quantum Matter, talks on April 7, 2015 and May 27, 2015.
25. J. Maldacena and D. Stanford, Physical Review D 94, 106002 (2016).
26. J. Maldacena, S. H. Shenker,  and D. Stanford, Journal of High Energy Physics 2016, 106 (2016).
27. A. Eberlein, V. Kasper, S. Sachdev,  and J. Steinberg, arXiv preprint arXiv:1706.07803  (2017).
28. A. M. García-García, B. Loureiro, A. Romero-Bermúdez,  and M. Tezuka, ArXiv e-prints  (2017), arXiv:1707.02197 [hep-th] .
29. N. Hunter-Jones, J. Liu,  and Y. Zhou, arXiv preprint arXiv:1710.03012  (2017).
30. J. Sonner and M. Vielma, arXiv preprint arXiv:1707.08013  (2017).
31. J. M. Magán, Phys. Rev. Lett. 116, 030401 (2016).
32. C. Liu, X. Chen,  and L. Balents, ArXiv e-prints  (2017), arXiv:1709.06259 [cond-mat.str-el] .
33. S. Ziraldo and G. E. Santoro, Physical Review B 87, 064201 (2013).
34. M. Rigol, Phys. Rev. A 80, 053607 (2009).
35. M. Rigol, V. Dunjko, V. Yurovsky,  and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
36. R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015)https://doi.org/10.1146/annurev-conmatphys-031214-014726 .
37. Y.-Z. You, A. W. W. Ludwig,  and C. Xu, Phys. Rev. B 95, 115150 (2017).
38. A. M. García-García and J. J. M. Verbaarschot, Phys. Rev. D 94, 126010 (2016).
39. J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker, D. Stanford, A. Streicher,  and M. Tezuka, Journal of High Energy Physics 2017, 118 (2017).
40. C. W. Beenakker, Reviews of modern physics 69, 731 (1997).
41. V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
42. Y. Y. Atas, E. Bogomolny, O. Giraud,  and G. Roux, Phys. Rev. Lett. 110, 084101 (2013).
43. S.-K. Jian and H. Yao, Phys. Rev. Lett. 119, 206602 (2017).
44. L. Vidmar and M. Rigol, Journal of Statistical Mechanics: Theory and Experiment 2016, 064007 (2016).
45. I. Kourkoulou and J. Maldacena, arXiv preprint arXiv:1707.02325  (2017).
46. W. Fu and S. Sachdev, Physical Review B 94, 035135 (2016).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters