Out-of-time-order correlators and Lyapunov exponents in sparse SYK

We use a combination of analytical and numerical methods to study out-of-time order correlators (OTOCs) in the sparse Sachdev-Ye-Kitaev (SYK) model. We find that at a given order of N, the standard result for the q-local, all-to-all SYK, obtained through the sum over ladder diagrams, is corrected by a series in the sparsity parameter, k. We present an algorithm to sum the diagrams at any given order of 1/(kq)n. We also study OTOCs numerically as a function of the sparsity parameter and determine the Lyapunov exponent. We find that numerical stability when extracting the Lyapunov exponent requires averaging over a massive number of realizations. This trade-off between the efficiency of the sparse model and consistent behavior at finite N becomes more significant for larger values of N.


Introduction
Understanding the underlying mechanisms of the spread of information in strongly interacting quantum systems is a challenging problem that connects quantum information, condensed matter theory and quantum gravity.Out-of-time order correlators (OTOCs) provide a useful measure of how fast a local perturbation spreads across a quantum many-body system under time evolution [1][2][3][4][5][6][7][8][9][10][11][12][13][14].Originally introduced as a semiclassical approach to superconductivity [15], OTOCs experienced a renaissance with the discovery of a deep connection between quantum chaos and black holes, where black holes have been shown to be fast scramblers [1,16] and a bound on the rate of growth has been conjectured [3].
The Sachdev-Ye-Kitaev (SYK) model [7,17,18], a system of N Majorana fermions with q-body interactions, i.e. q-local, has been extensively explored as a solvable model of quantum chaos and holography.At low temperatures, the system displays an approximate conformal symmetry and, as in generic black holes, it scrambles information at the fastest possible rate.Among its many variants, the sparse-SYK model [19,20] stands out as a simplified version with fewer interactions which allows one to probe finite N corrections via numerical simulations.Remarkably, the sparse-SYK model with sparseness parameter k and only k N interaction terms in the Hamiltonian, is still a fast scrambler [19,20] and it can be shown to approximately reproduce properties of the original SYK model while allowing to uncover other features as in the spectral form factor studied in [21].Additionally, the sparse-SYK model can serve as a replacement for the original SYK model when modeling traversable wormholes [22,23] including, remarkably, its simulation on a quantum processor [23].
In this work, we study out-of-time order correlators in the sparse-SYK model, both analytically and numerically.Even though the sparse Hamiltonian contains fewer terms than the all-to-all SYK, there are more diagrams contributing to the OTOCs at order O( 1N ).We investigate ladder diagrams contributing to the OTOCs and present an algorithm to sum the diagrams at any given order of 1/(kq) n .However, an analytic expression for the sum of the diagrams to all orders of 1/(kq) n is out of reach because the methods used in the all-to-all SYK are not applicable here.Another of our main results is the numerical study of the OTOCs.We investigate them as a function of the sparsity parameter k, with the goal of determining the Lyapunov exponent.We also explore to what extent the symmetry describing the large N behavior of the OTOCs in SYK is present in the sparse model.

Chaos and OTOCs
Scrambling in quantum systems can be quantified via out-of-time order correlators (OTOCs).In a chaotic system, OTOCs display an exponential growth at intermediate times characterized by a Lyapunov exponent λ L , and the fast scrambling property translates into a Lyapunov exponent that saturates the chaos bound λ L ≤ 2π/β [3].The definition of OTOCs can be motivated as the quantum analog of the classical characterization of chaos using Poisson brackets, where the brackets are promoted to commutators (or anticommutators for fermionic operators) where V and W are generic local Hermitian operators.The operators evolve according to the Hamiltonian H of the system via W (t) = e iHt W (0)e −iHt .The brackets denote the thermal expectation value ⟨•⟩ β = Z −1 tr[.e −βH ] at inverse temperature β = 1/T and Z = tr[e −βH ] is the partition function.Intuitively, this correlator measures how much an early perturbation V affects the later measurement of W .After the expansion of the commutator squared, we see the appearance of a four-point correlator that is out-of-time ordered where ρ = e −βH .In the SYK model, we typically consider W and V as two distinct Majorana operators.The above expression is the standard OTOC that is also referred to in the literature as 'unregularized' OTOC.In the study of quantum chaos at finite N , a slightly modification of this correlator have been considered, named 'regularized' OTOC in which we split ρ evenly between the operators.One of the reasons to work with the regularized version is because this version is less sensitive to finite-size effects, and it can be shown that the extracted Lyapunov exponent is the same as in the unregularized version [24].
From the OTOC, one could derive the Lyapunov exponent by fitting an exponential function a + b e λ fit t to the region of exponential growth.However, this method is known to lead to incorrect results in the SYK model for several reasons.The main subtlety is that, due to strong finite-size effects, the OTOC does not display a well-defined exponential growth window at intermediate times and finite size N .An alternative method was proposed [24], where it is assumed that the OTOC admits the following large N expansion based upon analytical arguments [5,11] for t ≲ 1 λ L log N , which means that F (t) obeys the rescaling symmetry This symmetry is expected to hold for any many-body chaotic model governed by ladder diagrams close to the semiclassical limit.

Sparse SYK
The Sachdev-Ye-Kitaev (SYK) model [7,17,18] is a toy model of lower dimensional quantum black holes consisting of a system of Majorana fermions with all-to-all interactions.The sparse-SYK model [19,20] is a variant of the original model, which we will refer to as all-to-all SYK, with the advantage of allowing for more efficient computer simulation while preserving important black hole physics behavior.The Hamiltonian of the sparse-SYK model with N Majorana fermions χ j , j = 1, . . ., N , and q-body interactions is defined as H = i q/2 j 1 <...<jq J j 1 ...jq x j 1 ...jq χ j 1 . . .χ jq . (2.6) The parameters x j 1 ...jq are either 0 or 1 and they can be defined in different ways leading to different sparse models.The couplings J j 1 ...jq are drawn from a Gaussian distribution with zero mean and variance given by1 where p is the fraction of terms in the Hamiltonian, i.e., the number of the terms such that x j 1 ...jq = 1 divided by N q , the number of terms in the all-to-all version.The parameter J has dimensions of energy and sets the energy scale of the theory.It will be convenient to define the parameter which is such that the Hamiltonian is a sum of kN independent terms.Due to the random nature of the couplings, physical observables are obtained after performing an average over different disorder realizations of the system.If we choose all x j 1 ...jq to be 1, we recover the original all-to-all SYK Hamiltonian H all-to-all = i q/2 1≤j 1 <...<jq≤N J j 1 ...jq χ j 1 . . .χ jq . (2.9) In this case we have p = 1 and we also recover the variance of the all-to-all SYK model One possible implementation of sparseness in the SYK model consists of taking x ijkl = 1 with probability p and x ijkl = 0 with probability 1 − p.The system becomes more sparse as p → 0 and we recover the all-to-all SYK when p = 1.This procedure, known as random pruning, is perhaps the simplest to implement.However, it makes necessary to take an additional average since x ijkl are treated as random variables and can potentially lead to disconnected clusters of Majorana fermions that do not interact with the rest of the system [20].
A different approach to characterize the sparseness is to use regular hypergraphs, in which the Majorana fermions are identified as vertices of a hypergraph whose hyperedges contain q vertices.The regularity condition imposes that each fermion appears exactly the same number of times in the list of hyperedges.The hypergraph is also uniform, meaning that all hyperedges contain the same number, q, of vertices.An important property of randomuniform-regular hypergraphs is that they are expected to be expanders [26].That is, they are expected to have good connectivity even though they are sparse.Note that the q-uniformity condition is related to the Hamiltonian being q-local but the regularity condition is imposed for convenience.A simple counting argument shows that the hypergraph should be kq-regular.Measures of connectivity of this type of hypergraph in the sparse-SYK context were studied in [22].An example of a regular hypergraph is given in Figure 1.

OTOCs in Sparse SYK
In [7,11], it was shown that, in the all-to-all SYK, to obtain the leading order in N and late time contributions to the OTOC, we should calculate ladder diagrams of dressed propagators.This is also true for the sparse model.In this section, we will present an analytic calculation of ladder diagrams in the sparse-SYK.We find that the standard contribution obtained through the sum over ladder diagrams is corrected, in the sparse-SYK, by higher order terms of order 1/(kq) n .In the all-to-all SYK the ladder diagrams can be summed exactly, giving an expression for the OTOC from which a Lyapunov exponent can be computed [7].The same ressumation does not go through in the case of the sparse model and an analytic expression for the ladder diagrams to all orders of 1/(kq) n cannot be obtained.However, we found an algorithm to sum the diagrams up to any given order of 1/(kq) n that we present in this section.The intricacies of the procedure are in Appendix A.

Ladder Diagrams
As a first step, we will analyze the combinatorics related to ladder diagrams in a sparse-SYK model with random pruning [19].In this model, each interaction term in the Hamiltonian comes with a factor, x abcd , which is either 0 or 1 with probability p.This is telling us that each term in the Hamiltonian is deleted with probability 1 − p.In this model, we normalize the couplings to have a p dependent variance, where after averaging, the multiplication of two couplings is only non-zero if their indices are equal (up to permutations of one another).As in the standard SYK model, the ladder diagrams will be the diagrams that contribute to the OTOC [5,7,11], at leading order in late time t.We will now proceed by showing how the ladder diagrams can be summed at zero temperature -although the OTOC is clearly a finite temperature quantity, the zero temperature calculation will be sufficient for showing the general procedure of ladder diagram summation.The zero temperature calculation can be extended to finite temperatures in a fairly straightforward manner by writing retarded propagators along the real time portion of a Schwinger-Keldish contour and Wightman functions between the two real parts of the contour [5,7,12].We find, as in [19], a series in both 1/N and 1/(kq).This means that at each order in 1/N we will find a series expression for a prefactor that is a proportionality constant for each diagram in 1/(kq).This series terminates at an order corresponding to the order of the ladder diagram being considered.At zeroth order in 1/(kq), we have the familiar diagrams of the dense model.But at higher orders, there are corrections organized in a 1/(kq) expansion.To leading order in N , we find that a ladder diagram of O(J 2n ) (think a ladder with n rungs) will contain a series in 1/(kq) up to order 1/(kq) n−1 , i.e.
The dashed line implies that the diagram is disorder averaged.For the rest of this section, dashed lines will be dropped from diagrams, and should be assumed.
Let us now write down the first few terms of the ladder diagram sum The first two diagrams are given by where G(τ ) is the Euclidean-time propagator, i.e. the time-ordered product Let us now focus on the combinatoric factor that appears in front of the integral, and is given generally by (3.3).First, let us look at the overall factor outside of the parenthesis in (3.3).The (−1) n comes from the addition of a fermion loop at each successive order.The (q − 1) n is a combinatoric factor coming from the number of ways Wick contractions can be done on the interaction vertices.Inside the parenthesis, we can find each term by considering the ways ⟨J ia 1 b 1 c 1 J ia 1 b 1 c 2 • • •⟩ can be contracted at leading order in N .Focusing on F 1 , we have ⟨J iabj J iabj ⟩ which gives us back (3.1) a priori, except now there is a sum over i, j, a, b which will give a factor of N q .We will also obtain a factor of p since, x iabj x iabj = x iabj and ⟨x iabj ⟩ = p.Ignoring the Wick contraction combinatorics, which amounts to the (q − 1) factor above, we have This is the 1 inside the parenthesis in (3.3).Every diagram in the sum will contribute the same 1, however the terms of order 1/(kq) and up will change based on the order of perturbation theory we are working in.
The two rung diagram will give There are now two possibilities to consider.The first is when ja ′ b ′ is not a permutation of iab.In this case, the second two terms in the cumulant expansion vanish while the first will give us a factor of 1 as above.When ja ′ b ′ is a permutation of iab, then the combinatoric pre-factor will be where the 3 is from the three different terms in the sum.The last equality can be seen from the identity p = kq!/N q−1 .So, we see that the one-rung diagram has two terms that contribute at leading order in N .One that is k independent and one that is of order 1/(kq).We can proceed in this way for higher order rung diagrams order by order and find a series in 1/(kq) at leading order in N .At each order, care must be taken to consider which indices are permutations of the others and in how many ways they can be made permutations of each other.The procedure is quite tedious; we have worked out the first six pre-factors and present them below for the reader.A general expression for the first two terms and the last term in the series is shown below (3.17)In Appendix A, we present a straightforward and general algorithm for calculating a given coefficient to arbitrary order.
We would now like to make an attempt at summing the ladder diagrams.Since we do not have an analytic expression for the ladder diagram to all orders in J 2 and 1/(kq), we choose to work at order 1/(kq) 2 .We can move between successive ladders (adding rungs to the zero-order diagram) through multiplication of a kernel K n , where K(τ 1 , τ 2 ; τ 3 , τ 4 ) = −(q − 1)J 2 G(τ 13 ) G(τ 24 ) G q−2 (τ 34 ). (3.20) We can now attempt to perform the sum over all ladders by viewing the integral operation in (3.18) as matrix multiplication.We obtain In all-to-all SYK, the sum of all ladder diagrams at early times can be shown to scale as F ∼ e λ L t /N .For large kq the sum of ladder diagrams is the same as the all-to-all case, as expected.However, it is unclear how the corrections of the form 1/kq will contribute to the exponential behavior of F.

Diagram Counting Intuition and Hypergraphs
The importance of ladder diagrams in all-to-all SYK is because they contain all O( 1 N ) contributions to the four-point function.As we showed in the previous section, in the sparse model there are additional contributions (compared to all-to-all) of the form 1 kq .This implies that sparse-SYK has more diagrams at O( 1N ).This may seem counterintuitive at first since sparse-SYK decreases the number of couplings, however, an analysis of how the couplings themselves change elucidates this problem.Note that we have redefined the variance in the pruning model to be where we have used k = p/N N q ≈ pN q−1 /q!.Therefore, the contributions of O( 1 N ) no longer depend on the number of couplings J ijkl , rather the O( 1 N ) counting moves to the pruning coupling x ijkl where n is any positive integer.Intuitively, the right equality follows because if x ijkl = 1, then the couplings with labels {i, j, k, l} exist in the Hamiltonian.Hence, a diagram vertex representing this coupling is non-zero, so we can insert these vertices an arbitrary number of times and the diagram will still exist in the theory.In other words, if couplings with indices {i, j, k, l} exist with probability p, then the probability of the same couplings contributing twice is not less likely (such as p 2 ), because said couplings already exist in the theory.Another version of this statement is that since x ijkl can only be 1 or 0, 1 n = 1 and 0 n = 0. Additionally, any non-zero term will have the same 1/kq contribution from the J's, specifically the nth order ladder diagram will have 1/(kq) n from the J's.Therefore, the order of any specific term in O(1/kq) is determined by how many p's (from the pruning couplings) contribute in that term.Hence, the 1/kq counting for any given term is also determined by the sparsity couplings x ijkl .In summary, this means that the O( 1 N ) counting moves to the x ijkl couplings, and since these behave differently from the J ijkl couplings, the contributions of O( 1N ) should change in the sparse pruning model.
An open question remains: we have performed this analysis using the random pruning method, yet the preferable numerical choice appears to be the hypergraph method, do we expect the ladder diagrams in the hypergraph method to yield the same expansion in 1/kq?While we cannot comment on the expected behavior for a fixed hypergraph configuration as any given diagram would be dependent on which couplings are present, we also find in the subsequent section that a fixed hypergraph is insufficient to extract the Lyapunov exponent using the known numerical methods for the SYK model.It is only upon averaging over hypergraphs that the Lyapunov exponent can be extracted accurately, it is important to note that averaging over random pruning and averaging over hypergraphs (for sufficiently large averages) will lead to the same behavior.This is because when averaging over hypergraphs, the hypergraph configuration is chosen randomly, in which case (upon sufficient averaging)  1: Time comparisons between all-to-all SYK and various sparseness (k's) in sparse-SYK using the hypergraph method over a range of N 's from N = 26 to N = 36.All simulations were run on a single processor with β = 1 and q = 4, averaging over 112 samples to account for time differences in individual runs.Times are rounded to the nearest second.
each coupling will occur will probability p by construction.This is analogous to the random pruning model and hence we expect the same behavior, albeit the hypergraph should require less averaging due to its higher connectivity.

Lyapunov numerics
Numerical methods are needed to study the Lyapunov exponent at finite q and finite coupling in the all-to-all SYK because no closed form is currently known in this regime2 .The need for numerical methods to study the growth of OTOCs is even more pressing in the sparse-SYK model, where the general analytic relation between the Lyapunov exponent, λ L , and the sparsity parameter k is unknown.As we showed in the previous section, the early time behavior of the OTOC contains 1/kq contributions and it is unclear how to account for these contributions in models for extracting analytically the Lyapunov exponent.
In this section we numerically simulate the behavior of the regularized OTOCs (3.2) at finite N in the sparse-SYK model.We employ the regular hypergraph method of generating sparseness, as it ensures greater connectivity for any given realization of the Hamiltonian.The qualitative behavior from varying k is shown in Figure 2, for which we see a decrease in the rate of exponential growth as k decreases.This behavior is physically expected, as the sparse model reduces the number of connections between fermions, and hence it decreases the sensitivity to operator perturbations.At this stage, the computational benefit from sparse-SYK is undeniable.Table 1 demonstrates that even for relatively low N , sparse simulations reap massive efficiency gains, an effect that further widens the higher N one chooses to probe.
If one wishes to capitalize on these efficiency benefits to quantify quantum chaos by means of a Lyapunov exponent, one needs a model for use in extraction.In the all-to-all of the all-to-all model (denoted 'a2a') and various sparsity parameters k for the sparse model.All curves use the same number of fermions (N = 30), number of fermions per interactions (q = 4), temperature (β = 1), and number of disorder realizations (D = 1000).For each k, hypergraphs are averaged over, a choice justified in the main body of the paper.Higher k's closely resemble the curves of the all-to-all case, however as k → O(1), i.e. the limit still expected to retain chaos, the exponential growth starkly decreases.model, the best known method for accurate extraction of the Lyapunov exponent [24] utilizes a symmetry of the OTOC expansion which is invariant under We know the sparse model will approach the all-to-all model in the large k limit, and hence this symmetry is expected to be recovered in said limit.Our aim is to determine whether the aforementioned symmetry still holds in the sparse-SYK model, even if only approximately, in such a way that it allows a quantifiable extraction of the Lyapunov exponent for a given k.

Extraction Methods
Typical numerical methods for extracting the Lyapunov exponent via the OTOC involve performing an exponential fit in the exponential regime of the OTOC [27][28][29][30].However, this method was found to be insufficient at finite N and q = 4 for the all-to-all model, where the OTOC at finite-size could only be matched to a simple exponential for a very short time interval.It has been shown [24] that a numerical derivative method exploiting the symmetry (4.2) is most effective in extracting the Lyapunov exponent.For each value of N , a fitted Lyapunov exponent λ fit (N ) was found as follows.Since a shift under the symmetry (4.2) leaves the OTOC invariant, we can pick a constant value of the OTOC, say F * , between all N 's simulated and use interpolation to find the time t * associated to each F * .It follows that the fitted Lyapunov exponent between two curves (assuming the symmetry holds) is which could then be plotted vs 1/N to see the effect as the system approached large N .Taking a second order fit in 1/N for λ fit , λ L is defined as the value when 1/N → 0 In our simulations, we use F * = 0.4.See Figure 3 for a concrete example of this procedure.We test the robustness of this choice by using different values in Appendix C. Due to the success of the numerical derivative method exploiting the symmetry (4.2), we choose to closely match the numerical tools employed by [24] in their extraction of the Lyapunov exponent.This includes utilizing dynamite [31], a Python package allowing fast simulation of quantum dynamics with Krylov subspace implementation [32].As stated in Section 2.1, a regularized OTOC has been shown to reduce finite-size effects, and here V and W are chosen as distinct single Majorana operators.Additionally, it is known for large N that the thermal expectation value involving the trace can be closely approximated by calculating the expectation value via Haar random states [24,33,34], i.e. for a Haar random state |ψ⟩ The finite N error in this approximation can be reduced through sufficient averaging over different initial random states, represented by the overline.It has been shown that averaging these states separately or simultaneously with the couplings J ijkl in the all-to-all case had no significant difference in error.Hence, in this work, we also employ simultaneous averaging of the couplings and initial random states.We find another source of averaging that is important to the extraction of the Lyapunov exponent not present in the all-to-all model, namely averaging over many hypergraphs (by  generating regular hypergraphs using different seeds).As we will explain below, the fixed hypergraph method does not produce the expected symmetry.Only after averaging over many hypergraphs does the desired symmetry becomes apparent (see Appendix C).Thus, only for an aggregate of hypergraphs the symmetry-based method could be used to extract the Lyapunov exponent.
In order to quantify the error of λ L , we simulate a large number of realizations ≳ 100, 000 and subdivide these realizations into subgroups, from which we extract the Lyapunov exponent for each subgroup.By performing statistics on many of these subgroups, we can derive a mean Lyapunov exponent λ, a standard deviation σ, and a coefficient of determination R 2 .We explain this procedure in more detail in what follows.

Averaging over Hypergraphs
The ultimate goal of this section is to glean whether the symmetry as described by (4.2) holds in the sparse-SYK model, and if so to extract Lyapunov exponents and determine their dependence on the sparsity parameter k.We can get a qualitative idea of whether the symmetry holds from Figure 4, where we plotted the Double Commutators of the sparse model in comparison to the all-to-all model over D = 1000 realizations.These results employ simultaneous averaging over hypergraphs, couplings, and random initial states.At this level of averaging, the inconsistent behavior in spacing between subsequent curves suggest that the symmetry breaks down as sparsity increases.This behavior admits at least three possible explanations.Either the symmetry itself does not hold in the sparse model, the sparse model is more susceptible to finite-sized effects, or the averaging over hypergraphs produces random fluctuations requiring a higher degree of averaging.
A quantitative approach must be employed to separate these possibilities.To test the impact of random fluctuations, we note that an average over a sufficiently large number of realizations must converge to a certain value, and therefore repeating the experiment using the same number of realizations should not change that value.We measure the dependence on the number of realizations by repeating the extraction of the Lyapunov exponent many times, where each extraction utilizes D g realizations.From these repetitions, we obtain a mean value for the Lyapunov exponent λ and a standard deviation σ for each D g .By varying D g , we can analyze how the error on λ is modified.The choice of how many repetitions to average over is completely arbitrary, we find that the choice yielding the smoothest results for Lyapunov extraction is to average the "maximum" number of times, i.e. if there are D tot total number of realizations, then we can make D tot /D g unique groupings of D g realizations and hence we average over all of these groups.When D tot /D g < 10 we stop increasing D g , so as to allow sufficient samples to be averaged over.Another possible choice is to choose a set number of samples to average over regardless of D g , we do this with 10 and 20 samples in Appendix C as a test of robustness.
The other quantitative tool utilized is the coefficient of determination R 2 , which represents in our case how well the numerical derivative data of a given group fits (4.4).It has the formula The term SS res refers to the sum of squares of residuals, SS res = i (x i − xi ) 2 , x i refers to our numerical data and xi refers to the prediction based upon our model from (4.4).This represents the total sum of variances between our model's prediction and the numerical data.SS tot on the other hand is the sum of the variance of the data, The ratio SS res /SS tot , also called the fraction of variance unexplained, represents the ratio between the "unexplained" variance (from the discrepancy between model and data) and the total variance.If a set of data implements the symmetry given by (4.2), λ fit will closely match the data, implying SS res /SS tot → 0 and therefore R 2 ≈ 1.Hence, by analyzing how R 2 changes with D g , we can see whether simply greater averaging is sufficient to recover the symmetry.
Due to the outputs of the mean Lyapunov exponent, associated standard deviation, and R 2 all roughly lying in the range [0, 1], we can plot these values together in the same plot,  which can be seen in Figure 5.These plots indicate that random fluctuations, caused in large part because hypergraphs must be averaged over, have a much larger effect for sparse models.It is not until far larger numbers of disorder realizations in the Sparse case that similar values for error and R 2 as compared to the all-to-all model are achieved.At this level of disorder groupings, the all-to-all model demonstrates convergence in the mean, standard deviation, and R 2 values with groups of roughly D g ∼ 10, 000, with R 2 → 1 suggesting the symmetry based model closely matches the data.As k decreases, the stability of this convergence slightly decreases, with somewhat higher fluctuations.However, given a large enough D g , such as D g ∼ 175, 000 for k = 4, we can also recover R 2 → 1 and some convergence in error and mean Lyapunov exponent.We illustrate this point in Figure 6.
Overall, we find that as k decreases, the number of disorder realizations goes up to (and exceeds) an order of magnitude greater than in the all-to-all case.This suggests the following interpretation: as k decreases, the structure/connections of the sparse Hamiltonian further and further diverge from that of the all-to-all model.To compensate for a given hypergraph being a poor representation of the all-to-model from which it is constructed, hypergraphs (and hence the connections/structure) must be averaged over if the symmetry is to be recovered.The subsequent section provides further evidence for this interpretation.
One important observation from Figure 5 is the appearance of the k = 50 average Lyapunov exponent being higher than that of the all-to-all model, this is surprising as increasing sparseness should decrease the chaotic properties.This is however an artifact of numerical error, specifically that the Lyapunov exponent approaches its asymptotic value from above.
At our level of disorder realizations, the standard deviation and R 2 (for all k) make it clear that full convergence has not yet been achieved, the k = 10 and k = 4 cases demonstrate that as one approaches convergence by increasing D g , the average Lyapunov exponent continues to decrease.They also demonstrate that for lower k, convergence is not achieved until higher levels of R 2 are achieved as compared to the all-to-all model.We expect that given enough disorder realizations, the extracted Lyapunov exponent for k = 50 would drop below that of the all-to-all model.A reasonable cause for why the standard deviation may asymptote above zero is due to finite N effects, and it appears for the k = 50 case this asymptote may be above that of the all-to-all case.It is also possible that some disagreement with symmetry-based models is due to finite N behavior acting more strongly in the sparse model.Hence, our next goal is to investigate whether this is the case, and if simulating larger N can lead to more accurate extraction of the Lyapunov exponent.

Is accurate extraction possible at larger N ?
It is natural to ask if going to larger N will allow us to extract the Lyapunov exponent more accurately, as it does for the all-to-all model.Unfortunately, we find that at fixed k, the answer seems to be no.In this section, we probe higher N 's in the sparse model to test if finite-sized effects can explain some of the disagreement with symmetry-based models.Our methods match those used in the previous section, except rather than comparing different k's we compare different ranges of N 's for the k = 4 model, which approaches the minimum limit of sparsity while still retaining chaos, and is correspondingly the fastest of our choices of k to simulate.We chose to retain the lower bound of N = 16 when extracting the Lyapunov exponent for all higher N 's as it led to better agreement with symmetry methods.In Appendix C we have also included an alternative analysis with a constant range of N 's (e.g. for N = 42 the range would be N = 22 to N = 42).
The corresponding plots for these simulations can be found in Figure 7.We find that not only do random fluctuations actually increase in size with higher N , but also the R 2 value has no better agreement as N increases.There is a likely explanation for these trends.Recall that the number of interactions scales as ∼ N q for the all-to-all model and as ∼ kN in the sparse models.Therefore, as N increases, the ratio of interaction terms between the all-to-all and sparse models grow as N q−1 .Our results in Section 4.2 show that the further one strays from the all-to-all configuration, the more averaging one must do to recover the symmetry.Thus, as N increases and the difference in interaction terms grows larger, the sparse model requires higher levels of averaging to recover the symmetry.
Hence we have demonstrated that for low k, there does not appear to be a simple method to consistently extract the Lyapunov exponent at finite N using currently known methods without averaging over massive number of realizations.Another approach one could employ is varying β, however it is unlikely this would lead to better behavior, certainly for β > 1 as fluctuations rise in this case, in part due to fewer values of N being suitable for extraction since the OTOC expansion being valid requires N ≳ βJ.

Efficiency and Lyapunov Exponents
The previous sections have made evident that higher degrees of sparsity (lower k) require greater averaging (and hence greater computation times) in order to accurately extract the Lyapunov exponent.We now seek to determine whether the computational benefit inherit in sparse-SYK is retained after the large amounts of disorder averaging that must be performed.
In this section, we present a rudimentary method to determine the relative efficiencies of sparse models as compared to all-to-all, granting that more sophisticated methodology could be employed to perhaps gain greater insight as to how they compare.Based upon our methodology, we find that the computational trade-off between inherent speed increase in sparse models and time increase due to greater disorder averaging roughly cancel out, meaning accurate Lyapunov extraction is no more efficient for low k models as it is for all-to-all models.We also present the numerically calculated Lyapunov exponents for our given range of k's.
Our rudimentary choice of calculating efficiency starts with choosing what qualifies as "accurate" Lyapunov extraction.Given our ranges of simulated D g , a goodness of fit of R 2 = 0.9 was our choice of where to begin considering the extraction as accurate enough.Since we do not have exact curves corresponding to how R 2 varies with D g , we instead chose to analyze the corresponding values of D g where the numerically calculated values of R 2 were within 0.005 of R 2 = 0.9.This gives us a range of D g values, interpreted to mean that on average, samples extracted with those values of D g will produce curves for λ fit with a goodness of fit of R 2 ≈ 0.9.For data that passes through R 2 = 0.9 multiple times, the largest continuous cluster is chosen.
We also need to quantify the time advantage sparse-models have over all-to-all models, which can be straightforwardly done by finding the time (on average) it takes to simulate a single disorder realization for our range of N 's (N = 16 to N = 36).The full list of times per N is provided in Appendix B, we present in this section just the total time for all values of N .By dividing the total time (in seconds) by the total time of the all-to-all model, one achieves a "time factor", representing the ratio of how much faster a given sparse model is per realization as compared to the all-to-all model.Multiplying the previously described ranges of D g by these time factors will yield the relative number of disorder realizations needed to achieve accurate (as deemed by R 2 = 0.9) Lyapunov extraction.For the calculated values of the Lyapunov exponent and associated standard deviation, we average over the data with R 2 > 0.9 (rounded within 0.01), which all achieve convergence within those ranges.Results of this analysis can be found in Table 2.For discussion regarding the appearance of the k = 50 extracted Lyapunov exponent being higher than that of the all-to-all model, refer to the end of section 4.2.
This analysis demonstrates that the necessity of greater disorder averaging required for accurate Lyapunov extraction within sparse-SYK roughly negates the computational benefit provided per realization by sparse models, since the relative efficiency is roughly comparable between all models.
Efficiency Table ( The time factor is the fraction of time (compared to all-to-all) spend per disorder realization.The D g range is the range of disorder groupings which yield a value of R 2 = 0.9 within 0.005.The relative efficiency is the range of D g multiplied by the corresponding time factor.The Lyapunov exponent and standard deviation is the average of the corresponding quantities over the range of D g 's with R 2 > 0.9 from Figure 5.

Conclusions
We studied the quantum chaotic behavior in the sparse-SYK model as diagnosed by out-oftime ordered correlators (OTOCs).We approached the problem through both analytical and numerical methods.In the large N limit, where the OTOCs are given by a sum of ladder diagrams, we found an algorithm for calculating the 1/kq corrections to the ladder diagrams.This allows us to obtain analytic expressions for ladder diagrams of arbitrary order expanded to arbitrary orders in 1/kq -see Appendix A for details.Importantly, the kernel in (3.18) obtained in this procedure is dependent on the ladder diagram order, a feature not shared by the all-to-all SYK.Thus, the summation over ladder diagrams needs to be truncated at some order in 1/kq, and the summation to extract the corrections to the early time exponential behavior of the OTOC cannot be performed with the same tools as in the all-to-all model.The numerical approach employed in this work followed the symmetry-based method in [24], previously applied to the all-to-all SYK, which is the only known accurate method for extraction of the Lyapunov exponent from finite N simulations.We found that sparse models with k ≲ 50 require far larger averaging of disorder realizations for accurate extraction as compared to the all-to-all model.In a fixed regular hypergraph sparse model, the extraction of the Lyapunov exponent led to inconsistent results, i.e., two fixed regular hypergraphs with the same level of sparseness resulted in different Lyapunov exponents even when we significantly increase the number of disorder realizations.This indicates that the symmetry (4.2) does not hold for a fixed hypergraph configuration.In order to recover the symmetry, we performed an additional average over regular hypergraphs with the same sparseness parameter k.We find that even with the extra average, the fluctuations in the Lyapunov exponent are still larger compared to the all-to-all SYK, and these fluctuations increase even more as k decreases.
As shown in Section 4, a reliable extraction of the Lyapunov exponent is still possible for low values of k but a very large number of realizations is needed.Therefore, there is a tradeoff between efficiency and variance in the extracted Lyapunov exponent in the sparse-SYK.One could argue that these fluctuations can be reduced through probing larger values of N .However, as we argued in Section 4.3, the opposite effect takes place for Lyapunov extraction in the sparse-SYK model.The reduction in computation time of a single realization is offset by the larger number of simulations needed to be performed.In our numerical tests, the offset of the two effects is roughly comparable, i.e., the computational resources required to achieve similar levels of accuracy for Lyapunov extraction are roughly the same between all-to-all and sparse models.
We remark that our numerical results are consistent with the previous statement in [20] that quantum chaos is present in the sparse-SYK model for k ∼ O(1), derived from level statistics [35].Interestingly, the analysis of [20] also revealed the existence of emergent symmetries in the sparse-SYK model.It is possible that these emergent symmetries are correlated with the large fluctuations observed in the Lyapunov exponent.
Our results suggest several directions to explore in the future.
• In [36] the authors calculate the Lyapunov exponent in a variant of the SYK model by solving the Schwinger-Dyson equations numerically.It would be interesting to investigate if this method can be an alternative way to extract the Lyapunov exponent of the sparse-SYK model.
• Our numerical results indicate that the early time behavior of the OTOCs is affected by the sparseness of the model.It would be interesting to further explore the summation of the ladder diagrams to determine the leading corrections to the early time exponential ansatz F ∼ e λt /N in the sparse-SYK model.
• A better understanding of the emergent symmetries found in [20], especially their effect at larger N , could also provide insights on how the sparseness affect the behavior of OTOC expansion (2.4).
Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper.URL: http://www.tacc.utexas.edu.

Appendices
A Sparse SYK Ladder Diagram General Algorithm

A.1 Simplifying Diagram Calculations
As seen from the calculations in Section 3, general ladder diagram calculations will heavily involve combinatorics, especially as their order increases.Additionally, since the kernel is dependent on the order of the diagram, finding a given ladder diagram to arbitrary order is a non-trivial procedure.However, as we will see from calculating the first several ladder diagrams, a pattern emerges.In this section, we will develop the notation and motivation behind this procedure.
Let us return to the example of the two rung ladder diagram General ladder diagrams require c = c ′ to be O( 1 N ).In this case, we notice the Gaussian moment expansion (3.9) simplifies to Recall from (3.1) that only if two couplings' indices are equal (up to permutation) may the multiplication of the two be non-zero after averaging.Let us adopt a convention of labelling vertices from left to right in the ladder diagram.At the order of 1/N the vertices connected by rungs have the same indices, so we can label them by the same vertex number.Our convention will involve replacing couplings by their corresponding vertex numbers, for instance J iabc → 1.In this new labeling scheme, (A.1) corresponds to Similarly, we can replace (3.9) with In our new convention, the term ⟨12⟩ J should be interpreted to mean that after averaging, ⟨12⟩ J = J 2 kq , only if vertices 1 and 2 have indices which are equal up to permutations of each other.Recall that the two relevant cases here are if vertex 1 and 2 are or are not permutations of each other.In the latter case (not permutations), ⟨12⟩ J = 0 (not permutations), so only the term ⟨11⟩ J ⟨22⟩ J contributes.Let us represent the sparsity couplings ⟨x abcd x abcd ⟩ = ⟨x abcd ⟩ as ⟨11⟩ x = ⟨1⟩ x , we found the pruning couplings contributed a factor of ⟨1⟩ x ⟨2⟩ x = p 2 .
In the other case, vertices 1 and 2 do have indices equal up to permutations.However since all vertices are equal, the averaging over x reduces to ⟨1122⟩ x = ⟨1⟩ x = p, so overall only a single factor of p is contributed from the x's.Recall that the total expression (in our new notation) is ⟨1122⟩ J ⟨1122⟩ J , so the factors from averaging over J's must also be considered.Overall, our goal is to find the numerical/combinatorial terms multiplying these p's since these are what determine the order in 1/kq, which are embedded in averaging over the J's.By using (3.9) and dropping the factors of J 2 (these can be pulled out across the entire expression), let us create a table which counts all relevant contributions: Let us repeat this process for the ladder diagram with three rungs in full detail.For this case we have Let us organize again by order in p, the case p 3 is given by no vertices being equal up to permutations.In that case, as was the case for n = 2, only the first term contributes.For order p 2 , the counting gets trickier.There are three cases which contribute, if vertices 1 and 2 are equal, vertices 1 and 3 equal, and vertices 2 and 3 are equal.For 1 and 2 equal, we have ⟨13⟩ J = 0 = ⟨23⟩ J , leaving a coefficient of 1 + 2 = 3.Each other case of order p 2 has this same numerical coefficient, making the total numerical coefficient 9.For order p, all terms contribute giving a numerical constant of 15.Since all we care about is the order of p, not the specific case of which vertex is a permutation of another, we notice that every term with coefficient 2 in (A.6) is the same up to relabelling of vertices, therefore they contribute at the same order of p.We therefore wish to label these terms by a factor which suggests that they all contribute equally at the same order of p.If the highest order of p a term contributes is p k , we will relabel these terms by the factor P k .In this way may rewrite (A.6) as Our current notation is still burdensome, as the expansion contains terms with both ⟨•⟩ J and ⟨•⟩ x , our goal will be to combine them.We can do this by notationally decomposing our average over x into groups of 2 vertices (like that which is naturally done from the Gaussian moment expansion for the J's).To see this concretely, let ⟨•⟩ d be a term in the decomposed average, such that ⟨aabb . ..⟩ x ≡ ⟨aa⟩ d ⟨bb⟩ d ⟨. ..⟩ d .Note that these decomposed terms are not a true average nor should they be interpretted as one, they can recombine/decompose freely (to conserve the relevant structure) but only after a total recombination of the ⟨•⟩ d 's do they represent something physically meaningful, namely it recovers the original average over x.
Their purpose is simply for notational convenience.The use of this becomes apparent when we let ⟨ab⟩ ≡ ⟨ab⟩ J ⟨ab⟩ d , meaning we have the following identities The second identity is true because the term is only non-zero when ⟨ab⟩ J ̸ = 0, meaning vertices a and b are permutations of one another.In that case, ⟨abab⟩ x = ⟨a⟩ x = p, meaning the highest order the term can contribute is p 1 , hence the term corresponding to P. The third identity can not be further simplified (ie, we cant set it equal to P) because if ⟨bc⟩ ∼ P, that would imply ⟨bc⟩ ⟨bc⟩ ∼ P 2 , which contradicts the second identity.The origin of this comes from equivalent vertices appearing in pairs in the ladder diagram.
Another way to see why we make this replacement is to note as stated earlier that a relabelling of indices will not change the way in which the counting at each order of p is done, therefore we can always choose to relabel terms such that they combine into the form (A.7), for example Note that the process of finding the coefficient at O(p k ) is not as simple as just reading off the coefficient of P k as written above, recall for example the term ⟨11⟩ ⟨22⟩ ⟨33⟩ ∼ P 3 contributed for every possible case of order p 2 .Therefore we need to multiply each coefficient by all the cases where it contributes.Let us define two vertices as "connected" (forming a "connection") if their indices are equal up to permutations.Additionally, we will define the number of "free" vertices as the number of vertices whose indices may take on independent values while still ensuring the term is non-zero.For instance, ⟨aa⟩ ⟨bb⟩ has two free vertices, since vertices a and b can both take on independent indices.On the other hand, ⟨ab⟩ ⟨ab⟩ is only non-zero if a and b are connected, meaning there is only one free vertex (say a, since b must have its vertices equal up to permutation of a).
For the third-order ladder diagram, there must therefore be one additional connection (involving two vertices) for the leading term to contribute at order p 2 .The term ⟨11⟩ ⟨22⟩ ⟨33⟩ has three free vertices (1,2, and 3), of which we choose two to be equal at order p 2 , hence its contribution will be multiplied by 3  2 .This is not the case for the term ⟨11⟩ ⟨23⟩ ⟨23⟩; vertices 2 and 3 already must be equal for the term to be non-zero.
The general counting procedure will therefore be as follows: create a table with the left column denoting the order of p, and the top row corresponding to the expansion of the average of the couplings in terms of P. For order p k , take the coefficient of P k and all higher powers of P and multiply them by the number of combinations required to form enough connections to be of order p k .For instance at order p, ⟨11⟩ ⟨22⟩ ⟨33⟩ must undergo 2 connections involving 3 free indices, so that means out of 3 free vertices choose 3 to be connected, which is 3 3 = 1.Similarly for ⟨11⟩ ⟨23⟩ ⟨23⟩ there are 2 free vertices (2 or 3 must be equal to the other to be non-zero, removing a free vertex) with one connection involving 2 vertices, implying a factor of 2 2 = 1.In general this pattern will hold.Another reason to see why their combinatoric factors should be 1 is that at order p, all vertices are equal so there is only one way to combine them all.So our table for the third order ladder diagram becomes We will continue to 4th order to further demonstrate how these combinatoric factors work.At order p 2 , the factor P 4 undergoes two connections, however this time there are two cases.With four free vertices, three vertices could be all connected (making 2 total connections) or two sets of two connections could be formed.The latter of which has 4  2 for the first choice and For this procedure to work for arbitrary order of ladder diagram, we must find a closed form for the expansion of averaged couplings into P's as done above, and find a general method for calculating the combinatoric factors.These will both be addressed in the next section.

A.2 General Procedure
Let us first address the combinatoric factors.Recall based on our definition of free vertices that the monomial P k has k free vertices.At order p k no connections need be formed, so the term does not gain a combinatoric factor (although still has a coefficient multiplying P k ).The term P k+1 has k + 1 free vertices, and must undergo 1 connection to be of order p k .Out of the k + 1 free vertices, two must be chosen to form a connection, hence the combinatoric factor becomes k+1 2 .For P k+2 , there must be two connections between free vertices at order p k .However we must choose how to partition these connections, they can either all be grouped together, meaning three free vertices are all permutations of one another (2 connections) or be separated into two different connections (1+1 connections).For any grouping of connections, there is always one more vertex than connection, which is what the combinatorics are performed on.So for the two connections grouped together, the combinatoric factor will be k+2 2+1 = k+2 3 .For two separate connections, we do the combinatorics out of all the free vertices with one connection k+2 1+1 and then do combinatorics on the remaining free vertices for the next connection k 1+1 making a total combinatoric factor of k+2 2 k 2 .
(A.12)So at order p 2 , we start with the combinatoric factors of P 6 .6 − 2 = 4 connections must be made.These connections can be partitioned into the following sets: {4}, {3, 1}, {2, 2}, {2, 1, 1}, {1, 1, 1, 1}.We add one to each of these sets to represent the possible combinatorics on the free vertices, so {5}, {4, 2}, {3, 3}, {3, 2, 2}, {2, 2, 2, 2}.We then need to find the combinatoric factors by taking the j = 6 free vertices and doing successive choose functions on each member of the set, this will give a factor of the form Notice that terms like 0 2 = 0 = 1 2 make those coefficients vanish, so those terms do not contribute.These factors effectively mean it is not possible to construct the sufficient connections with that partition set.Repeating this process for all P's along with using the coefficients from (A.12) yields Since ⟨112233445566⟩ contains 6 powers of J 2 abcd ∼ 1/p, at order p 2 to total contribution from p's is 1/p 4 ∼ 1/(kq) 4 , which matches our coefficient in the expansion from (3.11).
The general procedure for combinatoric factors should become clear from this.At order p k and ladder diagram order n, the combinatoric factor for P j is as follows: if k > j, then the factor is zero.If k = j, the factor is 1 (effectively zero choices of free vertices need to be made).If j > k, then j − k ≡ m number of connections need to be made.The number of distinct combinatoric factors is equal to the number of partition numbers of m, since m different connections must be partitioned into groups.For any given group of connections, the number of vertices is one more than the number of connections in said group.Create a list corresponding to each possible partition of m (so the number of different lists will be the number of partitions of m), and add 1 to each member of the list to represent the number of free vertices.For a given list m l , place each member m ls ∈ m l into the lower value of a choose function, • m ls , and multiply all of these choose functions together.The upper value will start with j (coming from the free vertices available from P j ), and then decrease by the number of free vertices connected m l1 , which will be of the form • • • .This process will repeat for all sets of partitions m l , and then for all values j from P j such that j > k.After multiplication by the corresponding coefficient of P j and summing all terms together, this will form the total factor contributing to p k .
The last step of our procedure is to find the expansion of a n (P) in terms of P's for arbitrary n.The first step is to form a recursion relation for ⟨11 • • • nn⟩ in terms of P's.There are 2 ways for vertices 1 and 2 to be paired, and n − 1 different vertices that may be paired with vertex 1.Since we have the ability to relabel vertices without impacting the expansion of P's, we arrive at where the underbrace refers to the number of unique vertices (before considering permutations).Since ⟨22 • • • nn⟩ involves n − 1 unique vertices which we are free to relabel, it is therefore equivalent to a n−1 (P).Additionally ⟨11⟩ = P, which leaves the second term.We can choose to define the second term as b k (P) and find its recursion relation by again expanding and relabelling Using (A.9) and our above expressions we can rewrite our recursion equations as a set of two coupled recursion relations These sets of equations can be solved to find

B Computational Specifications
All numerical simulations presented in this paper were done through the Texas Advanced Computing Center (TACC) on the Frontera supercluster [38].The supercluster is comprised of 8,386 nodes (of which a maximum of 200 were ever used per simulation), with each node contains 56 cores.For ease of access, compute node specifics are transcribed here:

C Fixed Hypergraph Configuration & Tests of Robustness
The analytical tools used in Section 4 can be applied to a fixed hypergraph configuration, results of which can be seen in Figure 8.For this set of figures the mean, standard deviation, and R 2 must be separated as their outputs lie in vastly different ranges.The lack of convergence and small value of R 2 demonstrate that a given fixed hypergraph exhibits complete disagreement with the symmetry based methods.These results remain true for various choices of hypergraph seeds, k, and D g .As a check that the results of this paper are robust to the various choices of analysis used in the main body of this work, we provide here alternative methodology to demonstrate our conclusions remain consistent.The first encountered choice which could make a significant change is which value F * is chosen in the analysis.The main body of the paper utilized F * = 0.4, however in Figures 9, 10 we demonstrate that the behavior exhibited in Figure 5 is consistent with F * = 0.3 and F * = 0.5 respectively.
The next natural choice of variation is the size of the groupings statistically analyzed for each D g , as it may be argued the proper choice would be to have a consistent number of samples regardless of D g .Hence, we provide plots for two different choices, the first being 10 samples (Figure 11) and the second being 20 samples (Figure 12).These plots are compared to Figure 5 from the main body, they exhibit a higher degree of fluctuation for smaller values of D g .This result is expected, as lower values of D g (meaning a given Lyapunov exponent is extracted with low number of realizations) will inherently have higher error than for higher values of D g .This inherent instability makes Lyapunov extraction more inconsistent, hence the choice of extraction used in the main body.However even given these choices, our main findings (namely R 2 → 1 can be achieved at low k with sufficiently high D g ) still hold.
Lastly, we analyze whether including the same lower bound N = 16 regardless of the maximum simulated N is a good choice.Instead, we could keep a set number of N in a given range (in this case, 11 different N 's) leading to a change in the lower bounds.Hence the ranges compared are N = 16 to N = 36, N = 18 to N = 38, N = 20 to N = 40, and N = 22 to N = 42, which can be found in Figure 13.As compared to Figure 7, these simulations actually exhibit worse correspondence with the symmetry based model, and exhibit higher amounts of random fluctuations.This behavior is consistent with our arguments from Section 4.3; by disregarding data points from lower N values, the remaining data (higher N 's) of the sparse models have higher discrepancy with the all-to-all model, hence requiring more averaging.This reduces the maximum D g (and hence the range) by 1 2 .As a result, higher values of R 2 are unable to be attained for the sparse models (i.e. for (b), (c), (d)), requiring probing of higher maximum values of D g .Within the ranges present, the behavior is consistent with Figure 11.

Figure 3 :
Figure 3: Plot visualizing the extraction of λ f it in all-to-all SYK (q = 4, β = 1, 116430 disorder realizations) from the double commutator C(t) via a numerical derivative.(a) Plots of the C(t) versus time for N = 16 to N = 36.To calculate the numerical derivative, a horizontal cross section F * must be chosen, here F * = 0.4.The symmetry (labelled within the grey box) describes at what times subsequent curves should intersect F * .(b) After the numerical derivative between subsequent curves is calculated via 4.3, λ f it is plotted versus 1/N and fitted to a 2nd order polynomial (in grey box) to extract the large N Lyapunov exponent, here λ L = 0.79.

1 R 2 Figure 5 :
Figure 5: Many realizations are subdivided into groups sized D g , from which a Lyapunov exponent and R 2 is extracted.This process is repeated D tot /D g times, the average values of λ L (blue) and R 2 (green) are plotted, as well as the standard deviation of λ L , σ (orange).Resolution of plots all have ∆D g = 50.(a) all-to-all model.Exhibits convergence for plotted quantities; R 2 → 1 shows agreement with the symmetry.(b,c,d) k = 50, 10, 4 respectively.For similar values of R 2 , display similar values of error (even compared to (a)), however fluctuations increase as compared to all-to-all.

Figure 6 :
Figure 6: R 2 for the plots in Figure 5.As the model becomes more sparse, the number of realizations has to increase considerably to obtain reliable results.As in previous plots, N = 16 to 36, q = 4, β = 1.

Figure 7 :
Figure 7: Many realizations are subdivided into groups sized D g , from a Lyapunov exponent and R 2 is extracted.This process is repeated D tot /D g times, the average values of the λ L (blue) and R 2 (green) are plotted, as well as the standard deviation of λ L (orange).Resolution of plots all have ∆D g = 10.All plots are of k=4 sparse model, and have the same max D g .(a) Max value N = 36, scaled down version of Figure 5 (d), baseline of comparison.(b,c,d) Max value N = 38, 40, 42 respectively.All exhibit greater fluctuations and no significant change to R 2 for this range of D g .

Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :
Figure 8: Data analysis of Lyapunov extraction for a fixed hypergraph seed (meaning for each N , the Hamiltonian interaction configuration is fixed for each disorder realization).Random variations persist through random states used in each disorder realization.Extraction was performed on a sparse model with k = 4 between N = 16 to N = 36 with β = 1.All plots demonstrate high variance and lack of convergence, suggesting the numerical derivative method of extraction is unsuitable for a fixed hypergraph seed.Each D g is averaged over 10 times.(a) The mean fitted Lyapunov exponent λ L versus D g (b) The standard deviation of the fitted λ L 's, σ, versus D g (c) The mean R 2 value over each set of λ L 's.The very low values of R 2 suggest a very poor fit to the data.

Figure 13 :
Figure 13: A replication of Figure 7 (see corresponding caption for plot details) with changes in the lower bounds of N used to extract the Lyapunov exponent, namely (a) N = 16 to N = 36 (unchanged) (b) N = 18 to N = 38 (c) N = 20 to N = 40 (d) N = 22 to N = 42.Despite a majority of data overflowing the plot, the choice of ranges for the axes was maintained to demonstrate even further disagreement from Figure 7.

Table 2 :
Efficiency comparisons between all-to-all SYK and various sparse models.The total time is the average total time (over 112 samples) to simulate a single disorder realization.