Exact moments of the Sachdev-Ye-Kitaev model up to order $1/N^2$

We analytically evaluate the moments of the spectral density of the $q$-body Sachdev-Ye-Kitaev (SYK) model, and obtain order $1/N^2$ corrections for all moments, where $N$ is the total number of Majorana fermions. To order $1/N$, moments are given by those of the weight function of the Q-Hermite polynomials. Representing Wick contractions by rooted chord diagrams, we show that the $1/N^2$ correction for each chord diagram is proportional to the number of triangular loops of the corresponding intersection graph, with an extra grading factor when $q$ is odd. Therefore the problem of finding $1/N^2$ corrections is mapped to a triangle counting problem. Since the total number of triangles is a purely graph-theoretic property, we can compute them for the $q=1$ and $q=2$ SYK models, where the exact moments can be obtained analytically using other methods, and therefore we have solved the moment problem for any $q$ to $1/N^2$ accuracy. The moments are then used to obtain the spectral density of the SYK model to order $1/N^2$. We also obtain an exact analytical result for all contraction diagrams contributing to the moments, which can be evaluated up to eighth order. This shows that the Q-Hermite approximation is accurate even for small values of $N$.


Introduction
Although the study of strongly interacting quantum many body systems has a long history, many aspects still remain poorly understood. One of the difficulties is that the size of the Hilbert space increases exponentially with the number of particles which severely limits the scope of numerical studies. This is why analytical studies of even simplified many-body systems contribute significantly to our understanding of this problem. One such model is the Sachdev-Ye-Kitaev (SYK) model [1][2][3] which is a Hamiltonian system with an infinite range q-body interaction acting on a many-body Hilbert space of Majorana fermions. A similar model with complex fermions was introduced several decades ago in the context of nuclear physics where it became known as the two-body random ensemble [4][5][6][7][8][9]. The motivation for this model is that the nuclear interaction is mostly a two-body interaction with matrix elements that appear close to random. It was also known that the level spacing distribution of nuclear levels can be described by Random Matrix Theory [10][11][12][13][14] , while the overall shape of the nuclear level density does not resemble a semi-circle at all but increases exponentially as exp(c √ E) where E is the energy above the ground state [15]. The two-body random ensemble addressed both of these issues, and has been studied intensively since then [6,[16][17][18][19][20][21][22].
In previous works [42,43], two of us have studied both the thermodynamic and spectral properties of the SYK model for q > 2, and have clearly established that the short-range spectral correlations are given by random matrix theory which is a necessary ingredient for the model to be quantum chaotic and therefore to have a gravity dual with black hole solutions. Moreover, it was found, by an explicit analytical evaluation of the moments of the spectral density, that it grows exponentially for low energies, a typical feature of conformal field theories [44] and therefore of gravity backgrounds with a field theory dual [45,46]. One of the surprising results of these works is that the spectral density at finite N , even for relatively small N , is very close to the weight function of the Q-Hermite polynomials.
Rigorous results for the moments of a similar random spin model [47], and very recently for the SYK model itself [48], show that in the large N limit its spectral density converges to the weight function of the Q-Hermite polynomials only for q ∝ √ N while in Refs. [42,43] q was fixed and N was relatively small, so such a good agreement was not expected. Another surprising feature of the Q-Hermite approximation is that for low temperatures it reproduces exactly the SYK partition function which in this limit reduces to the Schwarzian action and it is 1/N exact [3,26,29,30]. This is again rather unexpected because this region is in principle controlled by high moments of order N/q 2 where deviations from the Q-Hermite result should be larger.
The main goal of the present paper is to study why the Q-Hermite approximation is so accurate. This question is addressed in two ways. First, by an analytical computation of the 1/N 2 corrections to all moments which enables us to obtain the density up to that order for any q , and second, by an exact analytical calculation of the finite-N moments up to order eight. We note that originally the SYK model was only formulated for even q. However, it also makes sense for odd q, when the Gaussian-distributed operator becomes the supercharge of a supersymmetric Hamiltonian [49][50][51][52][53]. Unless stated otherwise, our results are valid for both even and odd q, and for odd q they refer to the spectral properties of the supercharge.
We proceed by using the moment method for the spectral density. The 1/N 2 corrections are derived in two steps. First, we show that the 1/N 2 correction to the Q-Hermite result for each contraction diagram is proportional to the total number of triangular loops of the corresponding intersection graph. In the second step, we evaluate the sum over all diagrams. This is a graph-theoretic problem with combinatorial factors that can be determined from the exact expressions for the moments for q = 1 and q = 2. The moments can be summed into a 1/N 2 correction to the spectral density.
Finally, we note that other aspects of 1/N expansions in the SYK model were discussed in [54][55][56] but they do not overlap with the present work. This paper is organized as follows: In section 2, we define the SYK model, discuss the moment method and introduce the graphical representations for the calculation of the moments. In section 3 the 1/N expansion will be discussed, where we also obtain the 1/N 2 correction for a given diagram. In section 4 we obtain a general and exact formula to evaluate the contraction diagrams. In section 5 we derive triangle counting formulas, which in turn give us the total 1/N 2 correction to moments. After obtaining the 1/N 2 -exact moments, the corresponding correction to the spectral density is evaluated in section 6. In section 7 we compute the exact sixth and eighth moments to further clarify the properties of the approximations we made. In section 8 we comment on the nature of the obtained results. Concluding remarks and prospects for future work are discussed in section 9.
where the Γ α are defined in terms of 2 N/2 × 2 N/2 dimensional Euclidean gamma matrices as with anti-commutation relations 1 The subscript α represents an index set with q elements: α = {i 1 , i 2 , . . . , i q }, with 1 ≤ i 1 < i 2 < · · · < i q ≤ N . Hence α can have N q different configurations. The couplings J α are random variables distributed according to where J is a dimensionful parameter that sets the scale. Note we have included the factor (i) q(q−1)/2 to make the Γ α Hermitian also for odd q, in which case H(J α ) is interpreted as the supercharge of the so called supersymmetric SYK model [49].

Moments and Wick contractions
An object of central interest is the spectral density ρ(E): where · · · denotes the ensemble average over the Gaussian distribution of J α . After a Fourier transform of the δ-functions, we can write Hence we can equivalently study the moments TrH k . Due to the J α → −J α symmetry of the ensemble, all odd moments must vanish, and thus It will be convenient to factor out the dimensionality of the Hilbert space and study the 2p-th moment defined by and normalize the moments with respect to the second moment It is easy to show Since the average over the J α 's is a Gaussian integration, the 2p-th moment is given by the sume of all possible (2p − 1)!! Wick contractions among p pairs of Γ's. A Wick contraction of the form, say, where the Einstein summation convention is assumed, can be represented by diagram (a) in figure 1. We will call diagrams like figure 1 contraction diagrams, which can be equivalently drawn as rooted chord diagrams on a circle [57], and we will use the two terms interchangeably for such diagrams in this paper.
The matrices Γ α 's satisfy where there is no summation over repeated indices in the first equality, and c αβ = |α ∩ β| is the number of common elements between sets α and β. We can use equation (2.12) to calculate traces of products of Γ α k like in (2.11) by permuting Γ α k 's until every two Γ α k 's with the same subscript neighbor each other. For intersecting neighboring contractions we We thus see that commuting two operators gives rise to the suppression factor which will play an essential role in the calculations of this paper. Using these relations, the trace in the example (2.11) can be written as Generically, a contraction with p contraction lines can be written as where n c is the number of crossings in the contraction diagram, α i k , α j k belong to {α 1 , . . . , α p } and they label the contraction lines that cross each other.

Intersection graphs
The chord diagrams that contribute to the 2V -th moment all have V contraction lines/chords. An intersection graph for a chord diagram with V chords is defined as follows: • Represent each chord by a vertex.
• Connect two vertices with an edge if and only if there is a crossing between the two chords that these two vertices represent.
We denote by G a generic intersection graph, by V the number of vertices and by E the number of edges of an intersection graph. Therefore, in the notation of the previous section, V = p and E = n c . We give some examples of such diagrams in figure 2.
Motivated by the combinatorial factors that enter in the scaled moments (2.9), we define the following object associated with each contraction diagram and hence with each intersection graph, contributing to the scaled 2V -th moment, where c(G) = E k=1 c α i k α j k , the α i k α j k are all the edges in G and c α i k α j k = |α i k ∩ α j k |. It is clear that an intersection graph G completely determines the value of η G . With the above definitions, we have by Wick's theorem, where G i 's are the intersection graphs with p vertices.
Notice that in the language of intersection graphs, η defined in eq. (2.14) corresponds to a V = 2, E = 1 graph, that is, a single edge connecting two vertices.

Q-Hermite approximation
Generally, a contraction diagram cannot be reduced to an expression only involving η as we did in (2.13). The reason is that indices in more complicated contraction patterns cannot be treated as being independent. However, we obtain an important approximation if we nevertheless treat all crossings as independent: if a diagram has E crossings the result is then simply given by [43] This approximation expresses η G of an intersection graph G by a product of its edges. This approximation is in fact at least 1/N -exact, as will be discussed in detail in section 3 and appendix A. It is exact for chord diagrams where multiple crossings are indeed independent (having tree graphs as intersection graphs, see appendix D). The approximation (2.19) allows us to use the Riordan-Touchard formula [58,59] to sum over all intersection graphs [27,42]: where E i denotes the number of edges of G i . The unique spectral density that gives the moments (2.20) is the weight function of Q-Hermite polynomials [60]: with c N a normalization constant and E 0 a scale factor that drops out of the ratio (2.9). For this reason we refer to this approximation as the Q-Hermite approximation and also introduce the Q-Hermite moments We take equations (2.19), (2.20) and (2.21) as the first approximation to η G , the moments and the spectral density respectively, and we use this as the starting point to investigate further 1/N 2 corrections. We stress again that Q-Hermite approximation already contains many higher order terms in 1/N , although the approximation is only exact to order 1/N .

1/N expansion
The goal of this paper is to understand better why the Q-Hermite result discussed in the previous section is such a good approximation to the spectral density of the SYK model. This section is a step in this direction as we show that indeed there are no 1/N corrections to the Q-Hermite moments, and give an argument that the 1/N 2 corrections are determined by the total number of triangles in an intersection graph. In appendix A, we rigorously demonstrate this statement.
The scaled Q-Hermite moments M QH 2p /M p 2 only depend on η which has the 1/N expansion Keeping only the leading power in q at each order of 1/N , it can be shown that this simplifies to (see appendix B) The Q-Hermite moments thus have a nontrivial large N limit when q 2 /N is kept fixed. For q √ N we have η → 0 so that only the nested contractions contribute, which give the moments of a semi-circle. For q √ N we have to distinguish even and odd q. For even q we have η → 1 so that all contractions contribute equally which gives the moments of a Gaussian distribution while in the case of odd q we obtain η → −1 corresponding to the moments of the sum of two delta functions located symmetrically about zero. In the latter case, all scaled moments are equal to one (see eq. (C.9)).
To understand the corrections to the Q-Hermite moments we evaluate 1/N corrections at fixed q and p, The Q-Hermite result is obtained when all crossings between contraction lines are treated independently with each crossing contributing a factor η. Corrections of order 1/N occur when two crossed contraction lines have one common index. Since this involves a single crossing, this correction is the same for the exact result and the Q-Hermite result and we thus have that Corrections to the Q-Hermite result occur when the crossed contracting lines can no longer be permuted independently. Generally, this happens when intersection graphs have closed loops, and all vertices in the closed loop have at least pairwise common indices. If a pair of vertices does not have any common indices, they can be commuted or anti-commuted resulting in a loop that is no longer closed and is thus given by the Q-Hermite results. A closed loop of length k, thus differs by O(1/N k−1 ) from the Q-Hermite result. Therefore, for the O(1/N 2 ) correction we only need to consider the triangular closed loops. For a closed loop of three crossed contraction lines, say α, β and γ, let us consider the crossed pair βγ with one common index and let the crossed pair αβ also have a common index. This is a 1/N correction that is part of the Q-Hermite result, and thus contributes as q 2 /N to leading order in 1/N . Deviations from the Q-Hermite result to this closed loop occur when also the crossed pair αγ has a common index. Choosing this index of γ to be one from α gives a second factor q/N . The index can either be among the indices shared with β or not. We thus conclude that the 1/N 2 corrections due to a triangular diagram occur as The proportionality factor in (3.5) can be obtained from the simplest triangular intersection graph, whose value is referred to as T 6 , see figure 1 (c). This graph first occurs in the calculation of the sixth moment and can be calculated by keeping track of the combinatorial factors [42] (see section 7 for more details). From the large N expansion of T 6 and η we then find Therefore, this correction is 1/ √ N suppressed in the large-N limit with fixed q 2 /N . Since for the lowest non vanishing order, the triangles in the intersection graphs contribute independently, we arrive at the first main result of this paper: where T is the total number of triangles that occur in the intersection graph η G . This tells us the total 1/N 2 correction is obtained by counting the total number of triangles in all intersection graphs. In appendix A we will prove (3.7) by a calculation of the 1/N 2 corrections starting from an exact formula for all Wick contractions which will be derived in section 4.
In the double scaling limit, the Q-Hermite result for the moments only depends on q 2 /N . Therefore in this case the 1/N expansion is really in terms of this quantity only. Corrections to a term of order (q 2 /N ) k occur when we fix additional indices in a closed loop. Choosing the remaining indices gives a combinatorial factor of the form (with m an integer satisfying m N ) For completeness we also give the first three terms of the 1/N expansion of the Q-Hermite approximation of η G which follows from the 1/N expansion of η. For a contraction diagram with E crossings the we find

Exact result for the contraction diagrams
In this section we derive an exact analytical expression for all contraction diagrams contributing to the moments of the SYK model. We will use this result to prove (3.7) by an explicit calculation of the 1/N 2 corrections, which we defer to appendix A because this is a tedious calculation. The results of this section can also be used to obtain exact analytical results for low order moments and some examples are worked out in appendix E.
Since the phase factor c(G) that appears in the definition of η G (see eq. (2.17)) is dependent on the number of common elements in the index sets, it is natural to write the combinatorics also in terms of intersections of sets. Although c(G) is determined by intersections of pairs, the combinatorics will depend on intersections of arbitrary number of sets, so we introduce the objects Red region has cardinality c α 2 α 3 Red region has cardinality d α 2 α 3 which are the number of common indices in the vertices {α 1 , · · · , α l }, and which are the number of common indices in the vertices {α 1 , · · · , α l } that are not shared with any of the other α k . By convention α i and α j cannot label the same vertex if i = j.
Essentially, the c α 1 ···α k and d α 1 ···α k are the cardinalities of certain regions in the Venn diagram of {α 1 , · · · , α l }. Figure 3 illustrates the difference between c α 1 ···α k and d α 1 ···α k in the case of three index sets which occur in the calculation of the sixth moment. By the inclusion-exclusion principle, the two objects are related by and conversely, where stars in the subscripts indicate sums over the remaining indices, e.g. 5) and the same definition goes for the d α 1 ···αp * ··· * . Eq. (4.4) implies that c(G) can be written in terms of d's, hence we can write η G as The multiplicity factor M is the number of configurations that the index sets {α 1 , . . . , α V } can take given the values of the d α 1 ···α k . In general, a Venn diagram of V index sets is partitioned into 2 V regions by the boundaries of the index sets, hence the multiplicity factor is the number of ways to distribute N elements into 2 V regions, each region with its own cardinality. If the cardinality of each region is given by m i , then the multiplicity is given by the multinomial factor As an example, figure 4 explicitly shows the cardinalities of all eight regions partitioned by the boundaries of three index sets.
Our final result for the contribution of a given contraction diagram is thus given by, where and so on. The expression that appears in the denominator of the first factor for the multiplicity, is the cardinality of the union of all index sets, i.e. of all the circles in a Venn diagram like figure 4. One way to see this is by noticing that the indices of the sets with cardinality d k are shared by k of the α k and that k − 1 of them are not new. The expression (4.8) is the second main result of this paper. The general expression (4.8) is a sum over 2 V − V − 1 variables, which limits its practical applicability. However, it can be simplified in several cases. First, nested contractions, which correspond to isolated vertices in an intersection graph, just contribute a multiplicative factor of 1. So we do not need to include the sums over isolated vertices in (4.8), and this reduces the number of sums for these diagrams. A second simplification occurs if parts of an intersection graph are only connected by a single vertex (see appendix D). In that case the graph factorizes and each part can be evaluated separately by the formula (4.8).
The 1/N corrections are controlled by the term in (4.8). Hence we have arrived at a convenient starting point for large N expansions.

1/N 2 corrections to the moments
In section 3 we have seen that the 1/N 2 corrections for each graph are given by the number of triangles in an intersection graph. To obtain the 1/N 2 corrections to the moments, we  have to find the total number triangles in all intersection graphs contributing to the moment of a given order. If T i is the number of triangles in an intersection graph G i with E i edges, we have to evaluate In table 5 we give the numerical results up to 2p = 18. Both for even q and odd q strikingly simple patterns emerge: where T i and E i are the numbers of triangles and edges of the i-th intersection graph G i . These identities, which are one of the main results of this paper, will be proved in the second part of this section. To order 1/N 2 , the moments are thus given by for even q and by for odd q.
The proof of (5.2) and (5.3) is based on the following simple idea: on the one hand, the counting of the number of triangles occurring in intersection graphs contributing to the 2p-th moment is a graph-theoretic quantity, which is independent of the SYK model parameter q (except for the parity of q); on the other hand, the SYK model for q = 1 and q = 2 is exactly solvable, which means that for q = 1 and q = 2 we can obtain (M 2p − M QH 2p )/M p 2 to order 1/N 2 from the exact solutions. Then the proof follows by matching both sides of equation We first consider the simpler case q = 1, where all moments are known analytically [61].
we have and Since J α is Gaussian distributed, Tr(H 2p ) can be easily calculated. In fact we can easily recognize it as the p-th moment of χ 2 distribution with N degrees of freedom and the result is standard: For q = 1, the 1/N expansion of the Q-Hermite moments simplifies to The total 1/N 2 term is thus given by The second derivative can be calculated analytically (see appendix C) and is given by Subtracting this from the exact q = 1 result, eq. (5.9) we find This proves (5.3).
For q = 2 there is no compact formula for the 2p-th moment, however we can still compute the moments to 1/N 2 from the exact joint probability distribution of the coupling matrices. The computation is more involved and we give the derivation in appendix E, here we only quote the final result: Meanwhile the Q-Hermite result is given by Again the 1/N 2 sum can be computed using the techniques explained in appendix C: We finally find This proves (5.2).

Corrections to the spectral density
The moments, both for even q and odd q, satisfy Carleman's condition and hence uniquely determine the spectral density [62]. In the following two subsections we give the spectral density corrections for the even q and the odd q cases.

Spectral density for even q
We decompose the spectral density into the weight function ρ QH (E) of the Q-Hermite polynomials, determined by the Q-Hermite moments, plus a correction δρ(E) determined by the 1/N 2 correction to the moments, where ρ QH (E) is given by [43,60] and for even q we have 2 We normalize ρ QH (E) by This results in the normalization constant [60] After performing a Poisson resummation and ignoring certain exponentially small (in N ) terms [43], the spectral density away from |E| = |E 0 | simplifies to From this we deduce that at large N , It is simple to verify that the correction term gives the moments (5.4) consistent with the normalization of the ρ QH (E).
For any fixed value of energy E, E/E 0 is small since E 0 ∼ N , and the leading behavior of ρ QH is given by while the leading behavior of δρ is This is indeed a small correction in the point-wise sense both for large N at fixed q and in the double scaling limit in terms of the variable E/E 0 1. Unfortunately, this is not a small correction in the uniform sense, for example, if instead of a fixed E one looks at a fixed value of the scaling variable x = E/E 0 , then for x close to 1, the correction term δρ(E) ∼ 2 N/2 e −N/q 2 (6.11) becomes exponentially larger than the leading term (6.6) This prevents us from obtaining a meaningful correction to the free energy by integrating δρ. This is indeed consistent with the fact that the SYK partition function is 1/N exact in the low temperature limit [3,29,30,43] where it is dominated by the spectral density for E ≈ E 0 , so 1/N 2 corrections in this region must be spurious.

Spectral density for odd q
For odd q, η < 0, but the expression (6.2) for the Q-Hermite spectral density is still applicable. Following the steps of the even q calculation, it is straightforward to show that for large N and away from the edge of the spectrum, the spectral density (6.6) is given by [63] ρ QH (E) = c N cosh π arcsin(E/E 0 ) log |η| exp 2 arcsin 2 (E/E 0 ) log |η| . (6.13) The normalization constant can be determined from dEρ QH (E)dE = 2 N/2 . (6.14) In the large N limit, the integral can be evaluated by a saddle point approximation. Using that log |η| ∼ −2q 2 /N in this limit we find which gives exactly the leading order 1/q 2 correction to the zero temperature entropy [49]. In terms of units where the second moment is normalized to one, the correction to the spectral density with moments given by (5.5) is equal to In terms of physical units with M 2 = σ 2 this can be written as So also for odd q we find that the 1/N 2 correction to the spectral density is given by derivatives of its large N limit.
Correction terms in the form of δ-functions are not strange to random matrix theory: for example the 1/N correction of the Wigner-Dyson ensemble is proportional to δ-functions at the edges of the semi-circle [16].

Exact calculation of the sixth and eighth moment
In this section we give exact results for the sixth and eighth moment. More details can be found in appendix F.
The sixth moment was already calculated in [42]. All diagrams except for the rightmost contraction diagram in figure 1 coincide with the Q-Hermite result which allows the application of Riordan-Touchard formula For the sixth moment the Q-Hermite result is given by while the exact sixth moment is given by with In  The exact result for the 8th moment is given by It involves three new structures (see the intersection graphs in Table 3). They are still simple enough that they can be expressed as simple sums by inspection. However, they can also be derived starting from the general formula (4.8) and we give two examples in appendix F.1. The first structure corresponding to the square intersection graph (see table 3) is equal to The second structure corresponding to the square intersection diagram with one diagonal (see table 3) only differs by an additional phase factor The most complicated diagram is the one with 6 crossings corresponding to the rightmost graph in table 3. It is given by r .  The results for the sixth and eighth moment have been simplified using the convolution property of binomial factors. For example, for T 8 we initially obtain the result in the form of an 8-fold sum. The general result gives an 11-fold sum which can be reduced to this result as is worked out in detail in appendix F where we also list all contraction diagrams contributing to M 8 .
The results for the moments are also valid for odd q, even for q = 1. We have checked that the above expressions simplify to the q = 1 result in eq. (5.9) and are in agreement with moments obtained numerically from the exact diagonalization of the SYK Hamiltonian.
In figure 5 and figure 6 we show the N dependence of the sixth and eighth moment, respectively. We compare the exact result the Q-Hermite result to q = 1, 2, · · · , 8 and find that the two are close in particular for even q, even for small values of N .

The nature of the Q-Hermite approximation and higher order corrections
The results we have obtained are, at least superficially, contradictory. The 1/N 2 correction to the moments of the Q-Hermite approximation, one of the main result of the paper, does not improve the spectral density in a uniform sense. At the same time the exact computation of low order moments (up to 8th moment), show that the Q-Hermite approximation gives surprisingly accurate results even for relatively small N . Adding to this point, the 1/N 2 correction we have obtained in this paper only improves the Q-Hermite for 2p N/q 2 . At small N , it actually makes the approximation worse than using Q-Hermite results alone even for relatively small p.
All the above suggests that the Q-Hermite approximation should be understood as a re-summed finite N result, which happens to be 1/N -exact in the large N expansion, while the extra 1/N 2 corrections from triangle counting are strictly asymptotic. However, a resummation that is only 1/N exact does not imply that the Q-Hermite approximation must be so accurate at finite N . To achieve further understanding let us expand the Q-Hermite moments in powers of 1/N , where we have used the expansion of η and the identities given in appendix C to perform the sums i (−1) qE i E i and i (−1) qE i E 2 i . Comparing this to the difference and we observe that the Q-Hermite approximation is not only exact at order 1/N , but also almost exact at order 1/N 2 in the following sense: • The Q-Hermite approximation captures the term q 4 /N 2 , which is the only term of order 1/N 2 that survives in the scaling limit N → ∞ with q 2 /N = constant. We have already seen that this property holds to all orders in q 2 /N .
• Even for fixed q, the Q-Hermite approximation contains the dominant contribution to moments, that is, the leading coefficients of 1/N 2 in (8.1) and (8.2) go as p 4 q 4 , while the corrections in (8.3) and (8.4) go as p 3 q 3 , and we have p 4 q 4 p 3 q 3 already for relatively small p and q.
The second property is also likely to hold to all orders in 1/N .
Note that in the truncated form of the Q-Hermite result (8.1), the corresponding spectral density at each order is a Gaussian or a derivative of Gaussian with the same distribution width, so is the spectral density of the extra correction (8.3). Therefore the breakdown of the 1/N expansion for x = E/E 0 ≈ 1 discussed at the end of section 6.1 is really an artefact of the asymptotic nature of this expansion.

Conclusions and Outlook
We have obtained analytically an exact expression for the moments of the q-body SYK model. For any q, we have computed them explicitly up to 1/N 2 order. One surprising result of the calculation of the 1/N 2 order is that it allows a simple and beautiful geometric interpretation in the form of triangular loops of the intersection graphs. The 1/N corrections which are part of the Q-Hermite approximation are also geometric in nature and are given by the contribution from the edges. Our results can be generalized to higher orders in 1/N . Preliminary results for the 1/N 3 correction to the moments indicate they are also characterized by the geometry of the intersection graphs in a simple manner. From a more mathematical perspective, they generate a remarkable set of graph-theoretic identities. In particular, the 1/N 3 computation enumerates a particular linear combination of the last three geometric objects in table 3. It is one of the miracles of the SYK model that the 1/N 2 and even higher order corrections can be calculated analytically. On top of this, given that we have an exact expression for all contraction diagrams contributing to the moments, this might be an indication that the moment problem of the SYK model is completely solvable. In future work we hope to elaborate on this question.
The original motivation to carry out the moments calculation was to obtain a more accurate description of the spectral density and also, closely related, to understand better why the Q-Hermite approximation, which reproduces only the 1/N correction to the exact moments of the SYK model, is so close to the numerical SYK spectral density at least for even q. Even more remarkable is that the Q-Hermite spectral density exactly reproduces the temperature dependence of the free energy of the SYK model to leading order in 1/q 2 at all temperatures both for even and odd q. For the moment, we have at best partial answers to these questions.
We have found that 1/N 2 corrections to the moments, though exact, lead to a correction to the Q-Hermite spectral density which is only accurate sufficiently far away from the tail of the spectrum. As the spectral edge is approached, it gives unphysical results. This is indeed consistent with the fact that, close to the edge of the spectrum, the density is 1/N exact [29,30,64], so in this spectral region, the 1/N 2 correction must be an artefact of the asymptotic expansion in 1/N . In order to understand the reason for this unphysical behavior we first note that the in large N limit at fixed q, η → (−1) q and the Q-Hermite spectral density, the leading term of the expansion, tends to a Gaussian for even q and to δ-functions for odd q. Interestingly, the 1/N 2 correction to the spectral density can be expressed in terms of derivatives of this large N limit of the spectral density which strongly suggests that, in terms of a nonlinear σ-model for the spectral density, the 1/N corrections are an expansion about the trivial saddle point. We expect that the spectral edge of the Q-Hermite result is given by a non-trivial saddle point of this effective σ-model. Indeed a similar effect is observed in the calculation of 1/N corrections to the semi-circle law in random matrix theory [16]. It is not clear to us how to re-sum the asymptotic 1/N expansion so that it yields a vanishing 1/N 2 correction to the density close to the edge of the spectrum. In fact, we would need an expansion to all orders in 1/N to do that. An issue that further complicates the solution of this problem is the non-commutativity of the large p (order of the moment) and the large N limit. The main contribution for moments of order p N/q 2 comes from the edge of the spectrum, while when we first take the large N limit the main contribution to the moments resides in the bulk of the spectrum. The 1/N 2 corrections also share this nonuniform large N behavior.
Regarding the reason behind the unexpected close agreement between the Q-Hermite density and the exact spectral density of the SYK model, we have found that the leading 1/N 2 contribution in q to the 2p-th moment, which scales as p 4 q 4 /N 2 , is actually included in the Q-Hermite result which helps explain why this approximation, with a difference from the exact result that is also subleading in p, goes beyond its natural limit of applicability. However, this does not help explain why the full exact 1/N 2 correction to the density gives worse results than the Q-Hermite approach such as spurious 1/N 2 corrections to the density close to the ground state. In future work we plan to address some of these problems.

A Calculation of the 1/N 2 corrections
The starting point is the general formula (4.8) for η G , which we replicate here for easier access of reading: We already argued in eq. (4.11) that the orders in 1/N is controlled by the term Note that because of cancellations, the sum (A.1) may be of higher order in 1/N , for example if a diagram contains nested contractions. Hence the following four cases contribute to the order of 1/N 2 : • d 2 = 0, d k≥3 = 0; • d 2 = 1, d k≥3 = 0; • d 2 = 2, d k≥3 = 0; We will compute each of the cases and sum them up. Note the summation indices in (4.8) are d α i 1 ···α i k , which involves all the k-vertex structures in an intersection graph G. For example, the sum over d α i 1 α i 2 involves summing over all edges that connect arbitrary two vertices in G, and there are V (V − 1)/2 such edges. Hence, the edges we need to sum over are more than just the edges of G itself, and to aid the forthcoming computation, we complete the graph G by adding dashed lines between all vertices that are not connected by a solid line. In graph theory this is known as the completion of a graph. In figure 7 we illustrate two examples of such graph completion. For a completed graph, we denote by w 0 , w 1 , w 2 the numbers of wedges with 0, 1 and 2 solid lines, and by n 0 , n 1 , n 2 , n 3 the numbers of triangles with 0, 1, 2 and 3 solid lines (see table 4), which will be useful for the computations to come.
−→ −→ Figure 7. Graph completion for a 3-vertex graph and 4-vertex graph. The resulting graphs have edges between each pair of vertices. This is equivalent to the edge-2 coloring (i.e. that we use two different colors to color a graph) of complete graphs (i.e. graphs where each pair of vertices is connected by an edge).
A.1 Cases with d k≥3 = 0 In this case the general expression for contractions reduces to Note the second equality of (A.3) holds because we are working with the particular case d k≥3 = 0. To order 1/N 2 we have three cases, d 2 = 0, d 2 = 1 and d 2 = 2 which we will analyze next.
Case d 2 = 0. The result (A.2) for a contraction diagram simplifies to Case d 2 = 1. In this case, among the V 2 possible d α k α l we sum over, E of them occur in the set of edges of G (solid lines in a completed graph) and they give c(G) = 1, while V 2 − E of them are not in this set of edges (dashed lines in the completed graph) and they give c(G) = 0. So this case contributes to η G by the following expression: Case d 2 = 2. When one of the d α k α l = 2 then either c(G) = 2 or c(G) = 0 and in both cases the phase factor is equal to 1. This results in the contribution (A.6) The case with two of the d α k α l equal to 1 is more complicated because we have to distinguish the case where they share a common index and the case where they do not. If they do not share a common index the combinatorial factor apart from the phase factor is while when they share a common index we obtain The coefficient of the q 4 term is the same in both cases which simplifies the counting. We have pairs with c(G) = 0. Summing over the {d α k α l } also including pairs that share a common index, adding the −q 3 /N 2 contribution from (A.8) which was not accounted for in the previous counting, we obtain to order 1/N 2 where we have separated the contribution of pairs that share a common index, because it combines naturally with the contributions that will be discussed in the next subsection. The w 0 , w 1 and w 2 are defined in table 4, they will give c(G) = 0, 1, 2, respectively, and hence the signs in (A.10). Summing all the terms with d k≥3 = 0, i.e. eqs. (A.4), (A.5) and (A.9), the result simplifies to A.2 Cases with d 3 = 1 When d 3 = 1 to order 1/N 2 only the d 2 = 0 terms contribute. The result for a contraction diagram is then given by The second equality (A.13) is true because we are working in the particular case where d k =3 = 0, and we remind the reader that there is another summation implied by the " * " in the subscript. Either 0, 1, 2 or 3 of the edges of d α k α l αm can be part of the edges that occur in c(G), and their total number are n 0 , n 1 , n 2 or n 3 respectively, as defined in table 4. We thus obtain the contribution Including the contribution of the wedges in eq. (A.10), the total q 3 /N 2 contribution is given by From the graphical interpretation of these quantities (see Table 4), it is clear that The reason is that each wedge is contained in one and only one triangle, and the identities follow by counting in each type of triangle how many wedges of different types occur. This results in the simplification We have the obvious identities [65] n 1 + 2n 2 + 3n 3 = (V − 2)E, (A.18) The first identity can be seen as follows. If we have an edge, it can be combined with 0, 1, or 2 other edges to from a triangle with one solid edge, two solid edges, or three solid edges, respectively. In total there are E(V − 2) possibilities to combine any edge with the remaining vertices, which gives the right hand side of this identity. The left-hand side counts the same thing triangle by triangle. That is, a triangle with one solid edge occurs once, a triangle with two solid edges is counted twice because each of the solid edges can be taken as the starting point in E. For the same reason a solid triangle is counted three times in (V − 2)E, and hence the identity. By subtracting the two identities we find The total contribution of the wedge and triangles can thus be written as Adding the thus contribution to the d k≤2 contributions (see eq. (A.11)) we obtain Comparing this to the Q-Hermite result, the 1/N 2 correction simplifies to This proves (3.7).

B Scaling limit of η
To obtain the large N double scaling limit of η at fixed q 2 /N we express η in terms of the hypergeometric function as The double scaling limit of the binomial factor is given by The hypergeometric function u := 2 F 1 (−q, −q, N + 1 − 2q, z) satisfies the differential equation In the double scaling limit this simplifies to which is solved by The constant is fixed to c = 1 by the requirement that η = 1 for z = 1. For z = −1 we thus obtain the asymptotic double scaling limit

C Edge counting from the Riordan-Touchard formula
To evaluate the 1/N and 1/N 2 contributions to the Q-Hermite moments we need the sums They follow from the first and second derivatives of the Riordan-Touchard formula at η = 1 where the sum is over all contractions for the 2p-th moment, and E i is the number of crossing for the i-th chord diagram. Below we will show that The first two equalities were already shown in [57]. Here, we give the key ingredients of the proof. We start from the following integral representation of φ p , where This can be used to show that φ p (1) = (2p − 1)!!. It also follows To obtain the derivatives of φ p (η) we expand H(x, t) p about t = 0, This proves the first equality (C.3). The second equality (C.4) can be similarly proved.
The proof of the last two equalities (C.5) and (C.6) requires some more work. Using the integral representation (C.7) of φ p (η) we obtain where we have substituted x = √ 2y for the last equality. It is straightforward to show Hence, which proves (C.5).
To prove the last equality, we start with The second derivative can be expressed in terms of the integral representation (C.7) of φ p (η) as (C.18) After some tedious but straightforward manipulations, we get where we have used two additional integrals to obtain the last line.

D Cut-vertices and factorization
Since the subscript space of N q elements is isotropic we can always fix one index and the result of a diagram does not depend on this index. So summing over this index gives a factor N q .
We define a cut-vertex as a vertex that when it is cut, the graph becomes disconnected. A graph without any cut-vertex is called a non-separable or a two-connected graph. If we apply the reasoning of the first paragraph to a cut vertex, we immediately arrive at the theorem Theorem 1. If a graph G contains a cut-vertex, which separates G into subgraphs G 1 and G 2 , then η G = η G 1 η G 2 .
As an example, the following graphs all contain one cut-vertex (drawn in red): Theorem 1, for example implies, In this appendix we calculate the moments M 2p /M p 2 to order 1/N 2 for the q = 2 SYK model. Since we are interested in large N asymptotics, it is immaterial whether N is even or odd, and for technical simplicity we choose N to be even. The Hamiltonian for q = 2 model is given by This can be rewritten as [27] where x k are the positive eigenvalues of the antisymmetric matrices J ij , and c k , c † k are the annihilation and creation operators for Dirac fermions. We have Thus we can compute the ensemble average by averaging over the joint probability distribution of anti-symmetric Hermitian random matrices [66,67], Here c is a normalization constant, which we have chosen such that P (x 1 , . . . , x N/2 ) is normalized to unity. Note that the joint probability distribution has the parity symmetry x k → −x k for each individual variable, so we can we can restrict ourselves to the configuration with only positive s k by compensating with an overall factor 2 N/2 . This results in In particular where we have used x k = 0 and x 2 1 = x 2 2 = . . . = x 2 N/2 because of the parity symmetry and permutation symmetry of P (x 1 , . . . , x N/2 ). Hence, To isolate the leading orders in 1/N we now need to analyze the terms in that are leading orders in N . Again, because of the permutation symmetry of P (x 1 , . . . , x N/2 ), the expectation values on the right-hand side of equation (E.9) only depend on the partition of p into {m 1 , m 2 , · · · , m N/2 }. Therefore, given a partition with k nonzero m i 's, {m i 1 , m i 2 , · · · , m i k } (i 1 < i 2 < · · · < i k by convention), all the expectation values of the form have the same value and this results in a multiplicity factor N/2 k from choosing k-element subsets of {x 1 , . . . , x N/2 }. It is not too hard to convince oneself that so we can isolate the leading terms in 1/N from the binomial factors N/2 k , and there can be more multiplicity factors from permuting the x j 's in the same partition, but that does not bring any factors of N . Now it becomes obvious that the leading terms which are relevant recursive relations developed in [68]. After some algebra we obtain (E. 16) Hence, to the relevant orders, we have Note that the each leading term is just the number of nested contractions when re-expressing the W k in products of TrJ 2k to leading order in 1/N . Substituting the above results into (E.15), we finally arrive at

F Calculation of the eighth moment
The sixth moment was already discussed in [42] and in this appendix we only quote the final result. For the eighth moment we work out all contributions explicitly. Although originally the combinatorics for the eighth moment were obtained by inspection, in several cases we also show how they can be obtained from the general formula (4.8).
According to the Riordan-Touchard formula we have For p = 3 the right hand side is given by Only the η 3 term deviates from the exact result which is given by For p = 4 the Riordan-Touchard formula yields 14 + 28η + 28η 2 + 20η 3 + 10η 4 + 4η 5 + η 6 , (F. 6) where the coefficient of η E gives the number of diagrams with E crossings. For one and two For E = 3 we have two different contributions. The first class of twelve diagrams is shown in figure 11, where we can move the Γ α to three consecutive pairs with equal indices by three independent pair exchanges. The results of each of these diagrams is given by η 3 . The second class of eight diagrams with three intersections contains a structure we first saw in the calculation of the sixth moment (see figure 12). It is given by T 6 [42], see eq. (F.4) and we thus obtain For E = 4 also two different classes of diagrams contribution to the eighth moment. The two diagrams of the first class are shown in figure 13. The result depends on how many       gamma matrices the second and third contraction have in common. It is given by The second class of eight diagrams has one intersecting contraction which can be removed by a pair exchange, and three other contractions which have three intersections and contribute as T 6 (see figure 14). The total contribution of diagrams with four intersections is thus given There are four diagrams with 5 intersections which are shown in figure 15. The contribution of these diagrams cannot be decomposed into contributions we have seen before. The result of each diagram is given by The contribution of diagrams with five intersections to the eighth moment is thus given by (F.12) The most complicated diagram has six intersections, see figure 16. The result is given by Using the convolution property of binomial factors, this can be simplified to (F.14) The contribution to the fourth moment with six crossings is thus given by We have checked the above results in several ways. First, when the phase factor is eliminated the contribution of each diagram evaluates to one. Second, the moments agree with numerical result of the moments obtained from the eigenvalues of the SYK Hamiltonian at finite N . Third, for q = 1 the results agree with the exact analytical result in eq. (5.9).
In the next subsection we derive the result for T 6 and T 8 from the general formula (4.8). The result of the diagram (c) of figure 1 using the general formula (4.8) is given by which is the result obtained in section 7.