Large $N$ expansion of the moments and free energy of Sachdev-Ye-Kitaev model, and the enumeration of intersection graphs

In this paper we explain the relation between the free energy of the SYK model for $N$ Majorana fermions with a random $q$-body interaction and the moments of its spectral density. The high temperature expansion of the free energy gives the cumulants of the spectral density. Using that the cumulants are extensive we find the $p$ dependence of the $1/N^2$ correction of the $2p$-th moments obtained in 1801.02696. Conversely, the $1/N^2$ corrections to the moments give the correction (even $q$) to the $\beta^6$ coefficient of the high temperature expansion of the free energy for arbitrary $q$. Our result agrees with the $1/q^3$ correction obtained by Tarnopolsky using a mean field expansion. These considerations also lead to a more powerful method for solving the moment problem and intersection-graph enumeration problems. We take advantage of this and push the moment calculation to $1/N^3$ order and find surprisingly simple enumeration identities for intersection graphs. The $1/N^3$ corrections to the moments, give corrections to the $\beta^8$ coefficient (for even $q$) of the high temperature expansion of the free energy which have not been calculated before. Results for odd $q$, where the SYK `Hamiltonian' is the supercharge of a supersymmetric theory are discussed as well.


Introduction
The statistical fluctuations of nuclear levels have been successfully described by the Gaussian Orthogonal Ensemble (GOE). However, the average spectral density of the GOE is a semicircle which is very different from the Bethe formula [2]. In agreement with experimental observations [3], this formula predicts an exp(c √ E − E 0 ) dependence (with E 0 the ground state energy) on the excitation energy E. Moreover, the nuclear interaction is mostly a two-body interaction while for the GOE all the many-body states interact. To address these shortcomings, French and co-workers [4][5][6][7] introduced the two-body random ensemble which is now known as the four-body complex SYK model. One of the first results for this ensemble was that the level density is a Gaussian [6] which is closer to the expectation of realistic many-body systems than the semicircular behavior. However, in the nuclear physics community it was not realized that this model actually reproduces the Bethe formula [8]. One of the reasons for missing this opportunity was the custom [9][10][11][12][13] to study this model as the sum of a two-body and a four-body interaction (in the nuclear physics convention, a one-body and a two-body interaction). The reason is that the nuclear interaction was seen as a residual random four-body interaction on top of a mean field which can be represented as a two-body interaction. Note that the four-body interaction is a irrelevant term with regard to the two-body interaction [14][15][16]. It was also understood early on that the two-body random ensemble still has the level correlations of the Gaussian Orthogonal Ensemble [7]. However, it was also realized that the model was not ergodic [16][17][18][19] in the sense that ensemble average of a spectral correlator is not equal to the spectral average of a spectral correlator, the latter given by the result of the Wigner-Dyson ensembles up to much larger distances. For more discussions of the two-body random ensemble in nuclear and many-body physics we refer to [9,18,20,21].
In condensed matter, the model was introduced independently as a random quantum spin model [22] where a mean field is not natural. In this context, Sachdev and co-workers discovered [23] a remarkable property of this model, namely that its zero temperature entropy is extensive which then was identified as the black hole entropy [24]. This property is directly related to its exponentially large level density starting with the ground state region, which is a non-Fermi liquid. The states are characterized by highly entangled states [25] which are very different from the particle-hole excitations of a Fermi liquid.
In the past two years the interest in this model was rekindled because it is possibly dual to 1+1 dimensional gravity [8,. The properties that made this model attractive as a model of a compound nucleus are exactly those which are required for the existence of a black-hole dual: it is maximally chaotic [47] with spectral correlations given by random matrix theory [16,33,48], it has level density given by the Bethe formula which also implies that the specific heat is linear in temperature for low temperature, the zero temperature entropy is extensive showing that the low-lying states of the model are strongly entangled. What is particularly important in this context is the existence of a conformal limit [26,27] where the action of the SYK model reduces to the Schwarzian action [27,42,49,50].
Much of the recent progress on the SYK model was made possible because of the path integral formulation [22,26,27], from which one can derive an exact result for large N limit after averaging over the random interaction. This was not possible for a formulation that started from the generating function for the resolvent [51]. Although it was straightforward to average over the randomness, it resulted in a complicated theory that was not amenable to taking a large N limit. The disadvantage of the path integral which at the same time is its strength is that it provides access to the Green's function rather than the level density. Since our main interest has been in the level density and the level correlations [1,8,14,16,33] of the SYK model, we have used the moment method [52] which also proved effective in the early applications to nuclear physics. Two limiting cases for N fermions fermions with a random q body interaction were easily recognized: q ∼ N and q N . In the first case [18,53] the SYK model is in the universality class of the Wigner-Dyson ensembles with a semicircular spectral density. In the second case, the spectral density is a Gaussian [6]. This suggests the existence of a double scaling limit which converges to a spectral density in between a semicircle and a Gaussian. Indeed this happens when q ∼ √ N for N → ∞ [54]. This scaling limit which reveals itself in the path integral formulation of the SYK model was not noted before in the nuclear physics literature either. In the moment method, this limit arises naturally in Wick contractions when treating all intersections as independent which gives moments that are exact to order 1/N [8,33,38,54]. Remarkably, these moments turn out to give the spectral density of the weight function of the Q-Hermite polynomials with a nontrivial double scaling limit [8,38,54,55].
The 1/N 2 corrections to all moments can be calculated analytically as well [1], with a result that has as a simple geometric interpretation. The p-dependence of the 2p-th moment also turns out to be relatively simple. One of the main goals of this paper is to explain this p-dependence of the moments. Since the high temperature expansion of the free energy is the cumulant expansion of the spectral density, the extensivity of the free energy puts strong constraints on the moments. In fact, we will show that the p-dependence of the moments follows almost entirely from this condition, and that it is determined by a few low-order moments only. Secondly, we study the way large N corrections to the moments contribute to the free energy. We already know that the Q-Hermite moments give the free energy for all temperatures to order 1/q 2 [8,56]. In this paper we will show that 1/q 3 corrections to the free energy follow from the 1/N 2 corrections to moments. Thirdly, using the inter-relationship between the free energy and the moments, we obtain the 1/N 3 corrections to the moments which are responsible for the 1/q 4 correction to the high temperature expansion of the free energy. Using these results we find new enumerative identities for intersection graphs. Results for the supersymmetric SYK model [57][58][59] which can be derived in a similar way are also given in this paper. This paper is organized as follows. In section 2, we give a brief review of the SYK model including the moment method. The relation between moments and the free energy is discussed in section 3. In this section we show that at a given order in 1/N the p-dependence of the moments follows from a few low-order moments . We also show that large N corrections to the moments give large q corrections to the free energy. Results are obtained for both even q and odd q. As a new result we obtain the 1/N 3 corrections to the moments and 1/q 4 corrections to the free energy. In section 4, we use the results for the moments to derive new enumerative identities for intersection graphs. Some technical results are deferred to two appendices. In appendix A, we evaluate the high temperature expansion to the free energy from the results of Tarnopolsky [56]. We obtain the high temperature expansion of Tarnopolsky's result to all orders and show that it is a convergent series with no singularities on the positive real axis. Finally, we derive in appendix B the 1/N 3 corrections for q = 1 and q = 2 with an independent method.

Review of SYK model and moment method
In this section we introduce the Sachdev-Ye-Kitaev (SYK) model and the moment method.

The SYK Hamiltonian
The SYK model is a system of N Majorana particles with the q-body interaction represented by and γ α are the Euclidean gamma matrices with commutation relations: The factors of i in the definition of Γ α have been included so that H is Hermitian. The sum is over all N q q-particles states denoted by the collective index α = {i 1 , i 2 , . . . , i q }, with 1 ≤ i 1 < i 2 < · · · < i q ≤ N . The couplings J α are Gaussian distributed: where the parameter J sets a physical scale. 1 For even q the Hamiltonian H of the SYK model is simply given by With the large N scaling of the variance in (2.4) this Hamiltonian has a negative-energy ground state energy that is proportional to N for large N .
For odd q the operator H is still a well-defined Hermitian operator, but because it has a fermionic grading, it is not a Hamiltonian but rather the supercharge of a supersymmetric theory with the Hamiltonian [57]: This Hamiltonian is positive definite with a ground state energy approaching zero in the thermodynamic limit with the scaling of the variance as in (2.4).

Moments and spectral density
One of the reasons for which we are interested in moments is to study the spectral density ρ(E) of the SYK model: 2 where · · · denotes the Gaussian average using the distribution (2.4). By a Fourier transform, we can express ρ(λ) as where we have defined the k-th moment M k to be The third equality of eq. (2.8) used the fact that M k = 0 when k is odd, due to the J α → −J α symmetry of the distribution (2.4). Therefore, we may focus on the calculation of M 2p . Notice that since · · · is an Gaussian integration over J α , M 2p is given by a sum over all possible Wick contractions among 2p Γ's. By convention, we have factored out the Hilbert space dimensionality 2 N/2 , but M 2p still grows like N p as N → ∞, as can be seen from the N dependence of M 2 : Hence, to formulate a useful large N expansion for moments, we consider instead the scaled moments:M We distinguish the moments of the SYK operator H, which will be denoted by M 2p and the moments of the Hamiltonian H which will be denoted by µ p . So we have that for even q, (2.12) 3 Free energy, cumulants and high temperature expansions

Moments and cumulants
The partition function is given by with the free energy denoted by F . In terms of high temperature expansion of the free energy, we have The quantity κ n is called the n-th cumulant. Note that the summation starts from n = 1, this is because the energy is finite, and at infinite temperature we expect only the entropy to contribute to βF . Since free energy is extensive, we obtain for all n. We will see later that this has important consequences for the N -dependence of the moments.
Alternatively, the partition function can also be expressed in the moments µ n of H, One may ask why we do not consider the "free energy" e −βF := Tr(e −βH ) also for odd q so that we can have a uniform treatment of the moments for all values of q, instead of two separate cases eqs. (2.12) and (2.13). The problem with this "free energy" is that it is not extensive. Extensivity will turn out to be vital for the application of our method.
The relation between µ n and κ n is well studied [60]. Consider the following partition of an integer n: where by convention we demand the ordering k 1 ≥ k 2 ≥ · · · ≥ k l , 3 then we have where Pn denotes sum over all partitions of n. As examples, we list some low-order relations: (3.8) We remark that since moments are computed by contracting Γ matrices, all the N -dependence comes from counting the subscripts of those Γ matrices. This implies the µ n must be rational functions of N , and then eq. (3.7) tells us the cumulants κ n must also be rational functions of N . Hence there can be no factors such as log N or √ N in the large N expansion of moments or cumulants. For even q all odd cumulants vanish, but they enter in the calculations for odd q.

N dependence of cumulants and moments
As discussed, for even q we have M k = µ k and hence µ k = 0 when k is an odd number. It follows from eq. (3.7) that κ k = 0 when k is odd. We are interested in the scaled moments M 2p , so let us also define the scaled cumulants to bẽ According to eq. (2.10), (3.10) We also have the trivial identityκ Since κ 2 = µ 2 it is clear that eq. (3.6) is also valid for the rescaled moments and cumulants. Because of the N dependence eq. (3.10), the corrections of order 1/N k−1 to all scaled moments only receive contributions from cumulants up toκ 2k , which by themselves are completely determined by the moments up to order 2k due to eq. (3.7). Hence we conclude: To order N −k+1 , allM 2p are determined by a finite number of moments up toM 2l (l ≤ k), expanded to N −k+1 .
For example, the 1/N 2 expansion of all moments is completely determined byκ 4 and κ 6 which in turn are determined byM 4 andM 6 , while at order 1/N 3 the moments receive only contributions fromκ 4 ,κ 6 andκ 8 and are thus determined byM 4 ,M 6 andM 8 . This is very surprising, because for large p, the Wick contractions contributing to the 2p-th moment become rather complicated, whereas to calculate up toM 8 we only need to consider a small set of contraction diagrams. This will have important implications when it comes to the enumeration of intersection graphs, which is the subject of section 4.
The scaled cumulantκ 2k is O(1/N k−1 ), but after being rescaled back to κ 2k its leading term contributes to the thermodynamic limit of the free energy. The leading term ofκ 2k is determined by the 1/N k−1 corrections of the moments up to order 2k. We thus emphasize: The leading term of the 2k-th scaled cumulantκ 2k is determined by the momentsM 2l (l ≤ k), expanded to order N −k+1 .
This in particular implies that even if we want the complete information of only the thermodynamic limit (leading in 1/N ) of the free energy, we would still need all-order knowledge of scaled cumulants and hence of the scaled moments.
We will see below that the q-dependence of the leading term of the sixth cumulant follows from the 1/N 2 corrections to the fourth and sixth moment. The correction factor is simply given by 1 − 1/3q with no other corrections from higher moments. So the large q expansion of this cumulant terminates at this order.
The discussion in this section is general and also applies to the odd q case, where only some minor changes of notations are needed.

Explicit results to 1/N 3
We will derive the expansion ofM 2p to 1/N 3 in this paper. To the relevant order, eq. (3.6) can be explicitly written as (3.12) We emphasize again the fact that only a finite number of terms need to be considered on the right-hand side of the above equation is due to the extensivity of the free energy, see eqs.
(3.3) and (3.10). The first eight moments were calculated exactly in [1], and their expansion up to order 1/N 3 thus gives the expansion of all moments to this order.
Using the results of [1] we obtaiñ This results in the cumulants, which have exactly the leading order N dependence of eq. (3.10) required for an extensive free energy.
Substituting these results for cumulants back to eq. (3.12), we obtain the following large N expansion for the scaled moments to 1/N 3 : Trivially known from σ 2 All zero due to extensivity ? Figure 1. Relation between two expansions (for even q). The horizontal axes denote orders in 1/N , while the vertical axes denote orders in (scaled) moments and cumulants. The red dots are where the corresponding 1/N coefficients are known, and the gray region is unknown. It is interesting to see that on the moment side infinitely many coefficients are known, but they are all determined by the finite number of dots in the dashed triangle. This is ultimately because on cumulant side we easily know an infinite number of coefficients due to extensivity (the cyan region).
it agrees with (3.15), which we just obtained by a completely independent method. The discussion above suggests an interesting "map of knowledge" between the moment expansion and cumulant expansion, as shown in figure 1.

Free energy
We can substitute eq. (3.14) into eq. (3.2) to obtain the high temperature expansion of the free energy: Since σ 2 ∼ N we obtain in the thermodynamic limit where we have set J 2 = 2 q−1 /q. Using the 1/N corrections to the cumulants and the variance, it is straightforward to calculate 1/N corrections to the free energy. We have sufficient data to obtain terms up to 1/N 3 , each expanded to β 8 , but since they are of limited physical relevance, we have not written them down. To summarize, we have obtained the high temperature expansion of the free energy to order β 8 .
Using a completely different method, a recent publication [56] computes the large q ex-pansion to order 1/q 3 for the free energy at leading order of 1/N . 4 In addition to reproducing the large q expansion calculated by [56] to order 1/q 3 (see appendix A), we also obtain the 1/q 4 correction to the free energy at order β 8 which is given by Our results also show that there are no further large q corrections for terms up to order β 8 .
We also remark that, if all one wants is the high temperature expansion of the free energy to a certain order in β and 1/N , a general expression for the large N expansion ofM 2p such as eq. (3.15) is not necessary, as only equations such as (3.14) are used to calculate the high temperature expansion, where only a finite number of moments are needed. However, the scope of this paper is wider than just computing the high temperature expansion, and a general expression for allM 2p will be needed to solve the enumeration problem for intersection graphs, which will be discussed in section 4.

Odd q case
We can repeat the same calculation for supersymmetric SYK models. However in this case we need to define the scaled cumulants as because M 2 = µ 1 = κ 1 for the odd q case. This means Repeating the same calculation that led to eq. (3.12) we get There is no (2p − 1)!! factor this time. To calculateκ 2 ,κ 3 andκ 4 we will needM 4 ,M 6 and M 8 as we did for the even q case. They can be again calculated from the analytical results in [1] resulting iñ Using the relation between cumulants and moments, eq. (3.7), we obtaiñ (3.23) Together with (3.21), we conclude that the p-dependence of the moments is given bỹ The high temperature expansion of free energy has the form: More explicitly, we have in the thermodynamic limit where again we have set J 2 = 2 q−1 /q. We only displayed the coefficient of the leading term in 1/N , but the 1/N corrections can be calculated in a straightforward way as well.

Enumerative identities of intersection graphs
In this section we derive enumerative identities for intersection graph, some of which have not appeared in the literature. We start with a short review of the graphical calculation of the moments, and the details can be found in [1].

Graphical calculation of moments
The moments of the SYK model are given by the expectation value Because the probability is Gaussian, the average is given by the sum over all possible Wick contractions, which can be represented by rooted chord diagrams, and some examples of such chord diagrams are given in the first row of figure 2. For large N and finite p, the indices of the Γ α are mostly different, so that they can be commuted to pairs of Γ α Γ α = 1 (no implicit summation over α). Since in this limit all contractions contribute equally for even q, this results in a Gaussian spectral density. For odd q, when two Γ α and Γ β with no common indices anti-commute, the contractions are alternately positive and negative, leaving only one net contraction for all moments which results in the moments of two delta functions in this limit. We always consider the scaled momentsM 2p := M 2p /M p 2 , so that the variance cancels in the ratio, and the values of chord diagrams always refer to the values of Wick contractions normalized by M p 2 . To calculate 1/N corrections, we have to take into account that Γ α and Γ β commute or anti-commute depending on how many indices they have in common, where k is the number of common indices in α and β. Taking this into account, we obtain the value of two intersecting contraction lines: We can further translate chord diagrams into intersection graphs, which are obtained by representing each chord by a vertex, and connecting two vertices if and only if the two chords they represent intersect each other in the chord diagram. We give some examples of intersection graphs in the second row of figure 2. We will denote a generic intersection graph by G. An important approximation to the moments is when all intersecting contraction lines are treated as being independent. In the language of intersection graphs, if an intersection graph G has E edges, its value η G is approximated by This approximation gives moments that are correct up to order 1/N . The corresponding moments are the moments of the weight function of the Q-Hermite polynomials. That is why this approximation is known as the Q-Hermite approximation. The Q-Hermite approximation to moments is thus the sum of η E over all the (2p − 1)!! intersection graphs.
For more details of the graphical calculations including the graph-theoretic identities needed to sum all graphs we refer to [1]. To conclude this review, we remark that the method of [1] relies on a back-and-forth interplay between the moment expansion and intersection graphs: intersection graphs inform us on how to calculate moments, and moment calculations for the exactly solvable q = 1 and q = 2 SYK model prove enumerative identities about intersection graphs, which in turn feed back to the moment calculation for general q. In the current paper the relation is more one-way: we have obtained results for general moments in previous sections without relying on the enumeration of intersection graphs, and the graphs at best could play a minor role as a book keeping device. However, now we can use the results for the moment expansion to prove identities for intersection graphs, which will be discussed in this section.

Structure of contributions
In section 3, we have calculated the large N expansions ofM 2p to order 1/N 3 . We did not use any enumerative identities like the ones in [1]. Nevertheless, we see that various binomial factors arise in eqs. (3.12) and (3.21) simply from the relations between the moments and the cumulants. This suggests we can turn around and use the calculations we just presented to generate enumerative identities for intersection graphs. To see what type of graph-theoretic objects are to be enumerated, we first state the main result of the next two sections: where E is the numbers of edges, T is the number of triangles, and f 6 , f 5 , f 4 are the number of the four-point structures depicted in figure 3, in an intersection graph G. We can check this formula for a few nontrivial graphs that contribute to moments up toM 8 denoted by η, T 6 , T 44 , T 66 and T 8 in tables 1 and 2. Expanding the results for these graphs obtained in   [1] to order 1/N 3 we obtain, (4.6) Using the graphs in tables 1 and 2 one can easily verify that the above results satisfy eq. (4.5). Note that triangles (whose value is T 6 ) made their first appearance as a complete intersection graph for the sixth moment (table 1), and in the eighth moment, they become substructures of various graphs. The same is true for the four-point structures whose values are T 44 , T 66 and T 8 : they first appear as complete graphs for the eighth moment (table 2), and will become substructures for higher moments. Some examples are given in figure 4.
How do we proceed to prove eq. (4.5) in its full generality? A rigorous proof of eq. (4.5) to order 1/N 2 was given in [1], which is essentially a very tedious brute-force calculation. The  1/N 3 term can be obtained rigorously by the same method, but the calculation turns out to be extremely tedious and uninstructive. So instead, in the next two sections, we will present a less rigorous but more instructive method to obtain the result (4.5). As a warm-up we will start with the 1/N 2 term in the next section, and then continue with the 1/N 3 corrections.

Graphical calculation of the 1/N 2 term
We first take note of a rather trivial fact that if an intersection graph G is disconnected with two components G 1 and G 2 , then η G = η G 1 η G 2 . Moreover, this factorization property also holds in a less trivial situation, where an intersection graph can be made disconnected by cutting a vertex, which was proved in [1]. 5 To be concrete, let us first look at the 1/N 2 order corrections. We have argued that the 1/N 2 coefficient forM 2p is completely determined bỹ M 4 andM 6 , which only depend on the values of edges and triangles (η and T 6 ) according to table 1 and factorization. However, this fact appears rather mysterious if one thinks about the moments in terms of the Wick contractions where G i denotes the i-th intersection graph and all the intersection graphs have p vertices. An intersection graph can get quite complicated when p is large, for example imagine a complete graph with p vertices, which has p 2 edges. On top of that one then needs to sum over all these complicated graphs. How can one tame such a complex beast by only using M 4 andM 6 ? The only plausible way to reconcile these two seemingly conflicting view points is that the 1/N 2 coefficient of each η G i ought to be a function of only the number of edges and triangles in that intersection graph, not something more global or any other property of the graph. So the fact that we are only looking at a fixed order in 1/N allows a huge simplification to happen. We may summarize this plausible result as where A(E, T ) is some function of E and T . 6 Note that since η E is a function of E, the plausibility argument loses no generality by considering η G − η E instead of η G on the lefthand side of eq. (4.8), and this is motivated by the Q-Hermite approximation eq. (4.4). Let us pause here and summarize the reasons for the necessity for such simplicity: • The cut-vertex factorization property of η G , which together with table 1 impliesM 4 and M 6 depend only on the values of edges and triangles (η and T 6 ).
• The extensivity of the free energy/cumulants, which forces the 1/N 2 coefficient ofM 2p to depend on onlyM 4 andM 6 , for any value of p.
Now the task is to fix the explicit form of A(E, T ). It is sufficient to consider those special graphs G with only disconnected irreducible structures, so that there is a complete factorization. Since the relevant irreducible structures can only be edges or triangles, let us take G to be an intersection graph with E edges and T triangles where all the triangles are disconnected, and the edges other than the ones that make up triangles are also disconnected, which means there are E − 3T of them. This implies where T 6 is the value of the Wick contraction represented by a triangle. We now have We define By the expansion of T 6 in eq. (4.6), we have δT Using that η = (−1) q + O(1/N ), (4.13) we finally obtain Hence by considering special graphs we have fixed the form of A(E, T ) to be A(E, T ) = (−1) Eq (−8q 3 )T , which must also be true for a general graph. This proves eq. (4.5) to order 1/N 2 .

Graphical calculation of the 1/N 3 term
We may continue with this strategy to order 1/N 3 . SinceM 2p to this order is completely determined byM 4 ,M 6 andM 8 , by the same argument as in previous subsection, we conclude for any intersection graph G we have because E, T, f 6 , f 5 and f 4 count all the irreducible structures that appear for moments up toM 8 . Again consider an intersection graph G with E edges, T triangles, and the four-point structures counted by f 6 , f 5 and f 4 (see figure 3), where all the irreducible structures are disconnected, so we have T − 4f 6 − 2f 5 disconnected triangles and E − 3(T − 4f 6 − 2f 5 ) − 6f 6 − 5f 5 − 4f 4 = E − 3T + 6f 6 + f 5 − 4f 4 disconnected edges. Thus where T 8 , T 66 and T 44 are the values of the intersection graphs counted by f 6 , f 5 and f 4 , see table 2. We have that The above equations should be understood as definitions of the symbols δη, δT 6 , δT 8 δT 66 and δT 44 from the left-hand sides. From eq. (4.6) we see δT 8 , δT 66 and δT 44 are of order 1/N 3 and δT 6 contains both 1/N 2 and 1/N 3 terms. To order 1/N 3 the only mixed contribution is of the form δηδT 6 , while the other corrections only contribute by their leading orders. We thus obtain If we use the explicit expressions for the irreducible structures given in eq. (4.6), we finally find This fixes the form of B(E, T, f 6 , f 5 , f 4 ) and proves eq. (4.5).
To calculate the correction to the Q-Hermite moments we sum over all intersection graphs: It is now clear that the new objects to be enumerated at order 1/N 3 are (−1) Eq ET and (−1) Eq (5f 6 + f 5 − f 4 ). We will discuss their enumerations for the even q and the odd q cases in the following sections.

Even q case
In the same way thatM 2p to order 1/N 3 is determined byM 4 ,M 6 andM 8 , the Q-Hermite momentM QH 2p to order 1/N 3 is determined byM QH 4 ,M QH 6 andM QH 8 . These Q-Hermite moments can be calculated most easily calculated from the definitioñ 22) the large N expansion of η in eq. (4.6), and tables 1 and 2. Alternatively, the sum in (4.22) may be evaluated by the Riordan-Touchard formula [1]. We find the moments, (4.23) In the same manner as in eq. (3.14), we get the following Q-Hermite cumulants: (4.24) The p-dependence of the moments is then obtained from eq. (3.12): In fact, using the definition of the Q-Hermite moments (4.22) and we can match the coefficients of powers of q and 1/N in eqs. (4.25) and (4.26), and this already gives several enumerative identities for E:

Odd q case
We repeat the same calculation for odd q case. The Q-Hermite moments up to eighth order areM (4.34) Together with eq. (3.23) this gives the Q-Hermite cumulants Substituting this into eq. (3.21) and we find Again, givenM QH 2p = i η E i and the form of η in eq. (4.26), we have the following enumerative identities by matching the coefficients of the expansion in 1/N and q : These edge identities could also have been obtained by the analytical-combinatorial approach used in [61]. Now we calculate the difference betweenM 2p andM QH 2p : By matching with eq. (4.5), we obtain the following graded identities: These identities could not have been obtained by using the method in [61]. Finally we discuss in what sense the method presented in this paper is more powerful than the method developed in [1] to compute moments of the SYK model. The idea there is to solve forM 2p −M QH 2p for q = 1 and q = 2 models, where moments can be evaluated exactly, and then do a matching with (η G i − η E i ) expressed in terms of graph-theoretic objects. One can get a flavor of this method from appendix B. The old method works very well to order 1/N 2 , however it becomes problematic starting from order 1/N 3 . If we look at eq. (4.5), we see that at order 1/N 3 , on top of triangles, there are two new types of structures that need to be enumerated: (−1) Eq ET for the q 5 term and (−1) Eq (5f 6 + f 5 − f 4 ) for the q 4 term. If we compute the moments for q = 1 and q = 2 models, at best we can obtain the enumeration of a linear combination of the two above-mentioned new structures, that is, (−1) E (ET − 5f 6 − f 5 + f 4 ) for q = 1 and 2ET − 5f 6 − f 5 + f 4 for q = 2. However to recover the full q dependence, we need separate enumerations of (−1) Eq ET and (−1) Eq (5f 6 +f 5 −f 4 ). The method presented in this paper faces no such difficulty since the full q dependence is retained at every stage of the calculation, and thus is capable of enumerating the two new structures separately. Nevertheless, the old method provides an independent consistency check of the results obtained in the present paper, which will be demonstrated in appendix B.

Conclusions and outlook
We have established the relation between the 1/N expansion of SYK moments and the 1/N expansion of the high temperature expansion coefficients of the free energy. It turns out that Table 3. Some of the enumerative identities for intersection graphs generated at each order of 1/N . All summation symbols run from i = 1 to i = (2p − 1)!!.
to compute the high temperature expansion of the free energy to a certain order in β and 1/N , we only need to compute a finite number of moments to an appropriate order in 1/N . In particular, we have found in the thermodynamic limit the 1/q 3 correction to the coefficient of β 6 (there are no other large q corrections) and both the 1/q 3 and 1/q 4 corrections to the coefficients of β 8 (these are the only large q corrections). The leading order 1/q 2 results as well as the 1/q 3 results are in agreement with a large N calculation of the path integral formulation of the SYK model, while the 1/q 4 correction was not calculated before. This relation also allows for calculations of moments to higher order in 1/N , and we have pushed the calculation to order 1/N 3 . More surprisingly, we found that to a given order in 1/N , all moments are determined by a finite number of moments. This also explains the pdependence of the momentsM 2p obtained in a previous work [1]. One important consequence of this is that the SYK model generates elegant enumeration identities, at each consecutive order in 1/N , as discussed in section 4. We tabulate some of them in table 3. To the best of our knowledge, it seems that only the identities at order 1/N are present in mathematical literature [61], and our study of the SYK model suggests they are only the first layer of a hierarchy of identities.
It is clear that the procedure we presented can be extended to even higher orders in 1/N , and the most computationally burdensome part is to calculate and expand η G for the irreducible structures, for which we have a general and practical formula [1]. For the large N expansions of generic systems, there is often no powerful simplifying principle that allows the calculation to be pushed to higher and higher order easily. For SYK model, although the large N coefficients of the moments are still somewhat complicated, the enumeration identities generated at each order of 1/N (table 3) are incredibly simple. Does this simplicity of enumeration persist to higher and higher orders? If it does, does it imply there is something we can say about the moments to all orders in 1/N instead of calculating them order by order? We hope to clarify these questions in future studies. Another issue is that we are in a somewhat awkward situation that on the one hand we have the Q-Hermite expression for the spectral density [1], which is a resummed approximation that is very accurate and gives the leading order free energy for all temperatures, but Q-Hermite moments are only 1/N -exact; on the other hand, we have expressions for the moments that are exact to order 1/N 3 , but they only give the high temperature information of the free energy. It would be desirable to find an improved resummed expression for the spectral density that goes beyond the Q-Hermite approximation, which gives a free energy that is both 1/N 3 -exact and accurate at all temperatures.

Acknowledgments
Mario Kieburg is thanked for pointing out the analytical result for q = 1. Gerald Dunne is acknowledged for discussions on large order expansions and Antonio García-García is thanked for a critical reading of the manuscript. Y.J. and J.V. acknowledge partial support from U.S. DOE Grant No. DE-FAG-88FR40388. After submission of the paper, an interesting preprint appeared [62] that also uses chord diagrams to calculate a class of moments of the SYK Hamiltonian, but otherwise has no significant overlap with the present work.
A High temperature expansion of free energy from [56] In section 3.2.3 we claimed that our result for the free energy for even q reproduces an independent calculation [56] to order 1/q 3 , and in this appendix we justify this claim.
In [56] the large q expansion for the extensive (leading in 1/N ) part of the free energy is given: If we set J 2 = 2 q−1 /q as we did in section 3.2.3, we have simply Then by iterating eq. (A.3) we obtain Using this high temperature expansion of u, we obtain for the free energy, eq. (A.1), This is consistent with our result eq. (3.26) to order 1/q 3 at leading order in 1/N .
The advantage of the result of [56] is that it is valid at all temperatures, however our result, albeit only valid at high temperatures, gives the expansion to even higher orders in 1/q and 1/N .
As a side remark, the coefficients of the high temperature expansion of u are known analytically. If we write u = 1 π ∞ n=0 (−1) n (2n + 1)!2 2n a n β 2n+1 , (A.6) then a n are given by a n = 1 2 2n+1 This is a convergent series with a convergence radius of β c ≈ 1.33. The numbers a n are also the number of labeled rooted trees on 2n + 1 nodes with each node having an even number of children [63]. Since there is no singularity on the positive real axis, we do not expect a phase transition as β varies.
B.1 q = 1 case As discussed in detail in [1], for q = 1, H 2p = N α=1 J 2 α p 1, and because J α is Gaussian distributed, N α=1 J 2 α p follows a χ 2 distribution. Hence for the q = 1 SYK model we havẽ where again we have only kept the 1/N 3 term.
To evaluateM QH , we first expand η E to 1/N 3 . For q = 1 the first equation of (4.

B.2 q = 2 case
As explained in detail in [1], for q = 2 we havẽ where the brackets · · · on the right-hand side denote the ensemble average with the probability distribution [43,64], We can work out the values of the W i 's by employing a set of recursion relations for Selberg integrals developed in [65] resulting in To relevant order, we have (B.14) For q = 2 we finally arrive at where we have omitted all but the 1/N 3 term. The results up to order 1/N 2 can be found in [1]. To evaluateM QH 2p , the 1/N expansion of η in the first equation of (4.6) simplifies for q = 2 to