supp.pdf

# Sparsity independent Lyapunov exponent in the Sachdev-Ye-Kitaev model

Antonio M. García-García<EMAIL_ADDRESS>Shanghai Center for Complex
Physics, School of Physics and Astronomy, Shanghai Jiao Tong University,
Shanghai 200240, China  Chang Liu<EMAIL_ADDRESS>Shanghai Center for
Complex Physics, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai 200240, China  Jacobus J. M. Verbaarschot
<EMAIL_ADDRESS>Center for Nuclear Theory and Department
of Physics Astronomy, Stony Brook University, Stony Brook, New York 11794, USA

###### Abstract

The saturation of a recently proposed universal bound on the Lyapunov exponent
has been conjectured to signal the existence of a gravity dual. This
saturation occurs in the low temperature limit of the dense Sachdev-Ye-Kitaev
(SYK) model, $N$ Majorana fermions with $q$-body ($q>2$) infinite-range
interactions. We calculate certain Out of Time Order Correlators (OTOC) for
$N\leq 64$ fermions for a highly sparse SYK model and find no significant
dependence of the Lyapunov exponent on sparsity up to near the percolation
limit where the Hamiltonian breaks up into blocks. This suggests that in the
sparse case, the Lyapunov exponent also saturates the low-temperature bound. A
key ingredient to reaching $N=64$ is the development of a novel quantum spin
model simulation library that implements highly-optimized matrix-free Krylov
subspace methods on Graphical Processing Units (GPUs). This leads to a
significantly lower simulation time as well as vastly reduced memory usage
over previous approaches, while using modest computational resources. Strong
sparsity-driven statistical fluctuations require both the use of a vastly
larger number of disorder realizations with respect to the dense limit and a
careful finite size scaling analysis. Our results potentially broadens the
landscape of theories that may have a gravity analogue.

The exponential growth of certain out-of-time-order correlation functions
(OTOC) up to the Ehrenfest time in the semiclassical limit, at a rate given by
the leading classical Lyapunov exponent, is an early signature of quantum
chaotic dynamics. Their calculation in simple single-particle problems such as
a particle in a random potential [1] or kicked rotors [2] were landmarks in
the early development of the theory of quantum chaos. However, they are
notoriously difficult to compute quantitatively in many-body systems because
the region of exponential growth is relatively short and can be easily
overshadowed by other contributions unless the system is strictly within the
semiclassical limit.

A resurgence of interest in OTOCs in quantum chaos came quite unexpectedly
from quantum gravity. Heuristic arguments [3, 4] suggest that the dynamics of
a particle close to a black hole horizon is quantum chaotic. Later, these
ideas were put on a much firmer ground by showing that a universal bound on
the Lyapunov exponents in quantum chaotic systems at thermal equilibrium was
saturated in field theories with a gravity dual [5]. Shortly afterwards,
Kitaev [6] demonstrated analytically that, in the low temperature limit, this
universal bound on chaos was saturated in a simple model, now termed the
Sachdev-Ye-Kitaev model (SYK) [6, 7, 8, 7, 9, 10, 11], consisting of $N$
Majoranas [6, 12] with random $q$-body interaction in zero dimensions. The
quantum chaotic nature of the SYK model for longer time scales was confirmed
by a level statistics analysis [13, 14] and its gravity dual was identified to
be Jackiw-Teitelboim gravity [15, 16, 17].

The analytical tractability of the SYK model is one of its most appealing
features. Unfortunately, generalizations of the model with finite range [18]
or sparsified [19, 20, 21, 22, 23] interactions do not inherit this property.
This begs the question: is the saturation of the bound, that indicates the
possible existence of a gravity dual, a particularity of the dense SYK, or is
it present in more general settings? For the dense SYK, a recent numerical
calculation [24] of the Lyapunov exponent based on the Krylov subspace method
[25, 26] for up to $N=50$ on a GPU-system (and $N=60$ on a CPU-only system)
confirmed the analytical results [6, 12]. An important benefit of the
sparsified SYK model is that it may be easier to simulate on a quantum
computer [27, 28, 29] which potentially facilitates addressing questions that
cannot be answered with classical computers.

In this paper, we aim to calculate the Lyapunov exponent for a sparse variant
[19, 20, 21, 22, 23] of the SYK model where a large (to be defined shortly)
number of random couplings are set to zero. A key ingredient in our study is
the development of a highly-optimized GPU computing code, which implements the
Krylov-based algorithm for computing time evolution of qubit systems. This
allows us to reach up to $N=64$ Majoranas on single GPU systems. Our main
result is that the Lyapunov exponent of the sparse SYK model has no
significant dependence on sparsity, and agrees with the dense case, all the
way up to close to the percolation limit. This suggests the existence of
gravity duals in a much broader family of field theories.

Sparse SYK Model. Our Hamiltonian describes $N$ strongly interacting Majorana
fermions in zero spatial dimensions [6, 7, 30, 8, 31, 9, 32, 10, 12] with
sparse [19, 20, 21, 22, 33] random interactions of infinite range:

$H=\sum_{0\leq
i<j<k<l<N}p_{ijkl}\,J_{ijkl}\,\gamma_{i}\gamma_{j}\gamma_{k}\gamma_{l}\,.$ (1)

The Majorana operators $\gamma_{n}$ satisfy the Clifford algebra
$\\{\gamma_{m},\gamma_{n}\\}=\delta_{mn}$, and can be expressed as a tensor
product of Pauli matrices. The $J_{ijkl}$ are random numbers with a Gaussian
distribution of zero average and variance $\langle
J^{2}\rangle={3}/({8pN^{3}})$. The sparseness of the Hamiltonian is modeled by
the stochastic variable $p_{ijkl}$ which is sampled from the Bernoulli
distribution $B(p)$ with probability $p>0$. When $p=1$, we recover the dense
SYK model. Models with $0<p<1$ are called sparse SYK models. In principle,
$p\sim N^{-\alpha}$ with parameter $\alpha>0$. It was shown in Ref. [20, 19],
$\alpha=3$ is the relevant scaling to study the effect of sparsity because for
$\alpha>3$ connectivity in Fock space is broken in the large $N$ limit while
the effect of sparsity is largely irrelevant for $\alpha<3$. Therefore, it is
natural to define the sparsity strength $k=\frac{p}{N}\binom{N}{4}$. For the
comparison with the dense case, we will focus on the $k\geq 3$ region only
because, for sufficiently large $N$, it is computationally expensive to impose
a regularity condition on the vertex connectivity. The latter is necessary for
$k\sim 1$ in order to prevent the Hilbert space from splitting into separate
invariant subspaces of the Hamiltonian.

OTOC calculation and results. We now define the following regularized out-of-
time-order correlation (OTOC) function for the Hamiltonian Eq.(1),

$\displaystyle F(t)$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\textrm{Tr}\left[e^{itH}\gamma_{\scriptscriptstyle
N-1}e^{-(it+\beta/4)H}\gamma_{\scriptscriptstyle N-2}\right.$ (2)
$\displaystyle\left.e^{(it-\beta/4)H}\gamma_{\scriptscriptstyle
N-1}e^{-(it+\beta/4)H}\gamma_{\scriptscriptstyle N-2}e^{-\beta H/4}\right]$

with $Z=\textrm{Tr}\left(e^{-\beta H/4}\gamma_{\scriptscriptstyle
N-1}e^{-\beta H/4}\gamma_{\scriptscriptstyle N-2}e^{-\beta
H/4}\gamma_{\scriptscriptstyle N-1}e^{-\beta H/4}\right.$
$\left.\gamma_{\scriptscriptstyle N-2}\right)$ so that $F(0)=1$. Different
regularizations may lead to slightly different prefactors in the $1/N$
expansion of the OTOC which may be time dependent but the Lyapunov exponent
was recently shown [34, 35] to be independent of the regularization.

Figure 1: The OTOC Eq. (2) versus time $t$ for $N=14,\;22,\;30,\;38,\;46$ and
54 as indicated in the legend. Results are given for $\beta=0.5$ (upper),
$\beta=5$ (middle) and $\beta=20$ (lower) and sparseness parameter $k=3$
(dotted), $k=6$ (dashed), $k=9$ (dot-dashed) and dense (solid). The dependence
on $k$ in the region where OTOC decreases exponentially is quite weak.

A strong hint of what will be the main result of this paper, the independence
of the Lyapunov exponent on the sparsity, can already be seen from Fig. 1
depicting the OTOC dependence on the sparseness parameter $k$ for different
$N\in[14,54]$: in the low-temperature limit, the OTOC depends only weakly on
$k$ even for relatively small $k=3$ close to the percolation limit where the
Hamiltonian breaks up into blocks for most disorder realizations. We restrict
ourselves to $N{\rm mod8}=6$ because the approach to the large $N$ limit may
be more uniform taking advantage of the Bott periodicity of the SYK model
[36].

GPU-based Numerical Optimizations. The calculation of the OTOC in the
large-$N$ limit of interest requires the development of a novel quantum time
evolution library [37] called REAPERS, short for “a REAsonably PERformant
Simulator for qubit systems”, that implements highly-optimized matrix-free
[38] Krylov subspace methods [25, 26] on nVidia GPUs. Written in C++20, it
provides a programming interface similar to that of dynamite [39], but does
not depend on low-level libraries such as the Portable, Extensible Toolkit for
Scientific Computation (PETSc) [40, 41, 42] for matrix operations. Instead, we
start from scratch and implement optimized Computer Unified Device
Architecture (CUDA) kernels that compute spin operator actions on quantum
states more efficiently than those provided by PETSc, and carefully manage
object allocations and de-allocations to use as little video memory as
possible. The result of these performance optimizations is that we are able to
simulate sparse SYK systems with $N=62$ fermions in double precision floating
point and $N=64$ fermions in single precision, on single-GPU systems with the
80GB version of the nVidia A100 graphics card. Comparing this with the
previous state of the art CPU-based calculation of an $N=60$ (albeit dense)
SYK system using a 500-node supercluster [24], our hardware cost is far less
and we consume far less energy. We refer the reader to the Supplemental
Materials for the technical details of our optimization techniques and further
benchmark data.

Finite Size Scaling Analysis. The expectation for quantum chaotic systems [1,
2] is that for sufficiently short times below the Ehrenfest time the OTOC
decreases exponentially. At least for low temperatures, an analysis based on
the Schwarzian action [43, 44, 45, 46] in the dense case shows that the decay
only remains exponential up to around the Ehrenfest time, after which it
approaches zero with a decreasing exponent before finally turning into a
power-like decay for very long times. In the exponentially decaying domain,
the dependence on $t$ and $N$ is only through the combination
$\exp(\lambda_{L}t)/N$ [12, 47, 48] so that

$\displaystyle F(t)=g_{0}(t)-g\left(\frac{e^{\lambda_{L}t}}{N}\right)$ (3)

with

$\displaystyle
g\left(\frac{e^{\lambda_{L}t}}{N}\right)=c_{1}\frac{e^{\lambda_{L}t}}{N}+c_{2}\left(\frac{e^{\lambda_{L}t}}{N}\right)^{2}+\cdots\;$
(4)

and $\lambda_{L}$ the Lyapunov exponent. In order to extract $\lambda_{L}$, we
have to restrict the numerical calculation of $F(t)$ to the region in which
the decay obeys Eq. (4). For that purpose, we largely follow the method of
Ref. [24] for the dense SYK based on the rescaling symmetry, $t\to
t+\log{r}/\lambda_{L}$ and $N\to Nr$ where $r>0$. In a first step, we
determine the time $t^{*}(N)$ for which $F(t)$ drops to a certain value
$F_{0}<1$. The value of $F_{0}$ cannot be too large because that would not
capture the exponential growth but it also cannot be too small because the
OTOC no longer decays exponentially. We shall see that for values of $F_{0}$
between $0.75$ and $0.85$ the results are consistent with an exponential
growth. For the scaling behavior Eq. (4) we find to leading order in $1/N$,

$\displaystyle t^{*}(N)=\frac{\log
N}{\lambda_{L}}+\frac{1}{\lambda_{L}}g^{-1}((g_{0}(t^{*}(N))-F_{0})/c_{1}).$
(5)

The rescaling symmetry requires that $g_{0}$ does not depend on $t$ in the
region of exponential decay. We will see in the Supplemental Material that
$g_{0}(t)$ depends only weakly on $t$. In that case, the second term can be
eliminated by differentiating with respect to $N$, resulting in

$\displaystyle\lambda_{L}=\frac{1}{Ndt^{*}(N)/dN}+O(1/N).$ (6)

In principle, the Lyapunov exponent can be obtained from the slope of
$t^{*}(N)$ versus $\log N$, but the slope has a residual $N$-dependence in the
time and size window at our disposal. Ideally, we fit observables for which
this residual $N$-dependence is minimized. In agreement with [24], our
numerical results suggest that the $1/N$ dependence of the inverse slope is
close to linear at low temperatures ($\beta\geq 5$) so that the Lyapunov
exponent is determined by

$\displaystyle\frac{1}{Ndt^{*}(N)/dN}=\lambda_{L}+\frac{\alpha_{1}}{N}+\frac{\alpha_{2}}{N^{2}}+O(1/N^{3})$
(7)

with $\alpha_{2}=0$. At high temperatures (say $\beta=0.5$), the $1/N$
dependence is fitted by a quadratic dependence with $\alpha_{1}=0$ except in
the dense case when the data are sufficiently accurate to use a three
parameter fit. An estimate for the Lyapunov exponent is given by the
extrapolation of $1/({Ndt^{*}(N)/dN})$ to $1/N\to 0$. Details of the fitting
procedure are left to the Supplemental Material.

Figure 2: The dependence of $t^{*}(N)$ as a function of $\log N$ for
$\beta=0.5$ (upper), $\beta=5$ (middle) and $\beta=20$ (lower) and values of
$k$ as indicated in the legends of the figures. For the numerical fit we
employ Eq. (8). For $\beta=20$, $k=3$ and $N=14$, 13 outliers that differ by
more than 30 standard deviations from the median have been excluded from the
averages [49].

Figure 3: Top: The Lyapunov exponent, obtained from $t^{*}(N)$ (bottom) and the derivative of $t^{*}(N)$ (top), as a function of the inverse sparsity parameter $k$ for $\beta=0.5$, $5$, and $20$ using two different cutoffs $F_{0}=0.75$ (dashed and open circles), $0.85$ (solid and solid disks). The dense case is represented by the point $1/k=0$. The lines stand for linear fits to the data points. See Table 1 for a summary of results. For $\beta=20$, the slope of the line which accounts for the $k$ dependence is consistent with zero.  $\lambda_{L}$ | $F_{0}=0.75,\;1/(N\partial_{N}t^{*}(N))$ | $F_{0}=0.85,\;1/(N\partial_{N}t^{*}(N))$ | $F_{0}=0.75,\;t^{*}(N)$ | $F_{0}=0.85,\;t^{*}(N)$ | Dense SD
---|---|---|---|---|---
$\beta=0.5$ | $0.86(3)-0.52(15)/k$ | $0.93(4)-0.35(20)/k$ | $0.84(2)-0.29(10)/k$ | $0.94(2)-0.20(7)/k$ | 0.86
$\beta=5$ | $0.52(1)-0.34(9)/k$ | $0.58(2)-0.19(14)/k$ | $0.50(1)-0.22(4)/k$ | $0.60(1)-0.22(6)/k$ | 0.59
$\beta=20$ | $0.33(4)-0.2(2)/k$ | $0.39(6)-0.2(3)/k$ | $0.30(2)-0.16(10)/k$ | $0.40(3)-0.2(2)/k$ | 0.24
Table 1: The Lyapunov exponents obtained by a linear fit to the data of Fig.
3. In the first two columns, $\lambda_{L}$ is obtained from Eq. (8), in the
next two columns $\lambda_{L}$ is obtained from Eq. (7). Dense SD stands for
the analytic result [12].

Alternatively, one can integrate Eq. (7), resulting in the equivalent
expansion up to logarithmic factors,

$\displaystyle t^{*}(N)=\gamma_{0}+\frac{\log
N}{\lambda_{L}}\left(1+\frac{\gamma_{1}}{N}+\frac{\gamma_{2}}{N^{2}}\right)+O\left(\frac{1}{N^{3}}\right)$
(8)

where $\lambda_{L},\;\gamma_{0}$, $\gamma_{1}$ and $\gamma_{2}$ are fitting
parameters. For $\beta\geq 5$ we set $\gamma_{2}=0$, but for $\beta=0.5$ we
use $\gamma_{2}/N^{2}$ as the correction term putting $\gamma_{1}=0$. As can
be seen from Fig. 2, this gives an excellent fit of the $t^{*}(N)$ data for
all considered temperatures, sparsity parameters and cut-off values. Fitting
$t^{*}(N)$ directly has the advantage that the errors are smaller. On the
other hand, fitting the numerical derivative $\partial_{N}t^{*}(N)$ has the
benefit of having one less fitting parameter at the expense of much larger
errors (see Supplemental Materials). We shall see the fitting results of both
methods are consistent though there is a significant systematic error. There
is also an issue of over-fitting which trades the $\log N$ dependence for the
$1/N$ dependence. For example using additive $1/N$ corrections instead of
multiplicative $1/N$ correction significantly changes the value of the
Lyapunov exponent.

Results for the Lyapunov Exponent. Our results for the Lyapunov exponents
$\lambda_{L}$ are shown in Fig. 3. In the range of temperatures we have
considered, $\beta=0.5,$ $5$, $20$, there is no significant dependence of the
Lyapunov exponent on the sparsity parameter $k$. The difference of
$\lambda_{L}$ for different values of $k$ is less than the discrepancies
between the two choices of $F_{0}$ and the two choices of the fitting methods,
both of which are a measure of the systematic error. Numerical results for the
fitted lines are given in Table 1.

We now compare our results with predictions for the dense case, either
analytical or based on the numerical solutions of the Schwinger-Dyson (SD)
equations. As shown in the last column of Table 1, for $\beta=0.5$ and
$\beta=5$ our numerical results are in good agreement with the large-$N$
prediction obtained by solving the SD equations [24]. For $\beta=20$, we find
a Lyapunov exponent that is consistent with the chaos bound of $2\pi/\beta$
[5]. Taking into account subleading finite temperature corrections (which can
be obtained by solving the SD equations) lowers the theoretical large-$N$
value to 0.24, which is still above our result but is in agreement with
previous numerical calculations in the dense case [24], where for $\beta=17.8$
the Lyapunov exponent was 0.36 versus 0.26 from solving the SD equations.
Therefore this discrepancy is not related to the sparsity of the model. We
stress that in order to reach this relatively low level of statistical
fluctuations, it is necessary to simulate a number of disorder realizations at
least of order $10^{4}$, which is several orders of magnitude larger than in
the dense case [33]. Surprisingly, unlike the dense case, the fluctuations are
not larger for low temperature and they are not reduced as $N$ increases,
which prevents us from including $N=62$ in this analysis despite the fact that
it is numerically accessible.

Conclusions and outlook. We have studied out-of-time-order correlators in a
sparse variant of the SYK model. After some careful data analysis, we have
shown that the Lyapunov exponent has no significant dependence on the sparsity
$k$ for all temperatures we have considered, and agrees with previous [24]
numerical results for the dense case. In the low temperature limit, the value
of the Lyapunov exponent for both the dense and the sparse SYK is above the
analytical prediction [12], which prevents us from explicitly confirming that
the bound on chaos [5] is still saturated for the sparse SYK model. However,
the independence of the Lyapunov exponent on the sparsity indicates that the
low-temperature bound will also be saturated in the sparse case.

A crucial part of our work is the development of an optimized GPU-based
quantum simulation library, which enables us to reach $N\leq 64$ Majoranas due
to drastic improvements in simulation speed and memory usage. The energy
consumption and the cost of hardware are both vastly smaller than equivalent
simulations on CPU-based systems. Natural extensions of this work includes
computing OTOCs in non-Hermitian SYK and sparse spin chains, such as those
employed in studies of many-body localization [50, 51].

###### Acknowledgements.

AMGG and CL were partially supported by a National Key R$\&$D Program of China
(Project ID: 2019YFA0308603), and a Shanghai talent program. JJMV acknowledges
support from U.S. DOE Grant No. DE-FAG88FR40388.

## References

  * Larkin and Ovchinnikov [1969] A. Larkin and Y. N. Ovchinnikov, Quasiclassical method in the theory of superconductivity, Sov Phys JETP 28, 1200 (1969).
  * Berman and Zaslavsky [1978] G. Berman and G. Zaslavsky, Condition of stochasticity in quantum nonlinear systems, Physica A: Statistical Mechanics and its Applications 91, 450 (1978).
  * Sekino and Susskind [2008] Y. Sekino and L. Susskind, Fast scramblers, Journal of High Energy Physics 10, 065 (2008).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Black holes and the butterfly effect, JHEP 03, 067, arXiv:1306.0622 [hep-th] .
  * Maldacena  _et al._ [2016a] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, Journal of High Energy Physics 08, 106 (2016a).
  * Kitaev [2015] A. Kitaev, A simple model of quantum holography (2015), kITP strings seminar and Entanglement 2015 program, 12 February, 7 April and 27 May 2015, http://online.kitp.ucsb.edu/online/entangled15/.
  * Bohigas and Flores [1971a] O. Bohigas and J. Flores, Two-body random hamiltonian and level density, Physics Letters B 34, 261 (1971a).
  * French and Wong [1970] J. French and S. Wong, Validity of random matrix theories for many-particle systems, Physics Letters B 33, 449 (1970).
  * French and Wong [1971] J. French and S. Wong, Some random-matrix level and spacing distributions for fixed-particle-rank interactions, Physics Letters B 35, 5 (1971).
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Gapless spin-fluid ground state in a random quantum heisenberg magnet, Phys. Rev. Lett. 70, 3339 (1993).
  * Benet  _et al._ [2001] L. Benet, T. Rupp, and H. A. Weidenmüller, Nonuniversal behavior of the $\mathit{k}$-body embedded gaussian unitary ensemble of random matrices, Phys. Rev. Lett. 87, 010601 (2001).
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Remarks on the sachdev-ye-kitaev model, Phys. Rev. D 94, 106002 (2016).
  * García-García and Verbaarschot [2016] A. M. García-García and J. J. M. Verbaarschot, Spectral and thermodynamic properties of the sachdev-ye-kitaev model, Phys. Rev. D 94, 126010 (2016).
  * Cotler  _et al._ [2017] J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker, D. Stanford, A. Streicher, and M. Tezuka, Black holes and random matrices, Journal of High Energy Physics 05, 118 (2017).
  * Jackiw [1985] R. Jackiw, Lower dimensional gravity, Nuclear Physics B 252, 343 (1985).
  * Teitelboim [1983] C. Teitelboim, Gravitation and hamiltonian structure in two spacetime dimensions, Physics Letters B 126, 41 (1983).
  * Almheiri and Polchinski [2015] A. Almheiri and J. Polchinski, Models of ads2 backreaction and holography, Journal of High Energy Physics 11, 1 (2015).
  * García-García and Tezuka [2019] A. M. García-García and M. Tezuka, Many-body localization in a finite-range sachdev-ye-kitaev model and holography, Phys. Rev. B 99, 054202 (2019).
  * García-García  _et al._ [2021] A. M. García-García, Y. Jia, D. Rosa, and J. J. M. Verbaarschot, Sparse sachdev-ye-kitaev model, quantum chaos, and gravity duals, Phys. Rev. D 103, 106002 (2021).
  * Xu  _et al._ [2020] S. Xu, L. Susskind, Y. Su, and B. Swingle, A Sparse Model of Quantum Holography (2020), arXiv:2008.02303 [cond-mat.str-el] .
  * Cáceres  _et al._ [2021] E. Cáceres, A. Misobuchi, and R. Pimentel, Sparse SYK and traversable wormholes, Journal of High Energy Physics 2021, 10.1007/jhep11(2021)015 (2021).
  * Cáceres  _et al._ [2022] E. Cáceres, A. Misobuchi, and A. Raz, Spectral form factor in sparse SYK models, Journal of High Energy Physics 2022, 10.1007/jhep08(2022)236 (2022).
  * Tezuka  _et al._ [2023] M. Tezuka, O. Oktay, E. Rinaldi, M. Hanada, and F. Nori, Binary-coupling sparse sachdev-ye-kitaev model: An improved model of quantum chaos and holography, Phys. Rev. B 107, L081103 (2023).
  * Kobrin  _et al._ [2021] B. Kobrin, Z. Yang, G. D. Kahanamoku-Meyer, C. T. Olund, J. E. Moore, D. Stanford, and N. Y. Yao, Many-body chaos in the sachdev-ye-kitaev model, Phys. Rev. Lett. 126, 030602 (2021).
  * N. [1931] K. A. N., On the numerical solution of equation by which are determined in technical problems the frequencies of small vibrations of material systems (1931).
  * Simoncini and Szyld [2007] V. Simoncini and D. B. Szyld, Recent computational developments in krylov subspace methods for linear systems, Numerical Linear Algebra with Applications 14, 1 (2007), https://onlinelibrary.wiley.com/doi/pdf/10.1002/nla.499 .
  * Jafferis  _et al._ [2022] D. Jafferis, A. Zlokapa, J. D. Lykken, D. K. Kolchmeyer, S. I. Davis, N. Lauk, H. Neven, and M. Spiropulu, Traversable wormhole dynamics on a quantum processor, Nature 612, 51 (2022).
  * Kobrin  _et al._ [2023] B. Kobrin, T. Schuster, and N. Y. Yao, Comment on Traversable wormhole dynamics on a quantum processor (2023), arXiv:2302.07897 [quant-ph] .
  * Jafferis  _et al._ [2023] D. Jafferis, A. Zlokapa, J. D. Lykken, D. K. Kolchmeyer, S. I. Davis, N. Lauk, H. Neven, and M. Spiropulu, Comment on ”Comment on ”Traversable wormhole dynamics on a quantum processor” ” (2023), arXiv:2303.15423 [quant-ph] .
  * Bohigas and Flores [1971b] O. Bohigas and J. Flores, Spacing and individual eigenvalue distributions of two-body random hamiltonians, Physics Letters B 35, 383 (1971b).
  * Mon and French [1975] K. Mon and J. French, Statistical properties of many-particle spectra, Annals of Physics 95, 90 (1975).
  * Benet and Weidenmüller [2003] L. Benet and H. A. Weidenmüller, Review of the k -body embedded ensembles of gaussian random matrices, Journal of Physics A: Mathematical and General 36, 3569 (2003).
  * Cáceres  _et al._ [2023] E. Cáceres, T. Guglielmo, B. Kent, and A. Misobuchi, Out-of-time-order correlators and lyapunov exponents in sparse syk (2023), arXiv:2306.07345 [hep-th] .
  * Tsuji  _et al._ [2018] N. Tsuji, T. Shitara, and M. Ueda, Bound on the exponential growth rate of out-of-time-ordered correlators, Physical Review E 98, 012216 (2018).
  * Romero-Bermúdez  _et al._ [2019] A. Romero-Bermúdez, K. Schalm, and V. Scopelliti, Regularization dependence of the OTOC. which lyapunov spectrum is the physical one?, Journal of High Energy Physics 2019, 10.1007/jhep07(2019)107 (2019).
  * You  _et al._ [2017] Y.-Z. You, A. W. W. Ludwig, and C. Xu, Sachdev-ye-kitaev model and thermalization on the boundary of many-body localized fermionic symmetry-protected topological states, Phys. Rev. B 95, 115150 (2017).
  * Liu [2023] C. Liu, github.com/cl91/REAPERS. REAPERS: a REAsonably PERformant Simulator for qubit systems (2023).
  * Coppersmith [1993] D. Coppersmith, Solving linear equations over gf(2): block lanczos algorithm, Linear Algebra and its Applications 192, 33 (1993).
  * Kahanamoku-Meyer and Wei [2023] G. D. Kahanamoku-Meyer and J. Wei, Gregdmeyer/dynamite: v0.3.0 (2023).
  * Balay  _et al._ [1997] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, Efficient management of parallelism in object-oriented numerical software libraries, in  _Modern Software Tools for Scientific Computing_, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen (Birkhäuser Boston, Boston, MA, 1997) pp. 163–202.
  * Hernandez  _et al._ [2005] V. Hernandez, J. E. Roman, and V. Vidal, Slepc: A scalable and flexible toolkit for the solution of eigenvalue problems, ACM Trans. Math. Softw. 31, 351–362 (2005).
  * Roman  _et al._ [2022] J. E. Roman, C. Campos, L. Dalcin, E. Romero, and A. Tomas,  _SLEPc Users Manual_ , Tech. Rep. DSIC-II/24/02 - Revision 3.18 (D. Sistemes Informàtics i Computació, Universitat Politècnica de València, 2022).
  * Maldacena  _et al._ [2016b] J. Maldacena, D. Stanford, and Z. Yang, Conformal symmetry and its breaking in two-dimensional nearly anti-de sitter space, Progress of Theoretical and Experimental Physics 2016, 12C104 (2016b).
  * Bagrets  _et al._ [2017] D. Bagrets, A. Altland, and A. Kamenev, Power-law out of time order correlation functions in the syk model, Nuclear Physics B 921, 727 (2017).
  * Bagrets  _et al._ [2016] D. Bagrets, A. Altland, and A. Kamenev, Sachdev–ye–kitaev model as liouville quantum mechanics, Nuclear Physics B 911, 191 (2016).
  * Altland and Sonner [2021] A. Altland and J. Sonner, Late time physics of holographic quantum chaos, SciPost Physics 11, 034 (2021).
  * Lam  _et al._ [2018] H. T. Lam, T. G. Mertens, G. J. Turiaci, and H. Verlinde, Shockwave S-matrix from Schwarzian Quantum Mechanics, JHEP 11, 182, arXiv:1804.09834 [hep-th] .
  * Yang [2019] Z. Yang, The Quantum Gravity Dynamics of Near Extremal Black Holes, JHEP 05, 205, arXiv:1809.08647 [hep-th] .
  * Rousseeuw and Hubert [2011] P. J. Rousseeuw and M. Hubert, Robust statistics for outlier detection, Wiley interdisciplinary reviews: Data mining and knowledge discovery 1, 73 (2011).
  * Basko  _et al._ [2006] D. Basko, I. Aleiner, and B. Altshuler, Metal–insulator transition in a weakly interacting many-electron system with localized single-particle states, Annals of Physics 321, 1126 (2006).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field heisenberg chain, Phys. Rev. B 91, 081103 (2015).

in 1,…,0 See pages \x, of supp.pdf