Sparse Sachdev-Ye-Kitaev model, quantum chaos and gravity duals

We study a sparse Sachdev-Ye-Kitaev (SYK) model with $N$ Majoranas where only $\sim k N$ independent matrix elements are non-zero. We identify a minimum $k \gtrsim 1$ for quantum chaos to occur by a level statistics analysis. The spectral density in this region, and for a larger $k$, is still given by the Schwarzian prediction of the dense SYK model, though with renormalized parameters. Similar results are obtained for a beyond linear scaling with $N$ of the number of non-zero matrix elements. This is a strong indication that this is the minimum connectivity for the sparse SYK model to still have a quantum gravity dual. We also find an intriguing exact relation between the leading correction to moments of the spectral density due to sparsity and the leading $1/d$ correction of Parisi's U(1) lattice gauge theory in a $d$ dimensional hypercube. In the $k \to 1$ limit, different disorder realizations of the sparse SYK model show emergent random matrix statistics that for fixed $N$ can be in any universality class of the ten-fold way. The agreement with random matrix statistics is restricted to short range correlations, no more than a few level spacings, in particular in the tail of the spectrum. In addition, emergent discrete global symmetries in most of the disorder realizations for $k$ slightly below one give rise to $2^m$-fold degenerate spectra, with $m$ being a positive integer. For $k =3/4$, we observe a large number of such emergent global symmetries with a maximum $2^8$-fold degenerate spectra for $N = 26$.

employed [5] to model quantum chaotic dynamics in a many-body context and also certain aspects of quantum magnetism [6].
More recently [7][8][9][10], a variant of these models based on N Majoranas [7], the so called Sachdev-Ye-Kitaev (SYK) model, has attracted a lot of attention as a toy model for holography, and for its potential to reveal novel insights in the dynamics of strongly interacting quantum matter. In the low temperature (strong coupling) limit, the SYK model shares the same pattern of soft breaking of conformal symmetry [11] by finite temperature and quantum (1/N) effects as that of Jackiw-Teitelboim (JT) gravity [12,13], a two-dimensional gravity theory with a dilaton in Anti-de Sitter space with non-trivial boundary conditions. This symmetry breaking pattern dictates low temperature thermodynamic properties [7,8,[14][15][16] such as a linear specific heat, and an exponential growth of low energy excitations. These are all expected features in field theories with a black hole gravity dual.
Another distinctive feature of these systems is quantum chaos [17]. Quantum chaos reveals itself in level statistics described by random matrix theory [18] and also in the exponential growth at the scrambling time of quantum corrections measured by certain out-of-time-order correlation functions, with a growth rate controlled by the Lyapunov exponent. Kitaev [7] found that this feature occurs in the SYK model and that, in the strong-coupling, low temperature limit, the Lyapunov exponent saturates a previously proposed universal bound on chaos [17]. Regarding level statistics, both the SYK model [16,[19][20][21][22] and JT gravity [23,24] are well described by random matrix theory [25][26][27][28][29] which indicates that the system is quantum chaotic at all times scales.
A natural question to ask is how the above features, which determine the existence of a quantum black hole dual, are robust to deformations of the SYK model. Typically, generalizations of the SYK model to higher spatial dimensions [30] or involving more Majoranas than the usual four-body interaction [31] share similar features. However, the addition of an integrable two-body interaction [32][33][34] prevents the saturation of the Lyapunov exponent. Moreover, in a certain range of parameters, the system is not quantum chaotic as spectral correlations are well described by Poisson statistics, typical of an integrable system.
Another example of a generalized SYK model in which quantum chaos may not occur is that of a two-site coupled SYK model that in the low temperature limit is dual [35] to an eternal traversable wormhole. It was shown in Ref. [36] that the traversable wormhole phase is not quantum chaotic. Quantum chaotic features are only observed for higher temperatures where the gravity dual undergoes a thermodynamic transition to a quantum two-black-hole background.
Another plausible deformation of the SYK model is to relax the requirement of infinite range interactions. Indeed, in the context of condensed matter physics [37], interacting quantum dots describing realistic electronic interactions are qualitatively similar to the SYK model with complex fermions but with a Fock space geometry living on a Cayley tree rather than on a complete graph. The effect in the SYK model of a sharp-cut off in Fock space distances [38] induces a metal-insulator transition. However, not much is known about the requirements on the range or the form of the interactions that guarantees the existence of a gravity dual. Progress on this problem would not only bring a more detailed understanding on the conditions for a field theory to have a gravity dual but also it might be useful to identify systems to test experimentally holography predictions.
Here we study the properties of a sparse SYK model where some of the couplings are randomly set to zero with a probability 1 − p, where p ∼ k/N α , α > 0 and k is a positive real number. This model was first articulated in a talk given by Brian Swingle [39]. Our main aim is to characterize the maximum sparseness for which both the spectral density and spectral correlations are consistent with that of a gravity dual. Namely, the spectral density is still described by the Schwarzian prediction of the dense SYK model [7,8] and level statistics are still quantum chaotic [16,20] and therefore well modeled by random matrix theory. For that purpose, we have computed analytically the spectral density, and the partition function, by an explicit calculation of the moments of the Hamiltonian. We show that it still has a Schwarzian form, and therefore it is likely related to a gravity dual, provided that α ≤ 3 and, for α = 3, k ∼ 1 or larger. A study of spectral correlations confirms agreement with random matrix theory in this region of parameters which indicates the dynamics is still quantum chaotic at late time scales. For k = 1 extra symmetries and chiral symmetries emerge for some disorder realizations, and level statistics of the three Wigner-Dyson ensembles and the three chiral ensembles are observed for an ensemble of 26 Majorana fermions. For k < 1 we find a large number of emergent discrete symmetries as well as chiral symmetries leading to exact degeneracies in powers of 2.
We note that formally the SYK Hamiltonian is defined over random hypergraphs. As we shall see, there are not many mathematically rigorous results for generic random hypergraphs as a function of the degree of sparseness. The situation is different in the simpler case of random sparse graphs, usually termed Erdos-Renyi graphs [40] which can be cast as L × L matrices. There is a rather rigorous characterization [41][42][43][44] of the bulk spectral properties: level statistics consistent with the prediction of random matrix theory will occur if the fraction of nonzero matrix elements satisfies p ≥ L /L with > 0. In this region, the spectral density is given by the semicircle law [45,46]. These results are fully consistent with numerical [47] and analytical results [48,49] in the physics literature. Close to the edge of the spectrum, the spectral region related to the gravity dual, it was demonstrated rigorously [42] that, for p ≥ L /L 2/3 , spectral correlations are described by RMT. As far as we know, it is unclear whether this bound is optimal.
These findings cannot be directly applied to the sparse SYK model as its Hamiltonian is not represented by a graph but by a more complex random sparse hypergraph for which not many explicit results for the density or spectral correlations are available. An exception [50] is the spectral density of a SYK-like model that can be cast as a √ N-hypergraph. We refer to [51] and references therein for recent mathematical results about the conditions to observe the semicircle law in random hypergraphs. We are not aware of any level statistics characterization of random hypergraph in the mathematical literature. In the physics literature, we refer to Ref. [52,53] for an analytical calculation of the two-level correlation function in a fermionic model with infinite range interactions.
The paper is organized as follows: in Section II, we introduce the sparse SYK model, the mechanism to tune the sparseness and the regularity condition that makes the connectivity on the hypergraph uniform for each disorder realization. Section III is devoted to an analytical evaluation of the spectral density as a function of the degree of sparseness. We notice a striking equivalence between leading corrections to the moments of the density, due to the sparsity of the model, and leading 1/d corrections of the same quantity in the Parisi's U(1) gauge model on a d-dimensional hypercube [54]. Based on exact analytical results for low order moments, we propose two approximate analytical expressions for moments of any order. For one of them, we write down a closed analytical expression for the spectral density. In Section IV, these predictions are compared to numerical results resulting from the exact diagonalization of the sparse SYK Hamiltonian. In Section V, we turn to the study of the conditions for the existence of quantum chaos by an analysis of spectral correlations in both the bulk and the edge of the spectrum as a function of the degree of sparseness. Section VI is focused on the description of emergent global symmetries that only occur in the limit of strong sparsity. For a fixed number of Majoranas, these additional symmetries, that depend on the disorder realization, lead to spectral degeneracies and spectral correlations described by random matrix ensembles of different universality classes including those with chiral symmetry. Finally, in Section VII, we summarize the main results of the paper and list some problems for future research. The numerical implementation of the regularity condition is discussed in Appendix A. In Appendix B, we discuss two examples of emergent symmetries.

II. SPARSE SYK MODEL
We investigate the following Hamiltonian representing N strongly interacting Majorana fermions [7] with sparse q-body infinite range interactions. For q = 4, where we have used Euclidean Dirac matrices γ i to represent Majorana fermions. Dirac matrices satisfy the anti-commutation relations They are the same as the anti-commutation relations of Majorana fermions up to a factor of two, which will be absorbed in the definition of the variance of J i jkl in Eq. (3). The sparseness is implemented by the random variable x i jkl : x i jkl = 1 with probability p and x i jkl = 0 with probability 1 − p. We may think of the interactions to be defined on random hypergraphs: i = 1, 2, . . . , N labels the N vertices. A hyperedge connecting the i, j, k, l vertices is present if x i jkl = 1. The expected number of the hyperedges (and hence the number of terms in the Hamiltonian) is p N 4 . The hypergraph gets sparser as p gets closer to zero, and p = 1 gives the maximal number of hyperedges, N 4 , which results in the conventional "dense" SYK model. The couplings J i jkl is a Gaussian random variable with the distribution where J sets the scale of the distribution. We will set J = 1 for later numerical calculations. We focus on a probability p that scales as p ∼ N −3 for which it is convenient to define an order-one quantity k: We shall study other scalings p ∼ N −α , but unless stated explicitly, we set α = 3.
Although the numerical results of this paper will be restricted to q = 4, certain analytical results will also be available for other integer values of q > 0. For general q, we write the Hamiltonian as where α is a multi-index with q elements: so that α can take N q different values. The random coupling J α follows a Gaussian distribution and Γ α is the product of q Dirac matrices indexed by α: The hyperedge variable x α is defined analogously, the probability p now scales as N 1−q and When k is small, some of the random hypergraphs are disconnected. The disconnectedness makes the Hamiltonian split into a sum of sub-Hamiltonians defined on independent tensor subspaces. In this case, the spectral statistics become a superposition of statistics from different sectors possibly belonging to different symmetry classes. To mitigate this complication, we can consider regular hypergraphs only, that is, we can impose that every vertex has the same degree kq, namely each vertex is contained in the same number kq of hyperedges. This regularity condition is imposed as a set of constraints on x α : where α m means that, among the N q choices of α, we sum over those α that contain a given integer m. The value of k must be chosen so that kq is a positive integer. This regularity condition implies that every realization of the Hamiltonian contains exactly kN = p N q number of independent non-zero terms, as opposed to the case without the regularity condition where kN is only the number of nonzero terms on average. For ordinary random graphs, where each edge connects two vertices, graphs are almost surely connected provided they are regular and its vertex degree is larger than 2 [40,55,56]. For random regular hypergraphs, since a hyperedge connects more than two vertices, we expect connectivity to be more easily achievable and hence the vertex degree need not be as large. This is indeed true for any vertex degree kq > 1 for which any random regular hypergraphs will be almost surely connected [57]. However, as we will see, even with regularity condition, in the very sparse regime 1/q < k ≤ 1, spectral statistics can still be a superposition of independent spectra because there are new emergent global symmetries.

III. SPECTRAL DENSITY: ANALYTICAL RESULTS
We evaluate analytically the spectral density by an explicit computation of the moments of the sparse SYK Hamiltonian, The spectral density can be expressed as Since we have a Gaussian distribution of J α (in the notation of Eq. (5)), the calculation of the average requires to consider all possible Wick contractions. In the end, we will also need to average over the random variable x α .
After averaging over J α , the result depends on whether pairs of two factors Γ α are adjacent or not. In the former case we can use that while in the latter case, the Γ α 's can be made adjacent by using [16], where c αβ = |α ∩ β| is the number of indices that α and β have in common. An exact calculation of a generic trace requires us to keep track of correlations with other factors Γ. This is in general a challenging combinatorial problem but some low-order moments have been evaluated exactly [16,58] for the dense SYK model. The simplest Wick contraction in which Eq. (14) plays a role is out of which we will define an order-one quantity In general, a Wick contraction that contributes to M 2l is a trace of a product of 2l matrices Γ, whose subscripts form l pairs that are summed over. Repeatedly using Eqs. (13) and (14), we can move all the pairs of Γ's with the same subscripts next to each other and produce a purely combinatorial expression of the form where c α i α j = |α i ∩ α j | and crossings includes all those pairs of i, j for which α i and α j form a "crossing" configuration in the trace: The binomial factor in front of Eq. (17) is to normalize the sum to an order-one quantity, which can be alternatively understood as normalizing the moments M 2l to reduced moments M 2l /M l 2 as we shall see soon. There is an intersection graph representation [58] of the quantity defined in Eq. (17): Two examples of intersections graphs G for the dense SYK model. The left graph represents N q −2 α 1 α 2 (−1) c α 1 α 2 which is the η defined in Eq. (16); the right graph represents N q −3 2. If any c α i α j is in the crossings c α i α j , connect the vertices α i and α j by an edge.
We can now rewrite Eq. (17) in terms of intersection graphs: a Wick contraction represented by an intersection graph G that contributes to the 2l-th reduced moment has a value of [59] where is an edge of the graph G connecting the vertices α i and α j and c α i α j is the number of common indices in α i and α j . We give two examples of intersection graphs and what they represent in Fig. 1.
In this notation, the reduced moments for the dense SYK model can be written as where G are all the l-vertex intersection graphs representing Wick contractions. An important approximation to the dense SYK model moments is the so-called Q-Hermite approximation [19,20,50,58]: where η is defined in Eq. (16) and E(G) is the number of edges in G. Under this approximation, the dense SYK moments are where the second equality is the Riordan-Touchard formula [60,61]. The moments of the Riordan-Touchard formula are exactly the moments of the spectral density for the Q-Hermite polynomials [62]: where is the ground state energy. This is the reason why we called Eq. (21) the Q-Hermite approximation for the (dense) SYK moments.
When it comes to averaging over x α , the regularity condition matters. Without regularity condition, a generic averaging x α 1 x α 2 · · · x α l can be worked out by simply noting that x 2 α = x α , x α = p and that two x-variables are statistically independent if they have different subscripts. With regularity condition, two x-variables can be correlated even if they have different subscripts (essentially because, given the regularity constraint, the various x-variables are not extracted independently from each other), which makes the combinatorial problem much more difficult. For this reason, in this paper analytical results are only available for the model without regularity condition, and the regular model will only be studied numerically. However, since the regularity condition in Eq. (10) implements N constraints on N q otherwise independent variables, we expect that the regularity condition only modifies the moments by contributions of order 1/N q−1 which are subleading with respect ot those considered below.
Without regularity condition, the second moment is given by, which is the same as for the dense SYK model. To calculate the fourth moment we need to be careful with the average over x α variables because We write down the result in terms of the reduced fourth moment without regularity condition, where M 4,SYK is the fourth moment of the dense SYK model and η is given by Eq. (16). Similarly, the result for the sixth moment without regularity condition is, Likewise, the full expression for the eighth moment without regularity condition is, (29) where T 6 is the value for the triangle intersection graph (see Fig. 1): which in the Q-Hermite approximation Eq. (21) is given by T 6 ≈ η 3 . In principle, it is possible to compute higher order moments but the final expression becomes increasingly cumbersome. It is clear from the explicit calculations so far, and from a general proof to be given soon, that the leading moment is always that corresponding to the dense SYK model. Hence, in the large N limit, we already know the behavior of sparse SYK moments: they are the same as the large N limit of the dense SYK model. Two commonly taken large N limits are: 1. Fixed q and N → ∞: in this limit the global spectral density approaches a Gaussian.
2. Fixed q 2 /N and N → ∞: in this limit the global spectral density approaches the density function of the Q-Hermite polynomials with Q = η which for large N can be approximated as Q → exp −2q 2 /N [50].
However, we would like to understand how the large N limit is approached or, in other words, we would like to understand the sparse SYK model at large but finite N with q fixed. Moreover, we would also like to understand the form of the low energy excitations slightly above the ground state for the q = 4 model, which is not captured by the two above-mentioned global limits. The form of the spectral density in this infrared region is relevant to the type of gravitational theory the sparse SYK model might be dual to.
We now set out to study the finite N behavior of the sparse SYK moments. As a start, we would like to draw the readers' attention to the subleading terms of the moments.
(14 + 28 cos φ + 28 cos 2 φ + 20 cos 3 φ + 10 cos 4 φ + 4 cos 5 φ + cos 6 φ) Note we will slightly abuse the terms "leading" and "subleading" for the Parisi's model: the natural parameter for the large d expansion of M 2l /M l 2 is not powers of 1/d but d(d − 1) · · · (d − m + 1)/d l instead. Hence for the 2l-th moment, "leading" means m = l and "subleading" means m = l − 1. We see that if we apply Q-Hermite approximation to both the leading and the subleading moments of the sparse SYK model, and make the identification of η = cos φ, then the leading moments of the sparse SYK model Eq. at least based on the observation of the first eight moments of both models.
We will see now why Eq. (33) is true not only for the first eight but for all moments. In the sparse SYK model, an intersection graph G represents a value of where the notation η G x serves to distinguish it from its dense SYK counterpart η G , and reminds us of the fact that there is an extra averaging over x variables in the sparse SYK model. We shall distinguish the two cases in the actual drawings of the intersection graphs by annotating the vertices of the sparse SYK intersection graphs by (α i , x α i ), as opposed to by α i alone for the dense SYK intersection graphs defined earlier. See Fig. 2 for an example. If all the subscripts in Eq. (34) are different, x α 1 x α 2 · · · x α l will be equal to p l ; if two of the subscripts become equal, x α 1 x α 2 · · · x α l will be enhanced by a factor of 1/p, but the restriction on the summation will suppress the sum by N q , and the total effect is a 1/(kN) suppression. It is clear then that at leading order in 1/(kN), Eq. (34) is given by coinciding with the dense SYK model value η G . This proves that to leading order, the moments of the sparse SYK are exactly the same as those of the dense SYK. Hence, in the Q-Hermite approximation, they are given by the Q-Hermite moments Eq. (22). In the case of the Parisi model, it is already understood that the leading contribution is given by the Q-Hermite prediction [65]. Therefore, to leading order, Eq. (33) is proven, namely, the moments of both the dense SYK model (after Q-Hermite approximation) and the Parisi model are given by the Q-Hermite prediction. This is perhaps not too surprising, but we will see that the subleading correction in Eq. (33) arises in a much more subtle and surprising way.
The subleading order of Eq. (34) can be written as where there are l 2 sums corresponding to letting two out of the l subscripts be equal. One might worry about excluding the cases where even more indices are equal, but they are of higher order and do not enter into our consideration here. We can summarize the above results as For example, Fig. 2 gives an intersection graph G that contributes the sixth moment, and in the form of Eq.
From this example, it is clear that there is a natural graphical representation of the subleading calculations: 2. An intersection graph example of Eq. (37). Note the left-hand side and the first term on the right-hand side have identical graph but different labeling of the vertices, and hence the left represents η G x whereas the first term on the right represents η G .
1. Merge two vertices α i , α j of the intersection graph G, and let the merged vertices inherit the original edges.
2. There may be loops (an edge connecting a vertex back to itself, representing (−1) q ) and 2-multiedges (two edges connecting the same pair of vertices, representing (−1) 2c α i α j ) formed after step 1, delete all such loops and 2-multi-edges. Call the resulting graph G (α i ,α j ) .
3. The subleading contribution to η G x is given by where v(G) denotes the vertex set of G, and η G (α i ,α j ) is the value for which the intersection graph G (α i ,α j ) would represent a dense SYK model. When the Q-Hermite approximation is applied, this where E G (α i ,α j ) denotes the number of edges in the graph G (α i ,α j ) .
The "merge and delete" procedures applied in the example of Eq. (38). We show how to obtain the subleadingmoment intersection graphs from the leading-moment intersection graph in Fig. 2.There are three subleading graphs corresponding to the merging of α 1 α 2 , α 1 α 3 and α 2 α 3 , we have drawn two of them because merging α 1 α 2 and merging α 2 α 3 result in identical graphs. Fig. 3 illustrates an example of the application of these rules. The above "merge and delete" graphical rules to calculate the subleading moments, which result in Eq. (40), are exactly the same as the "averaged scheme" defined in [64] for calculating the subleading moments of the Parisi's hypercube model, except that the averaged scheme for the Parisi's model has an extra factor of 1/3. We can now conclude that Eq.
(33) holds for all moments. Such coincidence does not hold to the next order in 1/(kN), as is evident by comparing the eighth moments of the sparse SYK and the Parisi's model at higher orders.

B. The renormalized and subleading Q-Hermite approximations
The leading intersection graphs G introduced in the previous section can be summed by the Riordan-Touchard formula [60,61] after applying the Q-Hermite approximation, η G ≈ η E(G) , to both the leading and subleading terms, We will call Eq. (41) the subleading Q-Hermite approximation. One can easily check that only including the leading term results in a fairly large discrepancy with the exact result so the subleading term is an important contribution. We would like to get a grasp of how accurate the subleading Q-Hermite approximation is. In we are not yet able to find a Riordan-Touchard-like formula for the subleading terms. A related difficulty is then that we are not able to write down an analytic expression for the subleading spectral density. At this point we can simply remark that it would be worthwhile to find a Riordan-Touchard-like formula since it would solve the subleading problem of two models at one stroke. This difficulty prompts us to try a different strategy of approximating moments. We will only use the leading moment expression, but with a renormalized η, with the hope it can capture subleading effects beyond it natural range of applicability of relatively low moments. More specifically, calling this renormalized parameter η(k), the moments are given by, We dub this approximation the renormalized Q-Hermite approximation.
The lowest moment in which η starts to make an appearance is M 4 /M 2 2 . The renormalized Q-Hermite approximation Eq. We remark that this renormalized Q-Hermite approximation already fails to fully capture the subleading term of the sixth moment (except at η = 1). However, this approximation can be justified a posteriori: we shall see it is surprisingly close to the exact moments for certain ranges of N and k and to the resulting spectral density as well.
In order to gain a more quantitative understanding of the suitability of these approximations, we compare the subleading Q-Hermite, see Fig. 4, and the renormalized Q-Hermite approximation, see Fig. 5, with exact results for the sixth and eighth moments of the q = 4 model Eq. (1). We have observed that: 1. In the very sparse limit, k = 1, the subleading Q-Hermite approximation is the better approximation for N 60.
2. When N is relatively small, N 30, the subleading Q-Hermite approximation is the better approximation.
3. For larger k, such as k ≥ 3, the accuracy of the renormalized Q-Hermite approximation starts to catch up with that of the subleading Q-Hermite approximation, and rather surprisingly at first glance, beyond N = 40 its accuracy exceeds that of the subleading Q-Hermite approximation. This can be understood partly from the observation that for η = 1 the renormalization cancels the 1/(kN) terms exactly in case of the sixth and eighth moments. In fact, it can be shown [66], that this observation for η = 1 is true for all moments due to a result for edge counting of intersecting graphs. We will see in the next section that this results in a surprisingly good agreement between the renormalized Q-Hermite prediction and the numerical spectral density.
Finally, we note that the moments in Eq. (42) give rise to the same spectral density as in Eq. (23) but with η replaced by its renormalized version η(k): where c N is a normalization constant, and is the ground state energy.
We compare in next section this analytical prediction with the numerical spectral density from the exact diagonalization of the sparse SYK Hamiltonian Eq. (1).

C. Conditions for the existence of a gravity dual
A distinctive feature of the existence of a gravity dual in the context of the SYK model is that, for E sufficiently close to the ground state E 0 , the spectral density becomes, with γ a non-universal constant directly related to η. This is the result of the exact quantum path integral computation of the classical Schwarzian action [67] which is 1/N exact. The classical Schwarzian captures the soft breaking from conformal to S L(2, R) symmetry that characterizes both, the infrared limit of the SYK model and certain near AdS 2 backgrounds [7,8,11]. These symmetry considerations are enough to determine the effective low energy theory that is then quantized.
The analytical moment calculation that we have carried out indicates that, for the sparse SYK model with α < 3, corrections due to the sparsity of the Hamiltonian are subleading with respect to 1/N corrections which strongly suggests that the spectral density is still given by Eq. (46) and therefore it could still have a gravity dual. The case α = 3 is more interesting. The leading correction due to the sparsity of the Hamiltonian is of order 1/kN and therefore it modifies the expansion leading to the Q-Hermite approximation in the dense SYK. However, the analytical moment calculation earlier in this section, together with the comparison of the renormalized Q-Hermite approximation with numerical results supports that Eq. (44) provides a good description of the spectral density of the model for large but finite N and even relatively small k provided that 1/kN is small. In principle, this means that the expression for the spectral density Eq. (46) is still valid with γ = γ(k). This will be shown explicitly in Fig. 10, but it may be argued that we had to remove by hand the strong fluctuations of E 0 in order to clearly observe the edge of the spectrum which casts some doubts on the applicability of Eq. (46) and indirectly on the existence of a gravity dual.
We think that these concerns are unfounded. The fluctuations in E 0 are a direct consequence of the quantization procedure we have followed. Instead of picking up the classical low energy effective theory and then quantizing the gravitational degrees of freedom of interest, we are quantizing the full theory without suppressing other degrees of freedom which leads to strongly enhanced fluctuations. Moreover, the collective excitations that induce fluctuations in E 0 are also 1/N suppressed so we expect them to become a smaller problem if larger N could be explored numerically. Therefore, we believe that removing the fluctuations of E 0 , a degree of freedom of no direct interest in our analysis, is an approximation in line with that of first identifying the effective low energy classical action and then proceeding with the quantization [67]. This is specially true when we have strong evidence that the renormalized Q-Hermite approach provides a very good description of the spectral density in the bulk of the spectrum. Right: For the critical scaling α = 3, we observe similar features for k < 1. For k > 1, the density is qualitatively similar to the dense SYK model.
It is a quite exciting prospect that even a strongly sparse SYK model could have a gravity dual. If so, it may be possible to push this idea further and investigate specific conditions on the geometry of the Fock space which could be favorable to the existence of a gravity dual. More specifically, it may be possible to establish the minimum requirements on connectivity so that the spectral density has black-hole like features such a stretched exponential form, ∝ e a √ E−E 0 with a independent of energy, in the infrared limit.

IV. SPECTRAL DENSITY: NUMERICAL RESULTS
We compute the eigenvalues of the Hamiltonian Eq. (1) by exact diagonalization techniques. The resulting spectral density is very sensitive to the probability p ∼ k/N α . For α > 3, and a small value of k, we observe, see Fig. 6, a depletion of eigenvalues around to E = 0, and an increase of statistical fluctuations.
For α < 3 and k > 1, we expect the spectral density to be similar to that of the dense SYK model. We    The analytical result given by the black curve has the right width but its average position disagrees with the numerical results.
The situation gets better for larger k. In Fig. 9 we show the distribution of the first 10 eigenvalues (red curves) for N = 26 and k = 1 with (left) and without (right) regularity condition as well as the analytical result for the distribution of the smallest eigenvalues (black curve). The width of the distribution is much larger than the level spacing, and the distributions of the first 10 eigenvalues are almost identical. The analytical and numerical are clearly closer, and it is clear that a large number of small eigenvalues contribute to the tail of the spectral density.
These results are a strong indication that α = 3 and k ∼ 1 is the maximum degree of sparseness, or the connectivity in Fock space, that can support the existence of a gravity dual.
Despite the good agreement, we observe visible differences in the infrared part of the spectrum. The numerical result has a smooth tail while the renormalized Q-Hermite density predicts an edge. The reason behind the numerical tail is the strong fluctuations of E 0 for different disorder realizations. It is well known that disorder induces collective excitations in the spectrum, which blur the existence of spectral edges.
In order to study the tail in more detail we remove these collective excitations by dividing all eigenvalues of each realization by its largest eigenvalue. The spectral density of these renormalized eigenvalues for an ensemble of 1000 realization with N = 26 and k = 4 is shown in Fig. 10. It is also shown the Q-Hermite spectral density with fitted values for η = 0.129 and E 0 = 1.008. The fitted value of η is considerably less that the theoretical value of 0.164 (without the 1/(kN) correction) or 0.193 (with the 1/(kN) correction).
One might argue that η should be given by the value corresponding to the internal fourth moment Eq. (55), but it is actually quite a bit smaller. In the right panel for Fig. 10, we depict a magnification of the tail of the spectral density. There is an agreement with the Q-Hermite spectral density (red curve) almost to the square root edge.

A. Scale Fluctuations
For sparse matrices, the number of independent stochastic variables defining the Hamiltonian is kN.
Therefore the relative error in an observable is 1/ √ kN. If we decompose each eigenvalue into the ensemble average and a small deviation, we thus have that This corresponds to scale fluctuations of the eigenvalues. It is natural to introduce a stochastic variable ξ, that describes the scale fluctuations of the spectrum over different disorder realizations. The scale fluctuations follow from the variance of the second moment: For the contribution of scale fluctuations to the reduced moment we find On the other hand, we can evaluate the above moment exactly through Wick contractions and explicit trace calculation, and the exact result is, This results in This means that the Thouless scale is only The O(1/N 2 ) correction also includes the 1/N q contribution from the dense SYK model [64].
In Fig. 10, where the subscript "int" (internal) refers to the fourth moment where the contributions of the scale fluctuations have been eliminated. This give the internal fourth reduced moment where the last term is sub-leading. Indeed this gives a reduced value of η, but the fitted value of η is still considerably smaller.
It is straightforward to numerically calculate δE 2 i for an ensemble of sparse SYK Hamiltonians. In Fig. 11 we show δE i ≡ δE 2 i 1/2 versus the ensemble average E i of the i-th eigenvalue for N = 32 and various values of k. In particular, for larger values of k there is a linear dependence on E i confirming the above analysis. The slope of the curves versus 1/k is given in the right panel of Fig. 11 and compared to Eq. (53) (red solid curve). Except for the point at k = 3/4, the agreement is excellent.

V. SPECTRAL STATISTICS AND QUANTUM CHAOS
We now study the late time dynamics associated with the Hamiltonian Eq. (1) by a level statistics analysis. Spectral correlations are a valuable probe to describe the quantum dynamics for long time scales of the order of the inverse mean level spacing. Agreement with RMT signals that the dynamics is quantum chaotic while Poisson statistics corresponds to an insulator or an integrable system [68]. The bulk of the spectrum corresponds to the high temperature phase while the low temperature/strongly coupled region is related to the lowest eigenvalues of the spectrum. In principle, only the latter is related to the existence of a gravity dual.
In order to proceed, we obtain the spectrum of the model by exact diagonalization techniques. Since the matrix representation of the Hamiltonian is extremely sparse, the use of Lanczos's algorithm allows to reach up to N = 42 Majoranas. As already discussed, for sufficiently small k, it is useful to impose the regularity condition Eq. (10) so that all Majoranas live on a connected hypergraph. We discuss in Appendix A an efficient method for the numerical implementation of the regularity condition.
Except for the calculation of the form factor, the procedure of spectral unfolding is carried out by relatively low order < 5 polynomials.
Since our main goal is to establish the maximum sparseness consistent with quantum chaos, we will be mostly interested in short-range spectral correlators, such as the level spacing distribution, P(s), and the adjacent gap ratio. The former is defined as the probability to find two consecutive eigenvalues E i , E i+1 at a distance s = (E i+1 − E i )/∆ (with ∆ the average local level spacing). For a fully quantum chaotic system it is given by Wigner-Dyson statistics [69] which is well approximated by the so-called Wigner surmise that depends on the universality classes [68]. For the Gaussian Orthogonal Ensemble (GOE), Gaussian Unitary Ensemble (GUE), Gaussian Symplectic Ensemble (GSE) is given by: P W,β (s) = a β s β exp b β s β with β = 1, 2, 4, respectively. a β , b β are numerical coefficients [68].
For an insulator, or a generic integrable system, it is given by Poisson statistics, P P (s) = e −s . The adjacent gap ratio is defined as [70][71][72], for the ordered spectrum For a Poisson distribution, it is equal to r P ≈ 0.38 while for a random matrix ensemble it depends on the symmetry class, with r ≈ 0.53, 0.60, 0.67 for the GOE, GUE, GSE [73], respectively. The advantage of r over P(s) is that it does not require us to unfold the spectrum. For that reason, we will also consider the full distribution of the adjacent gap ratio ρ(r). An analytical Wigner-surmise for ρ(r) is available for different random matrix ensembles [73], with β = 1, 2, 4 for GOE, GUE and GSE respectively. The prefactor A β is a numerical coefficient and [73] (note the difference with Eq. (56)). We note that despite P(s) and ρ(r) are both short-range spectral correlators, that probe the quantum dynamics at times of the order and larger than the Heisenberg time, ρ(r) is a shorter range correlator than P(s). Therefore, we expect that deviations from RMT predictions will become more apparent in P(s).
We start the spectral analysis with the study of spectral correlations near the center of spectrum, usually called the bulk, corresponding to the high temperature phase.

A. Bulk
We define the bulk as the central part of the spectrum comprising 80% of eigenvalues unless other percentage is explicitly stated. Our first task is to determine the critical scaling p N 4 = kN 4−α for which the dynamics is quantum chaotic, namely, level statistics are well described by RMT. For that purpose, we first compute P(s) defined above for N = 26, k = 2 and different scalings of the probability p parameterized by α. The results depicted in Fig. 12 strongly suggest that the maximum sparseness consistent with quantum chaotic dynamics is approximately p ∝ 1/N 3 , namely, α = 3. This is in agreement with the prediction for Erdos-Renyi graphs adapted to random hypergraph represented by the Hamiltonian Eq. (1). We note that for α > 3, not only the tail is exponential, as for Poisson statistics, but also there is a peak for small s related to spectral degeneracies that we shall see soon are related to the presence of emergent global symmetries for sufficiently strong sparseness.
In order to see that α = 3 corresponds to the maximum sparseness, called from now on the critical scaling, we first study the level statistics for a larger k = 4 and different N's. The global symmetries of the SYK model depend on N [16,21], so a study of the N dependence in the sparse case will also provide useful information about the robustness of these symmetries against the sparsing procedure. We have found that the agreement with the RMT results corresponding to the different universality classes (GOE, GUE, GSE) is excellent, see Fig. 12. Moreover, the results for N = 26 and N = 30, both belonging to GUE, are almost indistinguishable. Both features provide convincing evidence that in the region k 1 and α = 3 the system is still fully quantum chaotic with not much difference with the dense case at least for short range correlations of few neighboring eigenvalues.
We note that this robustness of quantum chaos is remarkable. The dense SYK model has ∼ N 4 non-zero different entries while the sparse one only 4N. This is however the analytical prediction resulting from a heuristic extrapolation of the rigorous mathematical results [41,42], and numerical simulations [47], for random sparse graphs to hypergraphs such as the sparse SYK model: the dynamics is quantum chaotic and spectral correlation are described by RMT only for sufficiently large k.
We now turn to the study of the dependence of spectral correlations on k for this critical scaling (p ∼ k/N 3 ) to determine the minimum k = k c for which this agreement to RMT persists. The theoretical expectation for random graphs [47] is that k c 1. In the previous investigation of level statistics for α > 3, we have noticed the emergence of level degeneracies at least for some disorder realizations.
Qualitatively, the reason is that the quantum dynamics is very sensitive to the overall connectivity of the hypergraph. Therefore, for sufficiently small k, or large α, level statistics strongly depend on the connectivity of the disorder realization. As an example, for sufficiently small k ≤ 2, in some cases, we observe double degeneracy while in others realizations, the spectrum has a chiral symmetry E → −E. For some disorder realizations, both a double degeneracy and a chiral symmetry occur at the same time. We will study this phenomenon in more detail in later sections. For the moment, we remark that large sparseness allows extra symmetries and chiral symmetries to emerge. Without the regularity condition, this effect of sparseness is more pronounced because disconnected hypergraphs can be present; the regularity condition eliminates the disconnectedness and mitigates the complication of emergent symmetries, however symmetries can still emerge once sparseness is further increased.
As an indication of the effect of the regularity condition, in Fig. 13, we compare results for the distribution of gap ratios ρ(r) and P(s) with and without the regularity condition for the N = 32 sparse SYK model.
For k = 4 no difference is observed even in the tail of the distribution. For k = 2, the degeneracy only appears in some of the realizations of the non-regular case, which results in a large peak at the origin. (see green square in Fig. 13, left). Here the effect of the hypergraph disconnectedness is concretely at display: for N = 32, k = 2 without regularity condition, often enough a realization misses a fermion, say γ 32 . Hence it is really an N = 31 model in disguise, which is incidentally still in the GOE class [74]. This produces a 2-fold degeneracy because the extra symmetry γ 32 anticommutes with the chirality operator γ c = 32 i=1 γ i . However in this case fixing the γ c chirality, which we always do, is enough to eliminate the degeneracy.
What happens much less often, but still with a non-negligible probability, is that a realization can altogether miss two fermions, say γ 31 and γ 32 . In this case we have an N = 30 model in disguise, which is in the GUE class. This is a 4-fold degenerate situation: 2-fold from the extra symmetries γ 31 , γ 32 , and another 2-fold from the fact that the N = 30 model (GUE) has a time-reversal operator that anticommutes with the N = 30 chirality operator [16]. Fixing the γ c chirality only eliminates the former 2-fold degeneracy and this explains the degenerate data point in the left figure of Fig. 13. Therefore, to reach any firm conclusion from the study of spectral statistics, we have to classify the realizations of the Hamiltonian with all emergent symmetries taken into account (see next section). Therefore, for the study of the critical k for quantum chaos to occur, it is advantageous to rely on regular hypergraphs, which will also exhibit degeneracies and In agreement with the theoretical prediction, we find good agreement with RMT for sufficiently large k. As k decreases, we observe a bump for small spacings which suggests that the spectrum starts to develop a twofold degeneracy. For k = 1, not shown, the degeneracy is exact for some realizations. We shall see that this is due to additional global symmetries induced by the increased sparseness. The regularity condition is imposed. emergent symmetries but only for smaller k's with respect to those in the non-regular case.
In agreement with the theoretical expectation, see Fig. 14, deviations from RMT become more evident as k decreases. The tail becomes gradually exponential and more importantly, for k = 1.25, we again observe a peak in P(s) for very small s instead of the expected level repulsion P(s) → 0 as s → 0. By direct inspection of the spectrum, we have found that, even after the regularity condition is imposed, the peak is related to an emergent spectral degeneracy. As k → 1, an almost exact two-fold eigenvalues degeneracy occurs for some disorder realizations. The peak becomes again very large which prevents a meaningful spectral analysis without further processing of the spectra.
We postpone this analysis to later sections. For the moment, we just mention this degeneracy in the k → 1 limit is related to the existence of additional global symmetries, represented by commuting and anticommuting operators, induced by the sparseness of the Hamiltonian. Once they are taken into account, the level spacing distribution still shows level repulsion but deviates markedly from the RMT prediction. The asymptotic decay is indeed exponential as for Poisson statistics. However, strictly speaking, it is unclear whether the nature of the quantum chaos transition is quantitatively similar to that of the Anderson metalinsulator transition or an chaos-integrable transition. The route to integrability is highly non-universal. In many cases it is not properly a transition but rather a crossover at least from the point of view of spectral statistics. To be specific, harmonic oscillators are integrable and a rectangular billiard is also integrable but the spectral correlations are very different so the transition from chaos to integrability will depend on the integrable system. By contrast, an Anderson insulator has Poisson statistics and the transition can be typically characterized by critical exponents and the scale invariance of level statistics at the transition so it is largely universal. We will return to this point when we investigate the k = 1 case in more detail.
One disadvantage of P(s) is that it requires unfolding of the spectrum. This does not pose any problem for large k, but as spectral degeneracies start to appear for smaller k, it is more challenging to carry out the unfolding procedure. In order to further characterize the deviations from RMT, we investigate the average adjacent gap ratio r which does not require any unfolding and also provides information on the nature of very short-range spectral correlations. Taking the 60% of the eigenvalues around E = 0, we have found that even for k = 1. Although these results are not inconsistent with those from the level spacing distribution, it appears that deviations from RMT predictions are smaller for this correlator. A possible reason for this quantitative difference is that the gap ratio provides information of spectral correlations of shorter range than P(s). In order to confirm this prediction, we compute the full distribution of the adjacent gap ratio ρ(r) using the 90% of the eigenvalues. Results depicted in Fig. 15, are consistent with those of r . Agreement with the RMT prediction is excellent except for k = 1.25. The main difference being the large enhancement of ρ(r) for very small r at k = 1.25. By direct inspection of the spectrum, we associate this peak to an emergent degeneracy of the spectrum. Therefore, even considering only regular hypergraphs, it is not enough to remove these spectral degeneracies related to new global symmetries of the system. The regularity condition only shifts its appearance to even smaller values of k ≈ 1.25.
In summary, we have found that a sparse SYK model with N Majoranas is still quantum chaotic for sufficiently high energies, or temperatures, provided that the probability p ∝ 1/N α with α < 3. For α = 3, spectral correlations are still well described by random matrix theory for k 1. For k ∼ 1, we gradually notice deviations from this prediction. In the k → 1 limit spectral degeneracies are frequently observed which makes the spectral analysis difficult even if the regularity condition is taken into account.
However, the spectral region related to the possible existence of a gravity dual is the edge corresponding to the lowest eigenvalues, and not the bulk of the spectrum. We now move to the study of this region.

B. The edge
In this section we study the spectral correlations of the lowest eigenvalues relevant for the time evolution of the system in the low temperature limit. This is the only region that may be related to a gravity dual.
Technically, it is more challenging to reach firm conclusions because the small spectral window close to the ground state limits substantially the use of spectral averaging to diminish statistical fluctuations. Moreover, the rigorous mathematical results for sparse random graphs are less sharp for the edge of the spectrum as compared to the bulk. As was mentioned earlier, for a random graph, RMT spectral correlations at the edge of the spectrum and a semi-circular spectral density law occur for p 1/L 2/3 where L is the matrix size.
A naive translation of these result to the sparse SYK would lead to a critical scaling p 1/N 8/3 . However, we stress that the results for graphs are not necessarily applicable here and that, even if they are applicable, the bound on p to observe RMT correlation does not have to be optimal, namely, it may be that an even stronger sparsing, such as 1/N 3 , may still lead to RMT correlations at the edge of the spectrum. In order to proceed with the spectral analysis, we obtain only the 2N lowest eigenvalues by exact diagonalization using special techniques for sparse matrices based on the Lanczos algorithm which allows us to reach N = 42. For a given set of parameters, we carry out ensemble average until we have at least 10 4 eigenvalues.
We first investigate the critical scaling p N 4 = kN identified in the bulk of the spectrum. We study the dependence of level statistics on k with the goal to clarify whether spectral correlations are still quantum chaotic and, if so, to identify the approximate critical k = k c . Results, depicted in Fig. 16, show a gradual weakening of quantum chaotic features as k is reduced. An exception to this trend is P(s) for k = 4 and N = 34 which is closer to Poisson than that of k = 2. Presently, we do not have a clear understanding of this anomalous deviation. Results for the adjacent gap ratio, see Fig. 18, indicates that the spectrum, at least for very short range correlations, is quantum chaotic in the large k limit.
In the k ≈ 1 region no level repulsion is observed which indicates that the Hamiltonian is too sparse to sustain quantum chaotic features. It is important to note that, also in the tail of the spectrum, we observe degeneracies of the spectrum for k → 1 though not in all disorder realizations. For the analysis of the spectral correlations, we have removed them ad hoc. This will be justified in the next section by the existence of global symmetries that cause the spectral degeneracies.
A feature of criticality [75,76] is the weak or no dependence of spectral correlations on the system size N. Results depicted in Fig. 17 show a weak N dependence in the k ∼ 1 region. Furthermore, spectral correlations deviate strongly from the RMT prediction. This is consistent with the idea that k ≈ 1 is the maximum sparseness consistent with quantum chaotic features.
The calculation of the adjacent gap ratio r , see Fig. 18, confirms that the maximum degree of sparseness consistent with quantum chaotic features k = k c 1. For smaller k, the gap ratio deviates strongly from the random matrix prediction and approaches the Poisson limit. Our main aim here is identify the region of parameters for which quantum chaos occurs rather than the description and nature of the transition.
Although a transition to Poisson statistics and a critical region with an approximately size invariant spectral correlations are typical of metal-insulator transitions induced by disorder, further investigations , beyond the scope of the paper, would be necessary to reach a firm conclusion.

VI. EMERGENT SYMMETRIES AND QUANTUM CHAOS
Having identified the critical sparseness p N 4 = kN with k 1 to observe quantum chaotic features, we now focus on the limiting case k = 1. Depending on the value of N mod 8, the SYK model for even q is in one of the three Wigner-Dyson universality classes, while the SYK model for odd q is in one of the three chiral random matrix classes. In this section, we show that, for small k = 1, at least six of the ten RMT universality classes emerge from a SYK model for a single value of N in the GUE class. Since the joint spectral density of the superconducting ensembles [77] is of the same general form as the chiral ensembles, our observables cannot distinguish the two. The study of emergent symmetries requires a large ensemble and, although we show some results for N = 34, our main analysis focuses on N = 26 where we can easily generate a large number of disorder realizations with and without the regularity condition.
Even after the regularity condition is imposed, for some disorder realizations, we observe what appears   to be an exact two-fold degeneracy while for other realizations there is no degeneracy or only a quasidegeneracy. In addition, the spectrum for some disorder realizations has a chiral symmetry E → −E while for others there are both chiral and two-fold degeneracy. As an example, in Fig. 19, we depict results for different disorder realizations for N = 34 and k = 1 where only a spectral average is carried out in the central spectral window comprising 90% of the total number of eigenvalues ∼ 65000. Surprisingly, despite the fact that the symmetry for N = 34 is GUE, we observe for some disorder realizations GOE and GSE symmetry. For others disorder realizations, a spacing distribution resembling that of the superposition of two random matrix ensembles is observed.
For a more systematic study, we turn to N = 26 and k = 1 where more disorder realizations can be generated. We start our analysis with an ensemble of 5000 configurations without imposing the regularity condition. To determine if a spectrum has chiral symmetry, we compute[78] In order to investigate the scale to which quantum chaotic features persists, we turn to the connected spectral form factor [20,[79][80][81][82] for the unfolded spectrum, with Z = i e iλ i t−βλ i with λ i the unfolded eigenvalues and β the inverse temperature (only the β = 0 case will be considered). We have removed the disconnected part related to the one-point function. In order to reduce finite size effects, the sum of λ k is cut-off by a Gaussian factor e − λ 2 i 2W 2 (60) with a width W determined such that a significant fraction of the eigenvalues is included in the calculation.
For example, in the case of N = 26 with 4096 eigenvalues, we choose W = 500 or W = 1000. In agreement with previous spectral analysis [20,22], we have observed, see Fig. 21, for k 1, an excellent agreement with RMT predictions even for relatively short times. The smearing of the peak at t = 2π (the Heisenberg time) is a well documented finite size effect.
In Fig. 22  Hermite method [22]. For k = 4, we find a large ramp and saturation in excellent agreement with the random matrix prediction. The peak for short time is a well known finite size effect [22].
have been rescaled to have an average spacing equal to 1. Despite the limitations on the ensemble average to reduce statistical fluctuations due to the different universality classes, we observe very good agreement with universal random matrix results. Deviations from RMT occur at t < 0.5 where we observe a large peak which should not be confused with the peak due to the disconnected part of the form factor. The width of the peak is of the order 1/W, but its area, which is responsible for deviations of the number variance from the RMT results, does not depend on W.
Next we consider realizations with a adjacent ratio of about 0.42 which show a pronounced peak in  [68]. It is simply given by with α, not to be confused with the scaling of probability introduced earlier, the fraction of the realizations in class 1, and the rest, 1 − α, in class 2. In Fig. 23 Since the spectral form factor agrees well with the universal random matrix results, we expect that also the nearest neighbor spacing distribution is given by RMT. In Fig. 24, we show the spacing distribution corresponding to the eigenvalues of Fig. 22. We have excluded realizations with spacings > 5 which actually occur quite frequently. We even have observed spacings of order 100 times the average spacing.
Including these realizations would shift the peak somewhat to the left, but otherwise the agreement with RMT is again good.
We now study the same parameters k = 1, N = 26 but imposing the regularity condition. In this case, about half of the realization have chiral symmetry, and about a quarter are doubly degenerate, but we did not observe higher degeneracies in our ensemble of 5000. For smaller values of k, below k = 1, the number of emergent symmetries rapidly increases. In an ensemble of 5000 disorder realizations, for

A. Origin of the Emergent Symmetries
We now investigate the origin of the emergent symmetries in the sparse limit. For the model without regularity condition, we can imagine an extreme sparseness k ∼ 1/N. Then a realization of the Hamiltonian typically involves only one product of four Dirac matrices, namely with no Einstein summation convention. Such a Hamiltonian has two energy levels with energies ±J if { j 1 , j 2 , . . . , j l } have even (odd) number of common elements with the set of subscripts of every term in the Hamiltonian. A simple example is the following: for N = 10, q = 4 and k = 0.5 with regularity condition, we can for example obtain a Hamiltonian of the form H = J 1357 γ 1 γ 3 γ 5 γ 7 + J 25610 γ 2 γ 5 γ 6 γ 10 + J 34610 γ 3 γ 4 γ 6 γ 10 + J 2789 γ 2 γ 7 γ 8 γ 9 + J 1489 γ 1 γ 4 γ 8 γ 9 .
The symmetries are, and all the operators generated by the above six operators. There is no chiral symmetry for this Hamiltonian.
Note that A 1 , A 2 , A 3 commute with each other and Given the above discussion, it becomes interesting to ask statistically how the number of symmetries and chiral symmetries scales with respect to N at different values of k. A precise study of this question is beyond the scope of the current paper, but we mention on the fly our preliminary numerical observations for spectra obtained imposing the regularity condition in the generation of the Hamiltonian: 1. If k < 1, the number of emergent symmetries grows quickly as N grows.
2. If k = 1, the number of emergent symmetries stays more or less constant (or grows very slowly) as N grows.
3. If k > 1, emergent symmetries only rarely occur and with a frequency that decreases rapidly as k increases.
The presence of emergent symmetries or chiral symmetries can alter the RMT symmetry class naively expected from the corresponding dense SYK model. For example the q = 4 dense SYK model does not have any chiral symmetry and hence always falls into one of the three Wigner-Dyson ensembles, whereas in the very sparse regime of the q = 4 sparse SYK model we see how chiral symmetry can emerge, and hence chiral ensembles can appear. The emergent symmetries can also alter the symmetry class in more subtle ways. We see that N = 26, q = 4 (regular) sparse SYK model, whose dense counterpart always lies in the GUE class, can have realizations in the GOE and GSE classes in the very sparse regime k = 1 (Fig.   24). To explain this we first briefly recap why the dense N = 26, q = 4 SYK model is always GUE. For any even N, there are exactly two independent symmetries for the dense SYK model: a unitary symmetry and an anti-unitary symmetry where K is the complex conjugation and C is the charge conjugation matrix such that Cγ i C −1 = ±γ T i . Since we always look at the eigenvalue statistics in a fixed quantum number sector of the unitary symmetry, for Then we can define a new anti-unitary symmetry which commutes with γ c , and T 2 = 1 or −1 depending on which Dirac matrices A contains. The former case gives us GOE and the latter gives us GSE. In Appendix B we give two concrete examples of this phenomenon. There also can be scenarios where the emergent symmetries give rise to degeneracies but do not change the symmetry class, such as the example shown in Fig. 26. Once these degeneracies are taken into account, so that we fix x i jkl and carry out disorder average over J i jkl only, we observe the following (see Fig. 26): P(s) in the bulk of the spectrum is well described by RMT but only for s < 1, for larger s, the agreement with Poisson statistics is excellent. By contrast, P(s) in the tail of the spectrum comprising the lowest 2N eigenvalues, shows excellent agreement overall with Poisson statistics. Results for the distribution of the gap ratio ρ(r) are qualitatively similar, the bulk of the spectrum is well described by RMT while the tail by Poisson statistics. There is no discrepancy with the level spacing results because the gap ratio provides spectral information of the shortest-range scale, a region where P(s) still agrees with GUE. The tail of the spectrum is close to the prediction for Poisson statistics though we observe a peak at small r likely related to some other emergent symmetry that we have failed to identify.
In summary, once the symmetries are factored out, it seems that even for k = 1, it remains some degree of level repulsion in the bulk of th spectrum that may indicate some remaining quantum chaotic features though deviations from the RMT prediction are very strong. By contrast, in the tail of the spectrum, the results are consistent with Poisson statistics. The latter suggests that the system may have a mobility edge at finite energy. It would be interesting to further characterize the exact nature of the transition though our main motivation here is only to determine the maximum sparseness for which quantum chaos is observed.

VII. CONCLUSIONS AND OUTLOOK
We have investigated the spectral density and spectral correlations of a sparse SYK model as a function of the degree of sparseness. We have identified the maximum sparseness strength consistent with a Schwarzian spectral density, once collective fluctuations are factored out, and quantum chaotic level statistics. These are features of a field theory with a quantum gravity dual. We have carried out explicit analytical calculations of the spectral density moments that have revealed a striking relation between the leading cor- We have fixed the non-zero x i jkl so that the system has a global symmetry that leads to a double degeneracy of the spectrum for all disorder realizations. This degeneracy is removed in the calculation of P(s) and ρ(r). rection due to the sparsity of the SYK Hamiltonian, ∼ 1/(kN), and the leading large d correction of the Parisi's model, a U(1) gauge theory on a d-dimensional hypercubic lattice. As the critical sparseness for quantum chaos is approached, we have noticed the emergence of novel global symmetries that not only induce spectral degeneracies but result in an ensemble that, for a single value of N, contains disorder realizations with level statistics well described by any of the three Wigner-Dyson symmetry classes, and the three chiral random matrix ensembles.
Our results raise some interesting questions: effectively, the sparse SYK Hamiltonian is represented by a sparse random matrix. Can the matrix defined in this way be relevant for matrix models describing quantum JT gravity? Is the critical sparseness to observe quantum chaos of relevance in the description of realistic interacting quantum dots [37]. Is there some explicit relation between Fock-space geometry and space-time so that these sparse SYK models have a natural gravity dual? Is it possible to characterize more generally the connectivity and regularity of a hypergraph so that we can establish the condition for quantum chaos and the existence of a gravity dual in terms of these parameters? About this last point, it would be interesting to study how the sparsity of the random hypergraph affects the early time diagnostics of quantum chaos: the OTOCs and the related diagnostics of operator growth. In particular, it would be interesting to clarify whether the high degree of sparsity has sharp effects on the growth of local operators built out of products of Majorana fermions [83]. It would also be interesting to push further the relation between the sparse SYK model and the Parisi's model to, among other things, to identify the role of the latter in the context of holography. We expect to address some of these problems and questions in the near future.
This emergent symmetry A makes the Hamiltonian belong to GOE through the mechanism described in section VI A.

GSE
We choose the Dirac matrices with the following subscripts to appear in the Hamiltonian This emergent symmetry A makes the Hamiltonian belong to GSE through the mechanism described in section VI A.