Spectral and thermodynamic properties of the Sachdev-Ye-Kitaev model

I Introduction

The insurmountable technical difficulties posed by the theoretical description of the many-body nuclear forces have led to many effective descriptions of nuclei to bypass the microscopic Hamiltonian. A crude assumption is to replace the nuclear Hamiltonian by a random matrix ensemble Wigner (1951); Dyson (1962a, b, c, d, 1972); Guhr et al. (1998) only constrained by global symmetries (the Wigner-Dyson ensembles). Surprisingly good agreement was found between spectral correlations of highly excited nuclei and the analytical predictions of random matrix theory for energy scales of the order of the mean level spacing. Despite of its success, this approximation has evident shortcomings. The nuclear-shell model suggests that nuclear interactions are well described by a mean-field potential plus a residual two-body interaction while in the random matrix approach higher many-body interactions are equally important. Moreover it was noticed that the spectral density associated to these high energy nuclear excitations did not follow the semi-circle law, the random matrix theory prediction, but it is better approximated by the Bethe formula Bethe (1936).

In response to these problems, a model of fermionic random $k-$body interactions of infinite range, the so called $k-$body embedded ensembles, was proposed more than forty years ago Bohigas and Flores (1971a, b); French and Wong (1970, 1971) as a more accurate stochastic description of nuclei. Although the interactions are random, the effective Hamiltonian is sparse and therefore deviations from the Wigner-Dyson ensembles were expected. Indeed numerical Bohigas and Flores (1971b) and later analytical results Mon and French (1975) show that, in line with the experimental data, the spectral density is Gaussian for sufficiently small $k$, instead of following the semi-circle law. By contrast, spectral correlations are still close to the random-matrix prediction Verbaarschot and Zirnbauer (1984) for sufficiently close eigenvalues. For more information on the model, especially in the context of nuclear physics and quantum chaos, we refer to Benet and Weidenmüller (2003); Gomez et al. (2011); Brody et al. (1981); Kota (2014); Kota et al. (2011).

Recently, similar models of fermions with $k-$body infinite-range interactions, called Sachdev-Ye-Kitaev models (SYK) Kitaev ; Maldacena and Stanford (2016); Polchinski and Rosenhaus (2016); Engelsöy et al. (2016); Almheiri and Polchinski (2015); Magán (2016); Danshita et al. (2016); Garcia-Alvarez et al. (2016); Bagrets et al. (2016); Sachdev (2015); You et al. (2016); Gross and Rosenhaus (2016), and originally introduced in the study of spin liquids Sachdev and Ye (1993), are being intensively investigated in a completely different context: holographic dualities in string theory Maldacena (1999). Based on the same pattern of conformal symmetry breaking, it has been speculated Kitaev ; Maldacena and Stanford (2016); Polchinski and Rosenhaus (2016); Engelsöy et al. (2016); Almheiri and Polchinski (2015); Jensen (2016); Cvetič and Papadimitriou (2016) that, in the infrared limit, the holographic dual of an Anti-deSitter (AdS) background in two bulk dimensions AdS2 is closely related to one of the variants of the SYK model, namely, a model of $N$ Majorana fermions Kitaev in zero spatial dimensions with random two body interactions of infinite range. Green’s functions Bagrets et al. (2016); Jevicki et al. (2016); Maldacena and Stanford (2016); Polchinski and Rosenhaus (2016), thermodynamic properties Sachdev (2015), such as the low temperature limit of the entropy, and also out of equilibrium features Maldacena and Stanford (2016) such as the exponential growth of certain out-of-time-ordered correlators are strikingly similar in both models. The latter, related to quantum corrections in the gravity dual Maldacena et al. (2015), is also a signature of quantum chaotic features. More interestingly, it is believed that the SYK model may describe the low energy limit of a higher dimensional gauge theory with a string theory dual still to be named. Very recent results Witten (2016) suggest that disorder is not strictly necessary for a gravity-dual interpretation.

Despite these advances, the description of many aspects of the SYK model dynamics still poses severe technical, both numerical and analytical, challenges. In closely related problems such as quantum chaos and disordered systems, the spectrum and level statistics provide a rather comprehensive description of the quantum dynamics without the need of the more expensive computation of eigenvectors. In the context of the SYK model, spectral correlations have so far been investigated in You et al. (2016), where level repulsion was found, typical of a disordered metal, though its strength changes with the number of particles $N$ modulo 8.

Here we aim to fill this gap by carrying out an extensive analysis of the spectral density, thermodynamic properties, and both short-range and long-range spectral correlations of the SYK model, with $N$ Majorana fermions.

Our main results are summarized as follows: we show analytically that in the $N\to\infty$ limit the fourth and sixth cumulant of the spectral density vanish which strongly suggests that it is Gaussian. However its tail at finite $N$, that controls the specific heat, is well approximated by the semi-circle law. Results from exact diagonalization, for up to $N=36$ Majorana fermions, are fully consistent with the analytical findings, including results for the entropy and the specific heat. Spectral correlations that test short range correlation as the level spacing distribution are in good agreement with the random matrix prediction. We find that, in agreement with You et al. (2016), the Bott periodicity of the Clifford algebra that governs the Majorana fermions labels the global symmetries of the model. However we have observed systematic deviations from the random matrix predictions, for sufficiently well separated eigenvalues, that suggest that the model is not ergodic for short times. The point of departure from the universal results of random matrix theory increases with $N$ which is a strong indication of the existence of a Thouless energy Altshuler et al. (1988); Braun and Montambaux (1995); Bertrand and García-García (2016) for the system.

This paper is organized as follows: in the next section we introduce the model and discuss its spectral density. The thermodynamical properties of the model are evaluated in section III. Spectral correlations are computed in section IV. We finish with concluding remarks and some ideas for future research in section V. Some technical details involving the calculation of the cumulants and the symmetry properties of the gamma matrices are worked out in two appendices.

II The spectral density

Kitaev recently introduced Kitaev a model of interacting fermions aimed to explore its potential as a gravity-dual. The Hamiltonian is given by,

where $\chi_{i}$ are Majorana fermions that verify

The fermions are coupled by Gaussian distributed random variables $J_{ijkl}$ with probability distribution,

We note that Eq. (2) is the defining relation of an Euclidean $N$-dimensional Clifford algebra. Many interesting features of the model are a direct consequence of Clifford algebra properties. For instance, the Bott periodicity of the Clifford algebra suggests that the global symmetries of the Majorana fermions, that to some extent control the spectral properties of the model, are sensitive to the arithmetic nature of $N$. We shall see that this is indeed the case when we study level statistics later in the paper. It will also be helpful for our first objective: to derive analytical results for the many-body spectral density.

We will follow the strategy of Mon and French Mon and French (1975) of evaluating moments of the spectral density. In this model, this is again facilitated by noticing that the Euclidean Clifford algebra in $N$ dimensions of the Majorana fermions Eq. (2) is shared by Euclidean Dirac $\gamma$ matrices. Therefore it is possible to employ the full machinery developed in that context to compute the trace of a large number of Majorana fermions, a key part in the calculation of energy moments. We leave the details of the calculation to appendix B. Here we just define the moments, sketch the main steps of the calculation, and give the final expression as a function of the number of particles $N$. Since the Gaussian disorder distribution is an even function, all odd moments will vanish. From now on we will focus only on the even ones:

where $p=1,2,3\ldots$, $\langle\ldots\rangle$ stands for spectral and ensemble average. The strategy to evaluate $M_{p}(N)$ is straightforward: we first perform the Gaussian average, equivalent to summing over all possible contractions according to Wick’s theorem, and then we evaluate each of these terms, involving the trace of products of $\gamma$ matrices, by using properties of $\gamma$ matrices in $N$ Euclidean dimensions.

Denoting the product of four $\gamma$ matrices by $\Gamma_{\alpha}$, we have that the moments are given by

The Gaussian average over the random couplings $J_{\alpha}$ of the Hamiltonian (1), denoted by $\langle\cdots\rangle$, is equal to the sum over all possible contractions. In the limit $N\gg 2p$ almost all $\Gamma_{\alpha}$ have no overlapping indices so that they commute. Because of

we find that in this case all $(2p-1)!!$ contractions give the same contribution resulting in the moments

These are the moments of a Gaussian distribution resulting in a Gaussian spectral density. We have evaluated the exact analytical result for $M_{4}$ and $M_{6}$. This requires the evaluation of diagrams that are subleading in $N$. For that purpose it is helpful to note that when we have common $\gamma$ matrices in $\Gamma_{\alpha}$ and $\Gamma_{\beta}$ they commute or anti-commute depending on the number of common $\gamma$ matrices. This results in large cancellations suppressing the contribution of intersecting diagrams. Following this procedure the first two non-trivial normalized cumulants, $\kappa_{4}$ and $\kappa_{6}$, are easily obtained as a function of $N$ from the moments $M_{2p}(N)$ (see Appendix B for details),

with large $N$ asymptotics $512\times 11/N^{2}$ where from now on we set $J=1$.

For higher moments the combinatorial problem becomes increasingly difficult and the final expressions are rather cumbersome. However these few cumulants already contain interesting information.

As we have seen above, for $N\to\infty$ the normalized cumulants vanish for orders $8p\ll N$. This is a distinctive feature of a Gaussian distribution. Therefore the average analytical spectral density converges (non-uniformly) to a Gaussian of zero average and variance equal to $6/N^{3}$.

We note that a Gaussian spectral density is expected for models with an entropy $S=Nf(E/N)$ in the large $N$ limit. The only requirement is that $f$ is a smooth function that has a maximum. Gaussian behaviour in the central part of the spectrum, assuming a maximum at $E=0$, results after expanding $f$ around the maximum.

In Fig. 1 we compare the analytical predictions Eqs. (8-9) of the normalized fourth and sixth cumulants with numerical results obtained by using exact diagonalization techniques. The agreement is excellent.

In Fig. 2 we depict the average spectral density for $N=34$, the largest size for which we can obtain numerically the full spectrum, with the analytical prediction, a Gaussian distribution with a variance that has been fitted to the data. Here the agreement is good but we observe clear deviations in the tail of the density. The reason for that discrepancy is that corrections to the Gaussian distribution, as described by the moments above, are still of order one for $N=34$. We were unable to compute analytically the leading $N$ corrections to the Gaussian density of states. However, in the next section, we carry out a detailed numerical analysis of the tail of the average spectral density.

III Thermodynamic properties in the low temperature limit

Part of the renewed interest in the SYK model stems from the fact that its low temperature properties are similar to those of a gravity background that in the infrared limit is well described by AdS2 geometry. Typical features includes a finite entropy at zero temperature, a ground state energy that is extensive in the number of particles, and a specific heat linear in temperature but with a prefactor different from that of free fermions. There are already approximate analytical predictions Maldacena and Stanford (2016); Jevicki et al. (2016) in the literature for these observables. Exact numerical diagonalization of the SYK Hamiltonian Eq. (1) was employed in Maldacena and Stanford (2016) to compute the zero temperature entropy Maldacena and Stanford (2016). We are not aware of exact diagonalisation results for the specific heat or the ground state energy. In this section we address this problem by a detailed numerical study of the tail of the spectrum that controls the thermodynamic properties in the low temperature limit. We start with the ground state energy. The lowest eigenvalue of the SYK Hamiltonian, $E_{\rm min}$, is the ground state energy of the SYK model with $N$ Majorana fermions Eq. (1). Due to the fermionic nature of model we expect $E_{\rm min}$ to be proportional to $N$. In Fig. 3 we show the ensemble average of $E_{\rm min}$ versus $N$ and it indeed shows a nice linear asymptotic dependence on the dimension $N$.

From a careful fitting of the numerical data we find that the tail of the spectrum is well approximated by

which also determines the low-temperature limit of the partition function,

The low temperature limit of the SYK model is given by Maldacena and Stanford (2016); Jevicki et al. (2016),

where the ground state energy, $E_{0}$, the entropy $S_{0}$ and the specific heat coefficient, $c$, are all proportional to $N$. The prefactor $\beta^{-3/2}$ is an order one contribution coming from one-loop quantum corrections and $c_{0}$ is a temperature independent constant. Comparing to Eq. (11) we can make the identification

In Fig. 4 we depict $\log a(N)$ and $b(N)$ by fitting of the exact partition function computed numerically by exact diagonalization. The zero temperature entropy and the ground state energy are then obtained from Eq.(13):

The value of $S_{0}$ is in rough agreement with the result $\sim 0.23N$ obtained by Maldacena and Stanford Maldacena and Stanford (2016).

We now move to the calculation of the specific heat. In the very low temperature limit with $\beta J\gg N$ we can expand the partition function as

It would be tempting to also make the identification

but in the parameter range we are looking at it is not justified to expand the exponential. Rather, we determine the specific heat coefficient $c$ by directly fitting the $\beta$-dependence of the specific heat,

where the internal energy per particle, $U(T)$, is defined in the usual way,

Setting $J=1$ for convenience, and using the low temperature expansion of the partition function given in Eq. (12),

where the exponent $q$ that controls the one-loop quantum correction $\beta^{q}$ to the partition function is left as a free parameter rather than fixing it to the perturbative Kitaev ; Maldacena and Stanford (2016) prediction $q=3/2$.

In terms of the eigenvalues $E_{k,p}$ of the $p$’th member of the ensemble of SYK Hamiltonians, the specific heat per particle is given by

For a given realization of the random Hamiltonian, the fluctuations of the average energy,

give rise to significant finite size contributions to the specific heat which can be eliminated by performing the ensemble average relative to the average energy for each realization of the SYK Hamiltonian, i.e.,

For a large number of particles this procedure should be equivalent to the calculation according to Eq. (20). However, for the values of $N$ we work with, this finite size effect must be removed in order to obtain accurate results for the low temperature limit of the specific heat.

The finite size effects discussed in the previous paragraph decrease rapidly with the total number of particles. As an example we show in Fig. 5 the temperature dependence of the specific heat for $N=28$ (left) and $N=36$ (right). We show both the result where the specific heat is calculated according to Eq. (20) (red dots) and the result where we first calculate the specific heat for each realization of the Hamiltonian and then perform the ensemble average as given in Eq. (24) (blue dots). The curves are fits to the blue dots.

Except for $N=36$, where we have only $2000$ eigenvalues for each configuration and use a linear fit on a shorter fitting interval, we use cubic fits

In Fig. 6 we show the $N$-dependence of $q(N)$ (left) and $c(N)$ (right) which are fitted by a constant for $N\geq 28$ (see curves). This results in the following estimates for the exponent $q$ in Eq.18 that controls one-loop quantum corrections and the specific heat coefficient

The value of $q$ is consistent with the estimate $q=3/2$ Maldacena and Stanford (2016) from an analytical calculation of one-loop quantum corrections to the classical action. It is also in agreement with the semicircular form of the spectral density, see Eq. (10). Likewise the analytical estimation of the specific heat coefficient $c/N=0.396$ Maldacena and Stanford (2016) is also consistent with our numerical results.

We note that all the results of this section are based on the ansatz Eq. (10) for the density of states. The exponent $1/2$ of the prefactor was chosen because it gave the best fit to the numerical results. However there is an indirect theoretical justification for that exponent. In the recent literature on the SYK model there are several studies Kitaev ; Maldacena and Stanford (2016); Polchinski and Rosenhaus (2016); Jevicki et al. (2016) of the one-point temporal correlation function which is the Fourier transform of the strength function

where $E_{0}$ is the $N$-particle ground state energy, $|k\rangle$ are eigenstates with $N\pm 1$ particles and $\gamma_{\alpha}$ is an Euclidean $\gamma$ matrix. These results are based on perturbative semi-classical techniques that typically are valid only up to time scales of the order of the Ehrenfest time. However in Bagrets et al. (2016) a non-perturbative treatment of quasi-zero modes enlarged the time domain of applicability of the analytical results to scales shorter but of the order the Heisenberg time. Interestingly, it was found Bagrets et al. (2016) that, in an energy representation, the strength function for low energies $\propto\sqrt{E-E_{0}}$. In principle the strength function is unrelated to the many-body spectral density Eq. (10) because the former provides also information of the correlations between eigenvalues and eigenvectors. However, if the eigenvectors and the eigenvalues are uncorrelated, as is the case for the Wigner-Dyson random matrix ensembles, the strength function is proportional to the spectral density. Below we will see spectral correlations of the SYK model are well described by the Wigner-Dyson ensembles which justifies a posteriori the ansatz Eq. (10) for the tail of the spectral density.

In summary, we have shown that the spectral density of the SYK model is Gaussian in the limit of a large number of particles $N$ so it is qualitative different from the semi-circle law typical of random matrices. However for a fixed finite $N$, the tail of the spectral density is close to a semi-circle law while the center is Gaussian. The value of the zero temperature entropy and specific heat coefficient, obtained numerically from the tail of the spectrum and the low-temperature behavior of the partition function, are close to previously obtained analytical estimates Maldacena and Stanford (2016); Jevicki et al. (2016).

IV Spectral correlations

In this section we investigate eigenvalue correlations that provide valuable information on the dynamics of the system. We focus on long time scales of the order of the Heisenberg time $\sim\hbar/\Delta$ where $\Delta$ is the mean level spacing. Disordered metals, or quantum chaotic systems, are expected to be described by the invariant random matrix ensembles in this region. Physically, agreement with random matrix theory predictions indicates that an initially localized wave packet reaches the boundary of the sample for sufficiently long time scales. For a disordered insulator we expect level correlations to be described by Poisson statistics. Although in the literature on $k-$body embedded fermionic ensembles there are some reports of Poisson statistics for two-body random interactions in the dilute limit Benet et al. (2001), there is broad evidence from numerical and analytical findings Srednicki (2002); Verbaarschot and Zirnbauer (1984); Bohigas and Flores (1971b) that level statistics are very close to the random matrix theory prediction at least for short-range eigenvalue correlations.

As was mentioned in the introduction, the only previous study of spectral correlations in the SYK model You et al. (2016) investigated numerically the ratio of consecutive level spacings which only explores time scales of the order of the Heisenberg time. For shorter time scales, corresponding to energy scales beyond the mean level spacing, level statistics for the SYK model is yet an open problem. We shall see that level statistics in this region are well described by random matrix theory though deviations, that decrease with $N$, are systematically observed for larger spectral distances corresponding to time scales much shorter than the Heisenberg time.

The universality class for the spectral correlations is determined by the anti-unitary and involutive symmetries of the system. Since the SYK Hamiltonian does not have any involutive symmetries, the universality class is given by the Wigner-Dyson random matrix ensembles with a Dyson index $\beta_{D}=1$, 2 or 4. The first case is when the anti-unitary symmetry squares to one, the second case when there are no anti-unitary symmetries, and the third case when the anti-unitary symmetry squares to -1. The SYK Hamiltonian has two anti-unitary symmetries (See Table I)

which is equivalent to one irreducible anti-unitary symmetry, $C_{1}K$, and the unitary symmetry $C_{1}KC_{2}K$. Physically, the symmetries $C_{1}K$ and $C_{2}K$ are charge conjugation symmetries which are equal to the product of the “even” gamma matrices or “odd” gamma matrices, respectively (choosing the right labeling for “even” and “odd”). Therefore, $C_{1}KC_{2}K\sim\Gamma_{5}$ with $\Gamma_{5}={\rm diag}(1,\cdots,1,-1,\cdots,-1)$ in a chiral representation of the Dirac $\gamma$ matrices. In this representation the SYK Hamiltonian splits into two diagonal block matrices of equal size. If $C_{1}KC_{2}K=\pm\Gamma_{5}$, the charge conjugation matrix commutates with the projection on the diagonal blocks. If $(C_{1}K)^{2}=1$ it is possible Porter (1965) to find an $H$-independent basis for which the blocks become real, corresponding to a Dyson index $\beta_{D}=1$. Moreover, if $(C_{1}K)^{2}=-1$, it is possible to construct an $H$-independent basis for which the Hamiltonian can be arranged into quaternion real matrix elements corresponding to a Dyson index $\beta_{D}=4$. If $C_{1}KC_{2}K=\pm i\Gamma_{5}$, the charge conjugation matrix does not commute with the projection onto the blocks. Therefore we cannot use these symmetries to construct a basis for which the Hamiltonian becomes real or quaternion real. Since there are no unitary symmetries the matrix elements of the SYK Hamiltonian are complex corresponding to a Dyson index $\beta_{D}=2$. However, the symmetry $C_{1}K$ still can be used to show that both blocks have the same eigenvalues (see Kieburg et al. (2015) for a similar reasoning). We refer to Appendix A for all technical details.

For our study we employ the level spacing distribution $P(s)$ (29), the probability to find two neighboring eigenvalues separated by a distance $s=(E_{i+1}-E_{i})/\Delta$, and the number variance $\Sigma^{2}(L)$ (31), that describes fluctuations in the number of eigenvalues in a spectral window of size $L$ again measured in units of the mean level spacing $\Delta$. The latter, a long-range spectral correlator directly related to the two-point correlation function, gives information on the quantum dynamics for times scales of the order but much larger than the mean level spacing (Heisenberg time). We shall use it to investigate deviations from random matrix predictions. The former is more suited to study longer time scales $\approx\hbar/\Delta$ and also provides indirect information on higher order correlation functions.

We investigate level statistics numerically by an exact diagonalization of the upper block of the Hamiltonian (1) for $N\leq 36$. The first step in the spectral analysis is the unfolding of the spectrum Guhr et al. (1998), namely, to rescale the spectrum so that the mean level spacing is the same for all energies. This is a necessary condition to compare level statistics in different parts of the spectrum. For that purpose, for each $N$, we employ the averaged smooth staircase function (the integral of the spectral density) resulting from a fifth order polynomial fitting involving only odd powers, to unfold the spectrum. The spectrum rescaled in that way, which has unit mean level spacing for all energies, is ready for the level statistics analysis. We have observed that level statistics are similar for all energies. Except for $N=36$, where we have only obtained about $2\%$ of eigenvalues close to the edge of the spectrum, we have taken about $70\%$ of the eigenvalues around $E\approx 0$.

The level spacing distribution $P(s)$ is the probability to find two eigenvalues separated at a distance $s$ in units of $\Delta$ with no other eigenvalues in between:

In an insulator it is given by Poisson statistics: $P(s)=e^{-s}.$ By contrast, the random matrix prediction, that applies to a disordered metal and to a quantum chaotic system, is very well approximated by the Wigner surmise,

Level repulsion, $P(s)\to 0$ for $s\to 0$, is a distinguishing feature of extended states though its strength depends on the global symmetries of the Hamiltonian (1). For systems that admit a real representation of the Hamiltonian, due to time reversal invariance (or more generally due to an anti-unitary symmetry that squares to 1), $\beta=1,\;a_{1}={\pi}/{2},\;b_{1}={\pi}/{4}$. Similarly if the Hamiltonian only admits a complex representation, due for instance to the breaking of time translational invariance as a consequence of a magnetic field or flux, $\beta=2,\;a_{2}={32}/{\pi^{2}},\;b_{2}={4}/{\pi}$. Finally the case $\beta=4,\;a_{4}={262144}/{729\pi^{3}},\;b_{4}={64}/{9\pi}$ corresponds to systems with time-reversal symmetry and strong spin-orbit interactions leading to a doubly degenerate spectrum (or more generally to systems with an anti-unitary symmetry that squares to $-1$). It is typical of random matrices with quaternionic entries.

In Fig. 7 we plot $P(s)$ for $N=28,\;N=34$ and $N=36$. Excellent agreement with the random matrix prediction is found in all cases. As can be seen from Table I, $N=28$ belongs to the Gaussian Symplectic Ensemble (GSE) universality class ($\beta_{D}=4$), while $N=34$ belongs to the Gaussian Unitary Ensemble (GUE) universality class ($\beta_{D}=2$). We note that the $N$ dependence of the universality class was already reported in You et al. (2016), although it was not discussed that this was a simple consequence of two features of Clifford algebras: the existence of real, complex or quaternionic representations for different values of the dimensionality $N$ and Bott periodicity, namely, these representations follow a periodic pattern, in this case the Bott periodicity is $N\,{\rm mod}8$. An example of a period is: $N=36$: GSE, $N=34$: GUE, $N=32$: GOE, $N=30$: GUE, and so on.

In Fig. 8 we depict $P(s)$ for $N=16$ and $N=32$ both belonging to the Gaussian Orthogonal Ensemble (GOE) universality class. Even though the large difference in size, we do not observe important differences between the two cases. We will see in the following analysis of the number variance, a long-range spectral correlator that deviations from random matrix theory eventually occur for larger eigenvalue separations which indicates that the SYK model is not ergodic for sufficiently short time scales.

The number variance is defined as the variance of the number of levels $N$ inside an energy interval that has (in units of the mean level spacing) $L$ eigenvalues on average:

For a Poisson distribution typical of an insulator, different parts of the spectrum are not correlated, so the number variance is linear with slope one, $\Sigma^{2}(L)=L$.

The random matrix prediction, that also occurs in non-interacting Altshuler et al. (1988); Braun and Montambaux (1995) and strongly coupled Bertrand and García-García (2016) disordered metals below the Thouless energy, is that level repulsion causes, for $L\gg 1$, a slow logarithmic increase, usually termed level or spectral rigidity of the number variance:

with $c_{1}={2}/{\pi^{2}},\;c_{2}=c_{1}/2,\;c_{4}=c_{1}/4$, $d_{1}=d_{2}=2,\;d_{4}=4$, $e_{1}=-\pi^{2}/8,\;e_{2}=0,\;e_{4}=\pi^{2}/8$ and $\gamma=0.5772\ldots$ is Euler’s constant. In Fig. 9 we depict the number variance for several values of the system size, $N=28$, $N=32$ and $N=34$, each of them belonging to a different universality class: GOE for $N=32$, GUE for $N=34$ and GSE for $N=28$. For all universality classes we find an excellent agreement with the random matrix prediction for small $L$. However we observe systematic deviations for sufficiently large $L\gtrsim 30$. As $N$ increases the region of agreement with random matrix increases as well, namely, deviations are observed only for larger $L$.

In Fig. 10 we depict the number variance for two sizes ($N=22$ and $N=34$) belonging to the same universality but one matrix size much smaller than the other. The idea is to study finite size effects, related to mesoscopic fluctuations in the number variance. For small $L\leq 20$ the number variance follows the GUE prediction for both sizes. However for larger $L$, deviations from the random matrix result occur much earlier, and grow much faster, for $N=22$ than for $N=34$. An eyeball estimate suggests that the region of agreement with random matrix predictions scales approximately as $2^{N/8}$.

Several conclusions can be drawn from these results: a) the SYK model has spectral correlations similar to that of a disordered metal or a quantum chaotic systems even for energy scales much larger than the inverse mean level spacing, b) deviations for sufficiently large scales, suggest that, unlike a dense random matrix, the SYK model is not ergodic for sufficiently short time scales. This is expected as the Hamiltonian is rather sparse with only $\sim N^{4}$ non zero elements. This feature is also required for a gravity-dual interpretation where it is expected that, for times of the order of the Ehrenfest time $\sim\log 1/\hbar$, certain correlation functions grow exponentially at a rate controlled by the Lyapunov exponent of the system Maldacena et al. (2015), c) the fact that, as $N$ increases, deviations from random matrix occur for larger $L$ is a strong indication that the observed chaotic features persist in the thermodynamic limit. It also suggests the existence of the equivalent of a Thouless energy in the system related to the typical time necessary to explore the full available phase space.

V Outlook and conclusions

We have shown analytically that, in the limit of large number of particles, the SYK Hamiltonian has a Gaussian spectral density, although for a fixed finite number of particles, we have found numerically the tail of the density is well approximated by the semicircle law. Level statistics are well described by random matrix theory up to energy scales much larger, but still of the order, of the mean level spacing. Deviations from random matrix theory for larger energies, or shorter times, are an indication that the model is not ergodic for short times. Together with previous results, this a further confirmation that the SYK model has quantum chaotic features at any time scale. According to Maldacena and Stanford (2016), this is an expected feature in field theories with a gravity-dual. Indeed, we have numerically calculated the specific heat and the entropy and found that the low temperature thermodynamic properties of the SYK model are similar to those of a gravity background with a AdS2 infrared limit. To some extent, our work on the SYK model shows that a compound nucleus may have a gravity dual. Finally we mention a few venues for further research. It would be interesting to explore metal-insulator transitions in the model by reducing the range of the interaction from infinity to a power-law decay. Another interesting problem is to evaluate analytically the two level correlation function in the $N\to\infty$ limit by the replica trick by following the procedure of Verbaarschot and Zirnbauer (1984) for the $k$-body embedded ensemble. Similarly, the analytical evaluation of the leading finite $N$ corrections of the spectral density, by a careful evaluation of higher order $N$ moments, would provide a full description of the low temperature thermodynamic properties of the model. This is necessary step for a full understanding of the relevance of the SYK model in holography. We plan to address some of these problems in future publications.

Appendix A Construction of the γ𝛾\gamma matrices

The $\gamma$ matrices are constructed iteratively starting from the $\gamma$ matrices in two dimensions

to extend it to $d+2=N$ dimensions where $N$ is the even number of Majorana fermions. As we will see below, in this representation, the product of four gamma matrices is block diagonal.

We can construct two anti-unitary symmetry operators (Note that the gamma matrices in $C_{1}$ are purely imaginary while the $\gamma$ matrices in $C_{2}$ are purely real.)

where $K$ is the complex conjugation operator (we could have interchanged the labels of $\gamma_{1}$ and $\gamma_{2}$ so that $C_{1}$ would have been the product of the odd gamma matrices and $C_{2}$ the product of the even gamma matrices). They satisfy the symmetry relations

with $\mu=1,\ldots N$. Since the Hamiltonian is a sum of products of four $\gamma$ matrices, we have

In the above table, which was also given in the main text, we give the main properties of these anti-unitary symmetries.

Because of (37) we have that $[\Gamma_{5},H]=0$, with $\Gamma_{5}=i^{-N/2}\prod_{i=1}^{N}\gamma_{i}$, so that $H$ splits into two block-diagonal matrices of the same size. If $C_{1}KC_{2}K=\pm\Gamma_{5}$, then

In this case we have that $(C_{1}K)^{2}=(C_{2}K)^{2}=\pm 1$. If $(C_{1}K)^{2}=1$ it is possible to find an $H$-independent basis in which $H$ becomes real, and the corresponding random matrix ensemble is the Gaussian Orthogonal Ensemble (GOE). If $(C_{1}K)^{2}=-1$ the Hamiltonian is self-dual quaternion up to an $H$ independent unitary transformation which corresponds to the Gaussian Symplectic Ensemble. In this case the eigenvalues of $H$ are a multiple of the quaternion identity and are thus doubly degenerate.

If $C_{1}KC_{2}K=\pm i\Gamma_{5}$ the projection operator is given by

but because of the “$i$” this projection operator does not commute with $C_{1}K$ or $C_{2}K$. So there are no anti-unitary symmetries when $H$ is block-diagonal, and we are in the universality class of the Gaussian Unitary Ensemble. In this case the charge conjugation matrices anti-commute with $\gamma_{5}$,

so that $C_{1}$ and $C_{2}$ are block off-diagonal

then the anti-unitary symmetries (37) result in the relation

Because $A$ and $B$ are Hermitian and $c_{i}^{*}c_{i}=-1$ we find from the secular equation that $A$ and $B$ have the same eigenvalues.

Appendix B Calculation of the fourth and sixth Cumulant

In this appendix we calculate the normalized fourth and sixth cumulant for the Hamiltonian of the SYK model.

The normalized fourth cumulant is given by

We now to proceed to the calculation of $M_{4}(N)$. The Gaussian average is the sum over all pairwise contractions. Because $\Gamma_{\alpha}^{2}=1$ with $\Gamma_{\alpha}$ a product of four different gamma matrices we find that the nested contractions are given by

with the factor 2 corresponding to the two contractions 4a) and 4b) in Fig. 11. For the intersecting contraction, see Fig. (11) (4c), we have to evaluate the trace

with $q$ the number of gamma matrices that $\alpha$ and $\beta$ have in common. For the sum over $\alpha$ and $\beta$ we thus obtain (see diagram 4c) in Fig. 11

Note that, as a check of this result, that without the factor $(-1)^{q}$ the sum over $q$ just gives ${N\choose 4}$. The result $T_{4c}$ can be simplified to

This results in the normalized fourth order cumulant

B.2 The sixth order cumulant

In this subsection we evaluate the normalized sixth order cumulant which in terms of the moments is given by

Since $M_{4}(N)$ was computed in the previous section we focus on $M_{6}(N)$. The Gaussian integral for the sixth moment is again evaluated by summing over all pairwise contractions. In this case there are fifteen diagrams, and five of them are nested, see Fig. 11 (6a-e). The nested diagrams are simply given by $M_{2}^{3}(N)$. The next simplest class of diagrams are those where two neighboring Hamiltonians are contracted, while the contractions of the remaining factors are intersecting, see Fig. (11)(f-k). Their contribution to the sixth moment is given by

By a cyclic permutation of the factors in $\textrm{Tr}H^{6}$, it is clear that the diagrams in Fig. 11 6l-n are the same. If we fix the index of the second factor in diagram 6l, it is clear that by commuting the factors as

we obtain the same combinatorial factor for the sum over $\alpha$ and $\gamma$ as in diagram 4c. We thus find

The most complicated diagram is diagram 6o corresponding to the trace

The simplest way to do combinatorics is to think of $\Gamma_{\beta}\Gamma_{\gamma}$ as a product of 8 gamma matrices with $q$ gamma matrices in common while $\Gamma_{\alpha}$ share $l$ gamma matrices with $\Gamma_{\beta}\Gamma_{\gamma}$ and of those $l$ there are $l-m$ in the common factors. The result for this diagram is given by

Again, as a check of this result, if the phase factor $(-1)^{q+m}$ is put to one, we find $M_{2}^{3}(d)$.

Combining all contributions we find the normalized sixth cumulant

References