A Generalization of Sachdev-Ye-Kitaev

The SYK model: fermions with a $q$-body, Gaussian-random, all-to-all interaction, is the first of a fascinating new class of solvable large $N$ models. We generalize SYK to include $f$ flavors of fermions, each occupying $N_a$ sites and appearing with a $q_a$ order in the interaction. Like SYK, this entire class of models generically has an infrared fixed point. We compute the infrared dimensions of the fermions, and the spectrum of singlet bilinear operators. We show that there is always a dimension-two operator in the spectrum, which implies that, like in SYK, there is breaking of conformal invariance and maximal chaos in the infrared four-point function of the generalized model. After a disorder average, the generalized model has a global $O(N_1) \times O(N_2) \times \ldots\times O(N_f)$ symmetry: a subgroup of the $O(N)$ symmetry of SYK; thereby giving a richer spectrum. We also elucidate aspects of the large $q$ limit and the OPE, and solve $q=2$ SYK at finite $N$.


Introduction
The Sachdev-Ye model [1], as recently revived and simplified by Kitaev [2], possesses, for large N , three remarkable properties: conformal invariance in the infrared, solvability, and maximal chaos. While there are models that contain some of these properties, SYK is the first to have all three, as was recognized by Kitaev in a series of incredibly insightful seminars [2].
Broadly speaking, until recently two classes of large N theories have been studied: matrix models and vector models, in which the dynamical variables transform in the adjoint or fundamental representation of a local or global SU (N ) or O(N ) symmetry group, respectively.
Matrix models are closely related to string theories [3][4][5][6][7][8], with the most concrete realization being the duality between supersymmetric gauge theories and string theory in Anti-de Sitter space [9]. N = 4 super Yang-Mills is conformally invariant, and at large 't Hooft coupling the bulk gravity has black holes, so it should be maximally chaotic. However, it is not easily solvable. Vector models also have a long history and recently have been shown to be dual to interesting gravity theories. The critical O(N ) vector model is conformally invariant and solvable, and the bulk dual is higher spin Vasiliev theory [10,11]. However, it is integrable for large N , so it is not likely to be chaotic. Roughly speaking, matrix models are too difficult to be explicitly solvable, while vector models are too simple to have the same rich properties.
One would like a model that lies in between: one that is sufficiently complicated to be chaotic, while still simple enough to allow for direct analytic calculations for strong coupling. SYK is such a model. At large N , the dominant Feynman diagrams for matrix models are planar diagrams, whereas the dominant diagrams for vector models are bubble diagrams. The SYK model is dominated by a new class of Feynman diagrams, which have been referred to as sunset, or watermelon, diagrams. The SYK model may be just one example out of a much broader and new class of models. Past studies of large N models have been extremely fruitful for understanding both quantum field theories and string theories. One may hope that the study of SYK-like models will also prove productive.
SYK is a quantum mechanics model, living in 0 + 1 dimensions. While two-dimensional CFTs have been extensively studied and categorized, one dimensional CFTs have not. In fact, it has been argued that 0 + 1 dimensional CFTs with nontrivial dynamics do not actually exist [12]. SYK confirms this: the four-point function breaks the SL(2, R) conformal invariance [2,13,14], consistent with holographic studies of AdS 2 [15]. It appears that in one dimension a theory can at best only be "nearly" conformally invariant. In SYK the breaking of conformal invariance, to leading order in 1/N , is confined to a single dimension-two operator appearing in the OPE, so the power of conformal invariance is still largely applicable.
The SYK model consists of N 1 Majorana fermions χ i , with a q-body Hamiltonian with quenched disorder, The model has qualitatively similar properties for any choice of even q ≥ 4. The couplings J i 1 ,...,i q are independently chosen from a Gaussian, O(N ) invariant, distribution with zero mean and a variance proportional to J 2 N 1−q . When evaluating observables, say correlation functions, a disorder average is performed at the end of the calculation. For the purposes of correlation functions, at large N , the model is self-averaging for q > 2: randomly chosen, but fixed, J i 1 ,...,i q give the same results as disorder averaged J i 1 ,...,i q . One can alternatively think of the J i 1 ,...,i q as nearly static free bosonic fields; at leading order in 1/N , this gives the same connected correlation functions [16], and furthermore, allows one to gauge the O(N ) symmetry [17]. where a tower of massless modes appears in the bulk. In [14] it was noted that the bulk dual of SYK might be a string theory with the string scale comparable to the AdS radius, and thus non-local or stringy. But what the dual of SYK is, and the extent to which it is nonlocal, remains an open problem.
The goal of this paper is to generalize the SYK model. We would like to understand how large the class of such models is, and which features are generic and which are special to SYK. This paper will not add anything new to the bulk interpretation of SYK, but the dual bulk theory, whatever it is, should be able to incorporate this more general class of models.
Two seemingly important ingredients in SYK are: (a) 0 + 1 dimensions, where the fermions are dimensionless, thereby ensuring that any product of fermions is a relevant perturbation, and (b) quenched disorder, which plays an important role in the solvability at large N . The generalization we explore is one in which there are f flavors of fermions, χ a i , where i = 1 . . . N a and a = 1 . . . f , with a Hamiltonian, We prove that for any choice of parameters: f , N a 's, q a 's, there is always a dimension-two operator in the spectrum. In SYK, the dimension-two operator is responsible for both the breaking of conformal symmetry in the four-point function and for maximal chaos. The same properties hold for the generalized model.
An instructive case to study is the generalized model with all N a equal to N/f and all q a equal to q. It has the symmetry O(N/f ) × · · · × O(N/f ). The spectrum contains a tower identical to that of SYK with a qf body interaction, along with a new tower that appears with a degeneracy of f − 1. Indeed, this model is similar to SYK with N fermions and a qf body interaction, but the full O(N ) symmetry is broken and consequently more singlet operators exist, allowing for a richer model.
In Appendix A we consider the path integral for the generalized model. This provides an alternate way of computing the correlation functions, with the saddle point giving the Schwinger-Dyson equations for the two-point functions, and the leading 1/N fluctuations about the saddle giving the four-point function. In Appendix B we consider (1.2) with an additional scalar; a special case of this includes supersymmetric SYK [19]. Finally, Appendix C solves SYK for q = 2 at finite N . SYK for q = 2 is like N fermions with a random mass matrix; the randomness makes it nontrivial, though it is less interesting than q ≥ 4. This appendix can be read independently of the rest of the paper.

SYK
Let us recall the SYK model [2]. 1 It contains N Majorana fermions with the anticommutation relation {χ i , χ j } = δ ij . The action is, where the coupling J i 1 ,...,i q is totally antisymmetric and, for each i 1 , . . . , i q , is chosen from a Gaussian ensemble. The two-point function of the J i 1 ,...,i q is taken to be, (2.2) At leading order in 1/N , (2.2) is equivalent to the simpler normalization, The particular scaling with N in the choice (2.3) is in order to obtain a nontrivial large N limit, while the other factors are for convenience. One can consider SYK for any even q ≥ 2, with q = 4 being the prototypical case [2]. At q = ∞ there are some simplifications [14].
The case q = 2 is simplest, and is equivalent to an O(N ) vector fermion with a random mass matrix, although in many ways it is qualitatively different from the q > 2 models. We solve the q = 2 SYK at finite N in Appendix C. At zero coupling, the Euclidean two-point function T χ i (τ )χ j (0) ≡ G(τ )δ ij is given by, where the factor sgn(τ ) (sgn(τ ) = 1 for τ > 0 and sgn(τ ) = −1 for τ < 0) accounts for the fermion anticommutation. To leading order in 1/N , the Schwinger-Dyson equations for the two-point function drastically simplify and are given by, The first of these, (2.5), is the standard equation expressing the two-point function in terms of the one-particle irreducible self-energy Σ(ω). The second equation, which is written in position space, is a special feature of SYK (see Fig. 1). At leading order in 1/N , the only diagrams that survive are nested sunset diagrams; all others are suppressed by some power of 1/N . These equations can be combined into a single integral equation; however, an analytic solution to this equation is not known. At strong coupling, |Jτ | 1 (equivalently, the infrared limit), one can drop the iω in (2.5), to get, One can verify that, is a solution to (2.7) provided one takes, The Fourier transform of G(τ ), given in (2.8), is useful in verifying this, where we defined, . (2.11) What is special to SYK is that the IR Schwinger-Dyson equations (2.7) are invariant under reparameterization of time, τ → f (τ ), with the propagator transforming as G(

A Generalization of SYK
The model we introduce is a generalization of SYK (2.1). It contains f flavors of fermions, with N a fermions of flavor a, each appearing q a times in the interaction, so that the Hamiltonian couples q = f a=1 q a fermions together. We continue to let the subscript on the fermion χ a i denote the site i ∈ {1, . . . , N a }, while the superscript a will now denote the flavor a ∈ {1, . . . , f }. Explicitly, the action is, where I is a collective index, I = i 1 , . . . , i q 1 , . . . , j 1 , . . . , j q f . The coupling J I is antisymmetric under permutation of indices within any one of the f families, and is drawn from a Gaussian distribution, where the disorder average J I J I is given by (q a − 1)! .
(2.14) It will be convenient to make the following definitions, The class of models (2.12) for large N is characterized by f − 1 independent continuous parameters 0 < κ k < 1, as well as the q k , which can be any positive integers provided that In the SYK model (2.1), one can generalize the action to have multiple interaction terms, with different q, each coming with its own independent disorder J i 1 ...i q . This sum of SYK Hamiltonians is just as solvable as SYK, however it is not especially interesting since in the IR the term with smallest q will be dominant. In the model (2.12), one can also consider generalizing the action to include sums of interaction terms. However, now the IR can be more interesting, since there can be multiple terms (with the same total q) that are equally important in the IR.

Two-Point Function
The free two-point function for each χ a i , is again given by (2.4). Away from the UV it will continue to be the case that the two-point function is diagonal in flavor and site space.
Denoting the two-point function for the flavor k fermion by G k (τ ), the self-energy for the flavor k fermion is (see Fig. 2 Making use of (2.14, 2.15), this simplifies to, We first determine whether there is an IR fixed point and, if so, the IR dimension ∆ k of the fermions of flavor k. In the IR, the two-point function should take the form, To find the normalization b k and dimension ∆ k , we first insert the above ansatz into (2.17) and take the Fourier transform, 3 Inserting (2.19) and (2.18) into the IR limit of (2.5), Σ k (ω)G k (ω) = −1, gives, The first equation is just the statement that the IR dimension of the coupling J I is zero. 2 The factors are as follows. The factor ( q a !) 2 comes from the square of the prefactor of the interaction term in (2.12). There is a factor of q k from the number of contractions with the ingoing fermion, and another q k with the outgoing fermion, and a (q k − 1)! from the contraction of the remaining flavor k fermions amongst themselves. There is also a factor of q a ! from contractions of flavor a fermions, for all the other flavors. 3 One should not confuse the usage of Σ as the self-energy with the usage of Σ as a sum.
Simplifying the second gives, Equating all the (2.22), for k ranging from 1 to f , gives f − 1 equations. Combined with (2.20), for any given choices of κ k and q k , we have a set of f equations for the f unknown dimensions ∆ k . 4 These equations have simple solutions in the limit of q a 1, as we show in the next section.

Large q k
If the number of fermions of flavor k appearing in the interaction (2.12) is large, q k 1, then from (2.20) we know that ∆ k 1. Let us assume q k 1 for all k. In this limit, (2.22) simplifies to b q a a = 1 2 κ k Q k ∆ k , with the solution, (2.23) Eq. 2.23 shows that a hierarchy in the q k 's for different flavors, or in the κ k 's, will lead to a hierarchy in the ∆ k 's.
The smallness of the dimensions ∆ a suggests one should be able to solve for the twopoint function at all energies. This was done for SYK at large q in [14]. Here we perform an analogous computation for the generalized model in the large q k limit. The two-point functions are taken to be, Taking the Fourier transform, for which we use the shorthand F, Inverting to get G k (ω) −1 , (2.5) allows us to identify, We have such an equation for every k. Thus, we can express g a in terms of g k for any a, k, Summing (2.27) for all k and using (2.28) gives, where the rescaled J is kept finite in the large q k limit. The solution to (2.29) is easily derived for finite temperature, namely with g a (τ ) = g a (τ + β) [14], where v is defined implicitly in terms of J . At zero temperature (2.30) becomes, which combined with (2.28) gives, where ∆ k is given by (2.23). Having the exact solution, we can take the IR limit J|τ | 1 to find the individual normalizations of the two-point function (2.18),

Graphical Solution
The two-point function in the large q limit can alternatively be found by summing an appropriate set of Feynman diagrams. We will show how this works in SYK. Due to large q combinatorics, the diagrams contributing to the self-energy that appear most often are like the ones shown in Fig. 3 (a), rather than those in Fig. 3 (b). The Feynman diagrams that are summed at large q can be characterized as those diagrams that, under a single vertical cut, break up into two tree diagrams. The self-energy can therefore be found recursively, as shown in Fig. 3(c). The equation corresponding to Fig. 3(c) is, where G 0 (ω) is the free two-point function (2.4). Rearranging (2.34) gives where the second equation is the inverse Fourier transform of the first. Letting we get, which is (2.27) for one flavor.

Spectral Function
The frequency space two-point function follows from (2.24, 2.32), Introducing a Schwinger parameter, and performing the τ integral in Eq. 2.38 gives, The spectral function (as defined in Appendix C by Eq. C.4) for the flavor k fermion is therefore, where ∆ k is given by (2.23). Since ∆ k 1, this is sharply peaked around small λ. If there is only one flavor, then ∆ = 1/q. This spectral function is for q 1. For q = 2, the SYK spectral function is instead a Wigner semicircle (C.5).

Effective Action
We have so far discussed the model (2.12) directly in terms of the fermions, finding the two-point function at large N through study of Feynman diagrams. It is useful to also consider the path integral approach. Employing the replica trick, one can carry out the disorder average, and then integrate out the fermions after the introduction of new (bilocal) For one flavor, this reduces to the effective action for SYK [2] (see [1] for an analogous expression for the SY model, and [18] for the Dirac fermion version of SYK). The large N saddle point of the action gives the Schwinger-Dyson equations for the two-point function found previously from Feynman diagrams. In particular, varying S ef f with respect to G k (τ 1 , τ 2 ), and assuming time-invariance, gives (2.17), while varying with respect to Σ k (τ 1 , τ 2 ) yields (2.5) for each flavor. The saddle point solutions are denoted by G k (τ 1 , τ 2 ) and Σ k (τ 1 , τ 2 ).
To leading order in 1/N , the free energy is given by the saddle of (2.42), Following [14], one can differentiate with respect to J to get, For large q a , G a (τ ) was found explicitly in Sec. 2.3. Also, since the partition function only depends on βJ, it follows that J∂ J = β∂ β . Thus for large q a , where v is defined in terms of J in (2.30). Up to the choice of normalization of the variance of the disorder, J I J I , this is the same as for SYK with N fermions and a f a=1 q a body interaction. So the entropies are also the same. In order to see a distinction, one must study the 1/N corrections.

Four-Point Function
The SYK model has an O(N ) symmetry after the disorder average. The bilinear primary In the UV, these have dimension 2n + 1. The IR dimensions of the operators are computed by summing a class of ladder diagrams. The four-point function of the fermions is then given by a sum over conformal blocks, one for each of these composite operators.

Dimensions of Composite Operators
We begin by reviewing and adding some detail to the computation in [2] for the IR dimensions of the SYK composite operators. The primary O(N ) invariant bilinear operators are, where the coefficients d nk are chosen so that the operators are primary. For instance, The general form of d nk will not be important for us.
We would like to compute the overlap between the state created by the composite operator O n acting at time τ 0 , and two fermions at times τ 1 and τ 2 , respectively. In other words, the three-point function, , which we will denote by v(τ 0 ; τ 1 , τ 2 ). If the two fermions just propagated without interacting with each other, this would be found by Wick contractions, Eq. 3.3 is the first diagram that appears in Fig. 4a. We must also include a sum over all the ladder diagrams in Fig. 4a. One can perform the sum by solving the equation (see where the kernel is the operator that adds a single rung (see Fig. 4c), Letting the composite have dimension h, in the IR the solution to (3.4) will take the form of conformal three-point function, For h > 2∆, the term G 0 χχO is much smaller than (3.6) in the IR, τ 12 1, so we can drop it in (3.4). Thus, (3.4) simplifies to (see Fig. 4d), take the eigenvectors to be, where ψ(∆) was defined in (2.11). A plot of g(h) is given in Fig. 5 for AdS 2 . So the dual of SYK has a tower of particles in the bulk, one for each solution to g(h) = 1 for h > 2∆. For large integer n, these have approximate masses m 2 n ≈ (2∆ + 2n +  Figure 6: The diagonal (a) and off-diagonal (b) components of the kernel (3.13, 3.14) for two flavors with q 1 = q 2 = 3. The coloring scheme is the same as in Fig. 2.

Generalized Model
We now generalize the calculation to the model with flavor (2.12). The operators (3.1) now have a superscript O a n to account for the different flavors, and the kernel is now a matrix in flavor space, K mn (τ 1 , τ 2 , τ 3 , τ 4 ), where m denotes the flavor of the incoming fermions on the left at times τ 1 , τ 2 and n denotes the flavor of the outgoing fermions on the right at times τ 3 , τ 4 , see Fig. 6. The off-diagonal component K kl has flavor k propagators along the rails, while the rung consists of q k − 1 flavor k propagators, q l − 1 flavor l propagators, and q a flavor a propagators for all a = k, l, where the combinatorial factor in front is, 5 The diagonal components of the kernel are similar, but with slightly different propagator powers and combinatorial factors. Using (2.14, 2.15) and simplifying we get the diagonal and 5 The factor of N q a a , for a = l, k, comes from the site index summation within the rung. For flavors k, l, there are only q k − 1, q l − 1 propagators in the rung, so those give factors of N There is then an additional factor of N l because the Feynman diagrams are built by adding the kernel to the left (see Fig. 4); so the l index will get summed over.
As in SYK, we must find the eigenvectors and eigenvalues of the kernel. Letting g be an eigenvalue, and v a (τ 12 ) the components of an eigenvector, (3.15) Following (3.8), an ansatz for an eigenvector is, v a (τ 12 ) = c a sgn(τ 12 ) with some coefficients c a . Since the eigenvector and propagators in the kernel only depend on time differences, (3.15) factorizes nicely under a Fourier transform, where the first term on the left is from the diagonal term in the kernel (3.13) and the second term is from the off-diagonal terms (3.14). Inserting the propagator (2.18) and evaluating gives, In order to eliminate the dependance on ω, we must choose, In SYK we know that (3.8) is a special case of (3.6) with 2α = 2∆ − h. Similarly here, we let 2α a = 2∆ a − h , (3.20) which is consistent with (3.19). Thus, (3.18) becomes an eigenvector equation for the matrix K, where the diagonal and off diagonal components of K are, . (3.23) In getting from (3.18) to (3.22) we made use of the product of normalizations of the propagators b q a a given in (2.21). If there is one flavor, K 11 reduces to (3.9). The next step is to find all h for which there is an eigenvalue g of K (3.22) that equals 1. This is in principle straightforward: for any q a , κ a in (2.12) ones solves (2.20, 2.22) to find the IR dimensions ∆ a of the fermions, then for fixed h one finds the f eigenvalues of K, and then for each of those eigenvalues solves for all h such that the eigenvalue is equal to one.
Aside from some special cases, we can not write a general and explicit answer for the h's.
However, it is easy to see that there will always be a dimension 2 operator in the spectrum.
For h = 2 (3.23) simplifies to, Inserting this into (3.22), one can easily verify that the following vector is an eigenvector of K with eigenvalue one. Verifying this requires using q a ∆ a = 1, and nothing else. Perhaps surprisingly, it is not even required that the ∆ a are actual dimensions: such an eigenvector has two nonzero components, In fact, verifying this requires no assumptions on ∆ a . The presence of these dimension-one operators suggests a symmetry. In fact, this symmetry is simple to see from the effective action (2.42). 6 One can rescale q a G a (τ 1 , τ 2 ), for any a = 1, while leaving the IR limit of (2.42) invariant.

Equal q a , κ a
A simple and instructive case is when all the q a are equal to some q, and all the κ a are equal, for all flavors a. The dimensions ∆ a are then, by symmetry, all equal to ∆ = 1 f q . The matrix K in (3.22) factorizes, where ρ(h) is given by (3.23) and is independent of the flavor. The eigenvalues of K are thus, where where n k ranges from 1 to N 1 , and gets rid of the flavor index on the fermions. There are now f N 1 sites; however, this is not the same as SYK with a qf body interaction, since the interactions are not all-to-all, being restricted to occur between particular qf sets of sites.

Four-Point Function
Having found the dimensions of the bilinear singlet operators O a n , the next step is to compute their OPE coefficient. The OPE between two fermions will include all the O a n and their descendants, and will take the form, where c a,b n are the OPE coefficients and C n (τ 12 , ∂ τ 1 ) = 1 + . . . is fixed by conformal invariance. The OPE coefficient can be extracted by computing the three-point function between the two fermions and O, which in Sec. 3.1 was labelled as v(τ 0 ; τ 1 , τ 2 ) and satisfied Eq. 3.4. Thinking of K as a matrix with indices (τ 1 , τ 2 ), (τ 3 , τ 4 ), the formal solution of (3.4) is, v a (τ 0 ; τ 1 , τ 2 ) = 1  it was unimportant in the IR. However, for finding the OPE coefficients one is interested in the UV, τ 12 1, so this term is essential.
We will also be interested in the four-point function. Defining the bilocal, and proceeding formally, one can perform a double OPE expansion on the four-point function, 7 g a (τ 1 , τ 2 )g b (τ 3 , τ 4 ) = 1 N a N b n,e c a,e n c b,e n C n (τ 12 , ∂ τ 1 )C n (τ 34 , ∂ τ 3 ) 1 |τ 13 | 2h n,e . (3.33) The right-hand side is a sum of conformal blocks, given by hypergeometric functions of the conformally invariant cross ratio, n,e c a,e n c b,e n x h n,e 2 F 1 (h n,e , h n,e , 2h n,e , x) , x = τ 12 τ 34 τ 13 τ 24 .

(3.34)
This is similar to two-dimensional CFTs, except here we have one cross-ratio instead of two.
At large N , the leading and first subleading in 1/N pieces of the four-point function are, The 1/N piece of the four-point function is found by summing ladder diagrams [2]. What we have is a slight generalization of what occurs in SYK, as the four-point function is now a matrix in flavor space. Starting with F ab 0 (τ 1 , τ 2 , τ 3 , τ 4 ) = δ ab (−G a (τ 13 )G a (τ 24 ) + G a (τ 14 )G a (τ 23 )) , (3.36) one uses the kernel (3.13, 3.14) to add rungs to the ladder. Summing all the ladder diagrams, OPE coefficients c n [14], where g(h) is given by (3.9) and h n are the solutions of g(h n ) = 1.
Eq. 3.38 is for h n > 2. There is an additional complication that occurs for the h = 2 block. One can notice that g(h = 2) = 1, and since h = 2 is part of the basis of eigenvectors used to invert 1 − K, this causes the four-point function to diverge in the conformal limit.
The h = 2 block must therefore be treated outside the conformal limit. Moving slightly away from the IR, the eigenvalue g(h = 2) gets slightly shifted away from 1, and so the h = 2 block gives a finite but large, and non-conformal, contribution to the four-point function. Since its prefactor is dominant, its growth controls the behavior of the finite temperature out-oftime-order four-point function used to probe chaos [39,40]. The growth of the h = 2 block occurs with a Lyapunov exponent 2πT that saturates the chaos bound of [41]. 8 In Sec. 3.1.1 we found that the generalized model always contains a dimension-two operator; assuming its OPE coefficient doesn't vanish, this implies that in the IR the generalized model, like SYK, both breaks conformal invariance and is maximally chaotic.
To compute the four-point function (3.37) for the generalized model with generic q a and κ a , one would need to repeat the procedure used for SYK, accounting for the additional complexity of having flavor. However, in the case that all the q a are equal and all the κ a are equal, it is simple to find the four-point function, and this is case we focus on.

Equal q a , κ a
If all the q a are equal to q, and all the κ a = 1/f , then the kernel matrix (3.13, 3.14) factorizes into a flavor-space matrix and a function of times, where K was defined in (3.27). By symmetry, the two-point functions are flavor-independent, we diagonalize (3.37) in flavor space, forming O T FO, to find, where, Here α 0 (q) is given by (3.38) and ρ(h) is given by (3.23) with fermion dimension ∆ = 1/2q.
To be clear, the h n appearing in (3.45) and (3.46) are the solutions of (2q − 1)ρ(h n ) = 1 and −ρ(h n ) = 1, respectively. Recall that we found in Sec. 3.1.2 that with two flavors (with q 1 = q 2 = q, κ 1 = κ 2 ), the spectrum of bilinear composite operators contains two towers: a tower that matches the 2q body SYK tower, and a new tower, see Fig. 7. The OPE coefficients c S n are for the 2q body SYK tower. Note that (3.45) is for h n > 2; as discussed before, the contribution of the h = 2 block diverges in the conformal limit.
The OPE coefficients c A n are for the new tower. Notice that this vanishes for the h = 1 operator. The OPE coefficients c 1 n and c 2 n , in terms of (3.45, 3.46), are given by, where, for simplicity of presentation, rather than writing c a,b n , we have explicitly separated the two towers. 9 A more intuitive way to think about this four-point function is to define the symmetric and antisymmetric combinations of the bilocals (3.32), The symmetric correlator probes only the SYK tower, and matches the SYK 2q body four-point function. The antisymmetric correlator,

Discussion
The SY model [1] involves all-to-all interactions between spins in some representation M µ,ν=1 J ij S µ i ν S ν j µ , with Gaussian-random couplings J ij . Writing the spins as products of two fermions, this becomes a four-fermion interaction. One of the key realizations of [1] was that the model is solvable in the double scaling limit, N → ∞, M → ∞, M/N → 0. It was recognized in [2] that a simpler model is one that avoids spins altogether and goes directly to the fermions, H = J ijkl χ i χ j χ k χ l . There is then a clear generalization to a model with a q-index coupling and a q-body interaction [2].
In this paper, we have made another straightforward generalization, involving f flavors of fermions with N a sites for each flavor and a f a=1 q a body interaction. Perhaps surprisingly, the model has an infrared fixed point for most choices of parameters N a , q a . We found a set of equations determining the dimensions of the fermions in the infrared, as well as the matrix determining the infrared dimensions of the bilinear singlet operators that are invariant under It was recognized in [14] that the SYK model simplifies in the limit q 1. Here we pointed out that in the q 1 limit, only a particular subset of Feynman diagrams need to be summed. For any even q ≥ 4, the SYK model has qualitatively similar properties. In the generalized model introduced in this paper, there are more parameters to vary, and one may wonder if there are corners of parameter space which either lead to simplifications or qualitative differences.
We have only begun exploring the parameter space, focusing on the symmetric case of an equal number of sites for each flavor, as well as interaction orders q a that are independent of the flavor. The main qualitative difference, as compared to SYK with a qf body interaction, is more singlet operators resulting from a symmetry that is a subgroup of the O(N ) symmetry of SYK. One feature we found, that holds for any choice of parameters, is the presence of a dimension-two bilinear singlet operator in the infrared. Another was the presence of a dimension-one operator; however, in the symmetric case considered, its OPE coefficient vanished. It would be good to better understand this operator.
Nontrivial and solvable models are both rare and valuable. It is now clear that the class of SYK-like models is much larger than just the SY model. Just how large this class is, if there are further generalizations, and the precise characterization of the Feynman diagrams, at each order in 1/N , are all still open problems. We may hope that exploring this structure will provide guidance towards understanding the dual string theory, if there is one.

Acknowledgements
We thank D. Anninos

A. Effective Action
In this appendix we compute the free energy (equivalently, the effective action) for the generalized model (2.12). The calculation is analogous to the one for SYK [2].
Employing the replica trick, instead of computing the disorder average of the logarithm of the partition function, one instead computes the disorder average of M copies of the system. Starting with (2.12), this is given by, Having done the disorder average, we see that there is a O( symmetry. We thus introduce the collective fields, by inserting delta functions, where Σ αβ a (τ 1 , τ 2 ) acts as a Lagrange multiplier. We insert into (A.2) such a delta function for each replica index pair α, β and each flavor a. This gives, Integrating out the fermions gives where As is standard in studies of SYK, one assumes a replica symmetric saddle point, G αβ a (τ 1 , τ 2 ) = δ αβ G a (τ 1 , τ 2 ), and so (A.6) becomes If there is only one flavor, we recover the SYK action [2],
Now consider the generalized model, (A.8), for two flavors with q 1 = q 2 = q and κ 1 = By symmetry, the saddle point is, As noted in Sec. 2.4, the saddle point equation is the same as the equation for SYK with a 2q body interaction. Let us now study fluctuations about the saddle, G a = G + |G| 1−q g a , and Σ a = Σ + |G| q−1 σ a . Expanding (A.13) to second order, Integrating out σ 1 , σ 2 gives, where g S /g A are the symmetric/antisymmetric combinations of g 1 , g 2 , (3.49). The correlators

B. Model with a Scalar
In this appendix we consider a model with a boson field. It is a slight variant of (2.12) and has the action, where I is a collective site index I = i, i 1 , . . . , i q 2 , . . . , j 1 , . . . , j q f , and q 1 = 1, and q = f a=2 q a . This has a similar interaction as (2.12), but the first flavor is with a boson instead of a fermion.
The boson is taken to be auxiliary, having UV dimension [φ] = 1/2, so the coupling J I has dimension 1/2. We restrict to only one boson, q 1 = 1, in order to ensure that the interaction is relevant. We also require q be even. 10 The only technical distinction between finding the IR dimensions for (B.1) compared with (2.12) is that the boson propagator is symmetric in time. The ansatz for the IR boson propagator is, where ψ(∆) is defined in (2.11), while for the IR fermion propagator it is, for k ≥ 2. In the IR, one drops the free propagator appearing in the Schwinger-Dyson equation (2.5), so for both the boson and the fermions one has Σ k (ω)G k (ω) = −1. The self-energy is again given by (2.16), with the disorder average normalization given in (2.14).
Repeating the several steps in Sec. 2.2 gives the following equations, The For q 2 = 2 this gives the dimensions ∆ 1 = 2/3, ∆ 2 = 1/6 found in [43,44]. 11 Intriguingly, the difference between the boson dimension ∆ 1 and fermion dimension ∆ 2 is 1/2 for any q 2 , as would be implied by supersymmetry. However, we have not checked that the model for general q 2 is supersymmetric.

C. Random Mass Matrix Fermions
The simplest SYK model is for q = 2: fermions with a random mass matrix. For this case, all computations can be performed exactly, without the restriction of being near the fixed points or working at large N . In this appendix, we solve the q = 2 model. For infinite N this is trivial, while for finite N it is slightly more involved but follows from standard 11 We thank Yu Nakayama for sharing his results and explaining the supersymmetric model.
One can also find this directly by summing non-crossing rainbow diagrams (see Fig. 9), where G 0 (ω) is the bare propagator (2.4), while C n are the Catalan numbers, The Catalan numbers are the number of different ways n + 1 factors can be completely parenthesized; here the parentheses are the rainbows. Summing (C.2) gives (C.1). One can also write (C.1) as, where the spectral function ρ(λ) is the Wigner semi-circle, with support for |λ| < 2J, (C.5) 12 In [45] it was proposed that the q = 2 model satisfies the Eigenstate Thermalization Hypothesis.
At finite temperature, the frequencies in (C.1) should be viewed as the Matsubara frequencies, ω n = (2n + 1)π/β. Taking the discrete Fourier transform of (C.4) gives, where 0 < τ < β. In the limit of zero temperature βJ 1, we can evaluate the integral to obtain, where L 1 is the modified Struve function, and I 1 is the modified Bessel function. While both L 1 and I 1 grow exponentially, the two-point function (C.7) decays monotonically with |Jτ |.
The combination is sometimes denoted by M 1 ≡ L 1 − I 1 . We can do a strong coupling expansion of (C.7), where we see that the first term matches what was expected from the IR limit of the Schwinger-Dyson equations (2.8, 2.9).

Comments
One comment is that in summing the Feynman diagrams giving (C.2) it is important to work at finite temperature. Each of the diagrams individually has IR divergences: the Fourier transform of any of the individual terms in (C.2) will diverge in the limit of β → ∞. Of course, one could have chosen to regulate the IR divergence in some way other than working at finite temperature. However, finite temperature is natural. The point is just that the dimensionless expansion parameter is βJ.
Another comment is that aside from the implicit appearance of β in the Matsubara frequencies, (C.2) has no explicit β dependance. This is a property that is special to q = 2.
For SYK with q ≥ 4, one could solve the Schwinger-Dyson equations perturbatively around weak coupling, giving an expansion of the form, (ω n β) 2l , (C.9) with some coefficients g kl . One can derive recursion relations for g kl , but we have not found a way of solving them.
We have been considering Majorana fermions. One can instead study q = 2 with Dirac fermions, At leading order in 1/N , this gives the same two-point function (C.1). Working with Dirac fermions gives slightly more flexibility, as one can introduce a chemical potential. With no chemical potential, as in (C.10), one is at half-filling. Explicitly, consider a single free Dirac fermion H = ω 0 c † c. (Adding a chemical potential just corresponds to adding to (C.10) such a term for each fermion, with chemical potential µ = −ω 0 .) The finite-temperature two-point function is trivially, Since we are at finite temperature, fields have the time range 0 < τ < β. The two-point function is a function of the difference between two times, and so naturally has the range −β < τ < β. However, from (C.11) we see that for 0 < τ < β, G(τ − β) = −G(τ ). We can thus restrict to 0 < τ < β. The filling fraction Q is defined as the expectation value of the occupation number, We can choose the filling fraction by choosing ω 0 . It is clear that for any finite temperature, if ω 0 = 0, then there is no energy cost to being in the state |1 versus |0 , and so the filling fraction is 1/2. Note that the limits of T → 0 and ω 0 → 0 do not commute. If we set T = 0 at finite ω 0 (including ω 0 = 0), then we get zero filling: from (C.11), G(τ ) = θ(τ ).
Finally, the q = 2 model sums rainbow diagrams. There are many other models that sum rainbow diagrams. For instance, two-dimensional QCD [46] has the same rainbow diagrams, where the fermions are the quarks, and the disorder lines are the gauge field propagators.
Also, the recently studied three-dimensional U (N ) k Chern-Simons theory coupled to scalars or fermions also sums rainbow-like diagrams [47,48]. A simple large N quantum mechanics model that sums rainbow diagrams is the Iizuka-Polchinski model [49] (see also, [50,16]). The IP model has a harmonic oscillator in the adjoint representation of U (N ) plus a harmonic oscillator in the fundamental representation of U (N ), coupled through a trilinear interaction.
In the limit that the mass of the adjoint goes to zero, this is essentially the same as the model Eq. C.10, at leading order in 1/N . The reason we say essentially the same is because in the IP model the fundamental is effectively at zero filling. In other words, its free twopoint function is θ(τ ) as opposed to 1 2 sgn(τ ), and correspondingly, the infinite N two-point function after summing the rainbow diagrams is only the first term, I 1 , in (C.7). In addition, at subleading orders in 1/N , differences will arise between the model Eq. C.10 and the IP model. This is because the adjoint propagator will receive quantum corrections, whereas the J ij J ij "propagator" in Eq. C.10 is always a constant.

Four-Point Function
It is simplest to write the four-point function in frequency space. This is defined as, Written as a series in 1/N , ijkl + . . . . At leading order in 1/N there is a disconnected piece, At first subleading order in 1/N , the four-point function is a sum of ladder diagrams, like SYK for general q. However for q = 2 there is an extreme simplification, since the rungs only contain the disorder lines. Since there is no momentum exchange, summing the ladders in frequency space simply involves summing a geometric series, which gives, It is only necessary to establish the first term in the s-channel piece. The other term, as well as the t and u channels, follow from antisymmetry. Through a Fourier transform and analytic continuation of (C.14), one finds there is no exponential growth in the out-of-timeorder four-point function [16], and so the q = 2 model is not chaotic.

C.2. Finite N C.2.1. Dirac Fermion
We now compute the two-point function at finite N for the random mass matrix fermion, (C.10). For fixed coupling J ij , this is just N free fermions with mass matrix J ij , so the nontrivial part is to perform the disorder average. Specifically, Consider first the trivial case of N = 1. This is just a fermion with a random mass. Then (C.15) reduces to (C.4) with a spectral function, 2 .
(C. 17) So at N = 1 the spectral function is a Gaussian, while at N = ∞ it is the Wigner semicircle (C.5). The two-point function for N = 1 can also be found by summing Feynman diagrams, Writing the double factorial as a Gaussian integral, and interchanging the integral and the sum, we recover (C.17). Explicitly performing the integral gives the two-point function in terms of the complimentary error function, In the zero temperature limit we also find, We now move on to the case of general N , using the method of orthogonal polynomials to evaluate (C.15). We can write (C. 16 We now take linear combinations of the columns of (C.22), transforming it into a matrix with i, j element, φ j (λ i ), where φ j (λ i ) is a polynomial with lowest element 1 and highest element λ j i . The determinant (C.22) remains invariant under these operations. We can write the determinant as a sum of permutations of the integers from 0 to N − 1, (C.23) We choose the φ n such that, dλ φ n (λ)φ m (λ) e −λ 2 /2J 2 = f n δ nm . (C. 24) This then gives, The φ n will be proportional to the Hermite polynomials, defined as, and so we find that the spectral function is, Evaluating the sum gives, where we have used that f k = J 2k+1 √ 2π k!, which follows from (C.26, C.27). A plot of (C.30) is shown in Fig. 10.
An alternative way to write the two-point function is to perform the integral over λ in (C.4) before evaluating the sum over k appearing in the spectral function (C.29). After the introduction of a Schwinger parameter, the integration over λ yields a Laguerre polynomial. Using that the sum of the Laguerre polynomials is an associated Laguerre polynomial .
The leading term in 1/N , g (0) , reproduces what we found from summing the planar diagrams, (C.1).

C.2.2. Majorana Fermion
Here we compute the two-point function for the Majorana version of q = 2 SYK (2.1) at finite N (note that N must be even). This will be slightly different from the Dirac version studied in Sec. C.2. The two-point function is given by, The procedure is now similar to the Dirac case. We can write the determinant as a sum of permutations of the even integers from 0 to N − 2, The φ n are the same as in the Dirac case. The partition function now involves just the even normalization constants, For evaluating the two-point function, note that since the eigenvalues come in pairs,