Many-Body Chaos in the Sachdev-Ye-Kitaev Model

Many-body chaos has emerged as a powerful framework for understanding thermalization in strongly interacting quantum systems. While recent analytic advances have sharpened our intuition for many-body chaos in certain large $N$ theories, it has proven challenging to develop precise numerical tools capable of exploring this phenomenon in generic Hamiltonians. To this end, we utilize massively parallel, matrix-free Krylov subspace methods to calculate dynamical correlators in the Sachdev-Ye-Kitaev (SYK) model for up to $N = 60$ Majorana fermions. We begin by showing that numerical results for two-point correlation functions agree at high temperatures with dynamical mean field solutions, while at low temperatures finite-size corrections are quantitatively reproduced by the exactly solvable dynamics of near extremal black holes. Motivated by these results, we develop a novel finite-size rescaling procedure for analyzing the growth of out-of-time-order correlators (OTOCs). We verify that this procedure accurately determines the Lyapunov exponent, $\lambda$, across a wide range in temperatures, including in the regime where $\lambda$ approaches the universal bound, $\lambda = 2\pi/\beta$.

Many-body chaos has emerged as a powerful framework for understanding thermalization in strongly interacting quantum systems. While recent analytic advances have sharpened our intuition for many-body chaos in certain large N theories, it has proven challenging to develop precise numerical tools capable of exploring this phenomenon in generic Hamiltonians. To this end, we utilize massively parallel, matrix-free Krylov subspace methods to calculate dynamical correlators in the Sachdev-Ye-Kitaev (SYK) model for up to N = 60 Majorana fermions. We begin by showing that numerical results for two-point correlation functions agree at high temperatures with dynamical mean field solutions, while at low temperatures finite-size corrections are quantitatively reproduced by the exactly solvable dynamics of near extremal black holes. Motivated by these results, we develop a novel finite-size rescaling procedure for analyzing the growth of out-of-time-order correlators (OTOCs). We verify that this procedure accurately determines the Lyapunov exponent, λ, across a wide range in temperatures, including in the regime where λ approaches the universal bound, λ = 2π/β.
Understanding the origin of thermalization in strongly interacting quantum systems is of central interest in condensed matter, quantum information, and quantum gravity. Recent developments towards this goal have revealed striking insights on the relation between quantum chaos and the delocalization, or scrambling, of quantum information. This unification is provided by out-of-timeorder correlators (OTOCs), which take the general form W (t)V (0)W (t)V (0) for local operators V and W [1][2][3]. From an information theoretic perspective, these correlators determine the degree to which local information becomes hidden in nonlocal degrees of freedom, leading to the effective memory loss of initial conditions [2,4]. From the perspective of chaos, OTOCs measure the sensitivity of one operator towards a small perturbation induced by another operator at an earlier time [5,6]. In particular, for semiclassical chaotic systems, OTOCs are expected to exhibit a period of exponential growth analogous to the classical butterfly effect [7,8].
At the intersection between these two perspectives lies the discovery of a new form of quantum chaos in strongly interacting systems, known as many-body chaos. This phenomenon is characterized by OTOCs whose leading order behavior is given by e λt /N , where N is related to the number of degrees of freedom per site [5,9]. While such behavior was first anticipated in [9] and confirmed using holographic duality in [2], the first concrete model that exhibited many-body chaos was introduced by Kitaev following previous work by Sachdev and Ye [10][11][12][13]. Remarkably, at low temperatures, the so-called SYK model saturates a universal bound, λ ≤ 2πT , where T is the temperature of the system and we henceforth set k B = = 1 [5]. The saturation of this bound is known to occur in theories of quantum gravity and their holographic duals [6], and indeed a direct correspondence has since been established between the low temperature dynamics of the SYK model and a universal theory of near extremal black holes (i.e. Jackiw-Teitelboim gravity) [13][14][15][16]. More recently, a number of other models that exhibit many-body chaos have been studied; however, their rate of chaos is parametrically slower than the thermodynamic bound [17,18]. In parallel, there have been numerous proposals to measure OTOCs in coherently controlled experiments [19][20][21][22][23] and various suggestions for models that can be realized in these settings [8,19,[24][25][26][27][28].
FIG. 1. Typical data for regularized OTOCs,F (t) ≡ F (t)/F (0), as shown for βJ = 10 and system sizes N ∈ [12,60]. The early-time behavior is characterized by ∼ e λt and different system sizes are approximately related by a time translation symmetry, t → t + 1/λ log N . (b) Applying a finite-size rescaling procedure to the data, we determine λ as a function of temperature (points). Our results show excellent agreement with the theoretical predictions of the Schwinger-Dyson (SD) equations (dashed line), including in the regime where λ approaches 2π/β (blue).
A major difficulty with benchmarking these proposals and identifying further many-body chaotic models is the lack of reliable numerical tools to diagnose manybody chaos. In particular, to observe a period of clear exponential growth, the scrambling time must be wellseparated from other effects related to local relaxation that occur at early times [5]. However, for relatively small systems, there may be little or no separation between these timescales, or from late-time effects that occur after the scrambling time [24,29]. Understanding how to increase the separation between these timescales, as well as to account for the non-exponential effects, is thus crucial for numerical studies and future experimental realizations of many-body chaotic systems.
Here, we take steps to overcome these challenges by employing massively parallelized Krylov subspace methods and developing new extrapolation tools to characterize many-body chaos. Specifically, we compute correlation functions for the SYK model for systems of up to 60 Majorana fermions and leverage the model's correspondence with quantum gravity to interpret finite-size effects. We present three main results. First, we demonstrate that our numerical results for two-point functions, G(t) = W (t)W (0) , agree quantitatively with analytic predictions in two distinct regimes: (i) at high temperatures, our results match the mean-field solution of the microscopic model, and (ii) at low temperatures, our results are consistent with the full quantum dynamics of near extremal black holes. These latter results represent, to the best of our knowledge, the first direct numerical verification of quantum gravity correlators, and highlight the close connection between finite-size corrections and gravitational fluctuations.
Second, we introduce an extrapolation procedure for determining λ that takes into account higher-order terms in OTOCs. This procedure exploits an emergent rescaling symmetry for OTOCs, t → t + λ −1 log N , that applies to a large class of many-body chaotic systems close to a semiclassical limit. We verify that this procedure accurately determines λ as a function of temperature, including at low temperatures where λ ≈ 2πT (Fig. 1).
Third, we discuss the importance of the magnitude C 1 of the leading exponentially growing term in the OTOC, (C 1 /N )e λt [30]. This magnitude determines the separation between the scrambling time and the early-time, dissipative dynamics.
In particular, we find that the temperature dependence of the magnitude varies significantly between two types of OTOCs [7,23,31]: the regularized configuration, 1 4 , with ρ = e −βH , which is the quantity typically computed in field theory; and the unregularized configuration, F (u) (t) = W (t)V (0)W (t)V (0)ρ , which is more easily accessible in experiments. As the scrambling time is generally shorter in the unregularized case, this correlator is subject to more significant finite-size effects, which we corroborate through our numerics. Our analysis further leads to a mathematical proof demonstrating that λ is equal for the two configurations for the SYK model and most other known many-body chaotic systems. Thus, although λ is the same in the two configurations, it is more difficult to extract from the unregularized OTOC.
The SYK model and its gravity dual.-We begin with a brief overview of the SYK model, highlighting prior analytical results that will be relevant for our numerical study. Consider the SYK Hamiltonian given by [11,12]: Here χ i (i = 1, . . . , N ) are Majorana fermions which obey the anti-commutation relation, {χ i , χ j } = δ ij . J ijkl are random (real) coefficients sampled from a Gaussian distribution characterized by a standard deviation σ = J, which sets the microscopic energy scale in the system. We are interested in probing the out-of-equilibrium dynamics with two different types of correlators. In-time correlators reveal how excitations in the system relax towards equilibrium. In particular, we will consider the average imaginary-time Green's function, G(τ ), and its real-time cousin, G R (t), given by where τ (t) > 0 is imaginary (real) time, · · · β = 1 Z Tr · · · e −βH is a thermal average at inverse temperature β = 1/T , and the overline denotes the (quenched) average over disorder realizations. On the other hand, to probe chaos and scrambling of quantum information, we will consider out-of-time-order correlators. We will primarily focus on the regularized OTOC, where i = j, and ρ = e −βH , the imaginary-time evolution associated with the thermal ensemble, is distributed evenly among the four operators. Eventually, we will consider the difference between this correlator and the unregularized version, i.e. χ i (t)χ j (0)χ i (t)χ j (0)ρ .
In the large N , semiclassical limit, both in-time and out-of-time correlators can be exactly computed via a diagrammatic approach [11,12]. The average Green's functions are determined by the self-consistent Schwinger-Dyson equations. For the OTOCs, the leading order term in 1/N is computed by summing a series of diagrams known as ladder diagrams. This reveals a temperature-dependent Lyapunov exponent that, at low temperatures, approaches the universal bound, λ ≤ 2π/β ( Fig. 1(b)).
Beyond the semiclassical limit, the dynamics at low temperatures (i.e. βJ 1) is captured by an effective theory known as the "Schwarzian theory" [11,13,32,33]. The same theory also describes Jackiw-Teitelboim (JT) gravity, a simple quantum gravity description of two-dimensional Anti-de-Sitter space; this is the basis for the connection between our numerical calculations and quantum gravity.
Crucially, correlators in the Schwarzian theory are exactly computable [33,34], which will enable us to make quantitative finite-size comparisons for two-point functions G(τ ) and G R (t) outside of the semiclassical limit. However, for the four-point function, the expressions are more complicated, and we will compare numerics to the following ansatz This is valid for large N and t 1/λ log N [15,33,34]. An analogous series expansion is expected to characterize OTOCs for the SYK model at high temperatures (and any other model described by ladder diagrams) [17,30]; it is this series expansion that will later motivate our extrapolation procedure for determining λ. Numerical methods.-Having established the theoretical framework for understanding dynamics in the SYK model, we now discuss our numerical tools for computing correlation functions for finite-size systems. Our approach relies on a class of iterative methods known as Krylov subspace methods. These methods approximate action of the unitary operator U (t) = e −Ht by successively multiplying the Hamiltonian to an initial state.
Such matrix-vector multiplications are amenable to massive parallelization, allowing us to simulate systems of up to 60 Majorana fermions [35].
One caveat with our approach is that calculating a thermal average exactly with Krylov subspace methods would require a separate computation for each state in the Hilbert space. Fortunately, we can overcome this limitation by approximating the thermal average with a typical state [36]: where |ψ is a Haar-random state, and we have distributed the Boltzmann factor evenly to reduce numerical errors. For thermalizing systems (i.e. those that obey the Eigenstate Thermalization Hypothesis), the error in this approximation is exponentially small with respect to the system size and can be reduced further by averaging over initial states. In practice, we average simultaneously over initial states and disorder realizations of the coefficients in the Hamiltonian, J ijkl . We verify that this procedure yields the same averaged correlation functions as exact diagonalization within accessible system sizes [37]. Two-point functions-To begin probing the thermalizing dynamics of the SYK model, we compute the average Green's functions for both real-and imaginary-time evolution in the temperature range, 0 < βJ ≤ 100. Our results verify quantitative predictions in the regimes of The disagreement at earlier times is attributed to the difference in high-energy modes, which are cut off at the energy scale J in the SYK model and are unbounded in the effective action. (Inset) A salient feature in our real-time numerics is a non-monotonic trend with respect to temperature, as shown for tJ = 20. This behavior is captured by the Schwarzian action (dashed) and can be understood as a consequence of the square root edge of the energy spectrum.
analytical control (Fig. 2), as well as reveal effects due to the microscopic nature of the model that are difficult to calculate analytically.
In Fig. 3(a), we present numerical data for the imaginary-time Green's function for a system of 44 Majoranas. At high temperatures (βJ < 5), our data show excellent agreement with the semiclassical solution given by the Schwinger-Dyson equations. This implies that the high-energy dynamics are well-described by the semiclassical limit.
At lower temperatures, the difference between our numerics and the semiclassical solutions widens. To understand the origin of these corrections, we plot the full solution predicted by the Schwarzian action. This exhibits close quantitative agreement with our data at temperatures corresponding to βJ 50. We thus confirm that the Schwarzian action, or its corresponding gravity dual, accurately captures finite-size corrections away from the semiclassical regime.
A few additional remarks are in order. First, we note that the agreement with the Schwarzian action is only valid for system sizes larger than N ≈ 30 [37]. For smaller sizes, we observe additional finite-size corrections that are attributed to the discreteness of the energy spectrum. Such non-Schwarzian corrections are expected to dominate when the temperature approaches the energy of the level spacings, which corresponds to N ∼ log β (Fig. 2) [38,39]. Second, the agreement between the Schwarzian and our numerics does not hold at timescales shorter than the inverse of the microscopic coupling strength (i.e. τ J 1); specifically, the Schwarzian dynamics diverge as τ J → 0 while our numerics approach a finite value. This difference arises from the fact that the Scharzian action is the effective theory only at low energies (compared to J); for higher energies, the SYK dynamics are governed by the microscopic nature of the model.
We next turn to the computation of the retarded Green's function. As in the imaginary-time case, we observe that our data agree with the semiclassical solutions at high temperatures and with the full dynamics of the Schwarzian action at low temperature ( Fig. 3(b)). We note, however, that the early-time discrepancy with the Schwarzian action is extended to later times (βJ ∼ 10). This can be attributed to the longer timescale required for the phase cancellation of the high-energy modes in real time, as opposed to the direct suppression that occurs in imaginary time.
Working in real time further allows us to probe effects at later times compared to imaginary-time evolution, which was bounded by 0 < τ < β. Based on random matrix theory, we expect the late-time dynamics to be governed by the functional form of the spectral density at low energies, which is given by ρ(E) ∼ E 1 2 [32,40,41]. In particular, this square-root singularity leads to a powerlaw decay of the Green's function, with a power that de- The Lyapunov exponent λ is determined from the slope of t * (N ) with respect to log N . In practice, we approximate the slope for successive system sizes and extrapolate to N → ∞ [37].
pends on both the temperature and the timescale. Moreover, it implies a non-monotonic temperature dependence for the decay of the Green's function, in stark contrast to the monotonic dependence predicted by the semiclassical solution. Examining our data at late times, we verify such non-trivial temperature dependence consistent with the full Schwarzian solution ( Fig. 3(b) inset). This further corroborates the validity of the Schwarzian action at describing the full quantum dynamics in the low temperature regime.
Extracting λ.-To probe many-body chaos, we now compute regularized OTOCs (Eq. 4) for temperatures in the range 0 < βJ ≤ 56. In the large N limit, one expects a well-defined period of exponential growth, starting from the timescale at which the two-point functions decay and persisting until the scrambling time [5]. However, for the system sizes accessible in our numerics, we find that there is little separation between these timescales. As a result, fitting our data to an exponential yields poor agreement with the expected Lyapunov exponent [37].
By contrast, we determine that a novel extrapolation method provides a robust way of extracting the Lyapunov exponent. The intuition behind this method is as follows: For a large class of many-body chaotic systems, the full form of the OTOC in the semiclassical limit is given by a series in e λt /N (Eq. 5). Crucially, this series exhibits a rescaling symmetry such that scaling N → rN amounts to shifting the full curve by t → t + 1/λ log r. This symmetry can be shown explicitly for the Scharzian action, which governs the low-temperature dynamics of the SYK model, and is also expected to apply to the SYK model at high temperatures.
While the rescaling symmetry is exact in the semiclas-sical limit, we anticipate that it remains approximately valid even at finite sizes. This suggests the we can determine λ at a given temperature by attempting to collapse our data through finite-size rescaling of the form t → t + 1/λ log N . More specifically, we implement the following procedure, as illustrated in Fig. 4. First, we interpolate our data to find the time, t * , at which each curve crosses a fixed value, i.e.
where N corresponds to the system size about which we take the numerical derivative. Finally, we fit our results to a 1/N series, λ fit (N ) = λ 0 + λ 1 /N + λ 2 /N 2 + · · · ; the leading order term λ 0 corresponds to the extrapolated value for λ as N → ∞.
In Fig. 1, we present our results for λ 0 as a function of temperature. We observe excellent agreement with analytic predictions for all temperatures in the range 0 < βJ ≤ 56. This confirms our procedure as a robust method for characterizing many-body chaos.
A natural question to ask is over what range of temperatures we expect our procedure to remain valid. Our analysis indicates that there are, in fact, three relevant considerations. First, as discussed in the context of twopoint correlators, the temperature must be high compared to the energy associated with the level spacing, which corresponds to log βJ N . We account for this requirement by considering only system sizes where at least 20 eigenstates, on average, lie within ∆E = 1/β of the ground state. Second, the system must be sufficiently close to the semiclassical limit for the rescaling symmetry to approximately hold. It is known from the Schwarzian action that this condition corresponds to βJ N . Asymptotically this is a much stronger requirement than the former, yet for the system sizes relevant for our study (N ≈ 50) both requirements imply a low temperature limit of βJ ≈ 50. Third, there must be a sufficient separation between the scrambling time and shorttime dissipative dynamics. In the case of the regularized correlator, this condition is approximately βJ N , leading to the same temperature range as the semiclassical requirement. However, as we will discuss next, this condition is more stringent in the case of unregularized OTOCs.
Magnitude of OTOCs.-To understand the separation of timescales, we consider the normalized OTOC and define the magnitude of exponential growth C 1 as Having a well-separated scrambling time, t * ≈ 1/λ log(N/C 1 ), thus corresponds to N C 1 . Crucially, C 1 depends on both the temperature and regularization.
Previous work determined the lowtemperature limit of C 1 for the case of the regularized OTOC [11,15]. More recently, an identity was derived that allows for the calculation of C 1 for any model whose With regularization, we observe that the growth timescale increases significantly at low temperatures; whereas, without regularization, there is little temperature dependence between βJ = 1.8 and βJ = 56. These qualitative trends are related to the difference in the scrambling time for the two types of correlators. In particular, the scrambling time for the unregularized OTOC is highly suppressed at low temperatures, implying that the observed growth in the numerics arises from dissipative dynamics rather than to chaos. (Inset) Schematic of the two configurations, represented as a path in real (horizontal) and imaginary (vertical) time.
OTOCs are governed by ladder diagrams [30]. We apply this identity to calculate C 1 for the SYK model at all temperatures [37]. At high temperatures, the regularization scheme is irrelevant and the magnitude C 1 is given by an order one constant. At low temperatures, however, the temperature dependence between the two cases differs substantially: for the regularized OTOC C 1 ∼ βJ and for the unregularized OTOC C 1 ∼ (βJ) 3 . This implies that larger systems are necessary to have a well-separated scrambling time in the case of the unregularized correlator.
To validate this prediction, we numerically simulate unregularized OTOCs in the same temperature and size range as in the previous section. The results are shown in Fig. 5. We observe that the unregularized OTOC exhibits a much weaker dependence on temperature than the regularized OTOC. This is consistent with a lack of separation between the scrambling time and early-time effects at low temperatures; i.e. the growth is attributed to dissipative dynamics rather than to chaos. Furthermore, when we repeat our scaling analysis to extract λ using the unregularized correlator, we find that the extrapolation does not converge for temperatures below βJ ≈ 10 [37].
We emphasize that while the magnitude of the OTOC differs for the two types of regularizations, it can be shown rigorously that the exponent does not. The proof is given in Supplementary Information and applies to any many-body chaotic system whose OTOCs are governed by ladder diagrams [37]. Discussion and outlook.-In this work, we demonstrated that by employing massively parallelized Krylov subspace methods and careful extrapolation tools we can accurately capture the thermalizing and chaotic dynamics of the SYK model. Our results for two-point Green's functions represent a direct verification of the dynamics of quantum black holes in a highly fluctuating regime. Moreover, our finite-size rescaling procedure for extracting Lyapunov exponents leads to the first numerical evidence for the theoretical bound, λ ≈ 2πT . We also discussed the importance of the magnitude of OTOCs for finite-size simulations, revealing that unregularized OTOCs are generally subject to larger finite-size corrections than regularized OTOCs.
We anticipate that the numerical tools demonstrated here will open the door to a number of intriguing future directions. First, our numerical tools can be applied to variations of the SYK model (i.e. large q limit) for which the effective action (i.e. Liouville action) is known for all temperatures [41][42][43]. This will enable quantitative studies of finite-size corrections in the high-temperature regime, where the Schwarzian action is not valid. Second, our procedure for characterizing Lyapunov exponents can diagnose many-body chaos in other models beyond the SYK model; this is of particular relevance for experimental platforms that have constraints on the types of interactions and disorder that can be realized [27,44,45]. Third, we envision future numerical simulations to test more complex gravitational phenomena, including traversable wormholes [46,47]. Notably, the finite-size constraints of these simulations correspond to gravity far from the semiclassical regime and thus enable a platform for testing full quantum gravity dynamics.

References 12
NUMERICAL METHODS

Jordan-Wigner transformation
To simulate the SYK model, we represent N Majorana operators as N/2 spin-1/2 operators using the canonical Jordan-Wigner transformation: where 0 ≤ i < N and k = floor(i/2). The size of the Hilbert space is thus 2 N/2 . We note that the SYK Hamiltonian contains a single Z 2 symmetry, which leads to the conservation of the fermionic charge parity, i.e. P = ( i σ z i ) mod 2 [1]. To exploit this symmetry, we prepare initial states in one of the symmetry sectors and evolve under the relevant sector of the Hamiltonian, whose dimension is half of the full Hilbert space. We verify that this simplification has a negligible effect on correlation functions compared to initial states that span both symmetry sectors.

Benchmarking Krylov subspace methods
We estimate the numerical error of the Krylov subspace method by comparing to exact diagonalization (ED). In particular, we compute the maximum absolute difference of OTOCs between the two methods, for system sizes, N ∈ [12, 16, 20, 24] and for times, tJ ∈ [0, 20]. The results are shown in Fig. 1. In all cases, the absolute error is less than 10 −12 and is considered negligible for the purposes of this study.

Disorder fluctuations
Although the Krylov subspace error is negligible, our numerical data are subject to large fluctuations due to two sources of disorder: the Hamiltonian coefficients, J ijkl , and the initial state |Ψ 0 . Fortunately, the magnitude of both types of fluctuations is expected to decrease with system size. For the former, the SYK model is self-averaging, implying that as N → ∞ the correlation functions for a single disorder realization approach the disorder average. In particular, one expects the fluctuations to decrease as a function of the number of random J ijkl coefficients, i.e. polynomially with system size. For the latter, our method for approximating thermal averages with random pure states is expected to be accurate up to exponential corrections in the system size (Eq. 6 in main text).
With this qualitative understanding in hand, we numerically determine the magnitude of the fluctuations with respect to each type of disorder. Specifically, we calculate the fluctuations in the timescale t * such thatF (t * ) = 0.25 in two different ways: (a) by fixing J ijkl and calculating the standard deviation with respect to different initial states; and (b) by averaging first over initial states and then determining the standard deviation with respect to different realizations of J ijkl . These results are shown in Fig. 2. In general, the two types of fluctuations are on the same order of magnitude, and their magnitude increases dramatically at small system sizes and low temperatures. The size dependence is consistent with the aforementioned self-averaging behavior, and the temperature dependence is attributed to the relatively small set of low-energy states that contribute to the behavior of the low temperature correlators.
We hasten to emphasize that at the system sizes relevant for our study, both types of fluctuations are significant and extensive disorder averaging is required to obtain precise results (e.g. 100 disorder realizations even for N ≈ 40 Majoranas). As a result, while single curves have been obtained for 60 Majoranas, the primary results for this study were based on N ≤ 46 Majoranas, for which sufficient disorder averaging could be performed.

Fitting to a simple exponential
Numerous prior studies of many-body chaos have characterized Lyapunov exponents by fitting OTOCs to a simple exponential form, ∼ e λt [2][3][4][5]. In this section, we apply this fitting procedure to our numerical data and compare our results to the known theoretical values for λ. In particular, we perform least squares regression on each normalized OTOC curve,F (t) ≡ F (t)/F (0), using the fitting function, a + be λ fit t , within a range defined by F 0 ≤ 1 −F (t) ≤ F 1 . We then extrapolate the fits at different system sizes using a quadratic extrapolation, λ fit (N ) = λ 0 + λ 1 /N + λ 2 /N 2 . In Fig. 3, we illustrate this procedure and show the extrapolated values for λ using various values for F 0 and F 1 . It is clear that this approach does not converge to the exact values of λ for any temperature; indeed, the estimated values are approximately a factor of 2 smaller than the theoretical expectations.
A few remarks are in order. First, our fitting procedure differs slightly from other studies in the sense that we perform the fits for a fixed range in the magnitude of the OTOCs rather than a fixed range in time, i.e. t 0 ≤ t ≤ t 1 . We chose this approach because the growth of OTOCs occurs at different timescales, depending on the temperature and system size. To compare directly to previous work, we also tried fitting our data across a fixed range in time and found no improvements in the estimates for λ. Second, we observe that our fitting results depend sensitively on the choice of F 0 and F 1 , though for all choices of these parameters our results for λ were inconsistent with theoretical expectations.
In principle, one expects the best results using F 1 1, as the simple exponential form is only defined for the initial growth. For a more precise estimate of the range of validity, we turn to the semiclassical solution for F (t) at low temperatures, (3), which takes into account higher order terms. By plotting this full solution against the leading order exponential, we find that the exponential is a good approximation only up to F 1 ≈ 0.01. In finite-size numerics, probing such small magnitudes is challenging for several reasons. First, the absolute size of numerical fluctuations is approximately constant at all times, and thus the relative size of fluctuations compared to the signal is enhanced for small F 1 . Second, one requires the timescale for which the exponential reaches F 1 to be much longer than the dissipation time. As the former timescale scales as log F 1 N and the latter timescale is constant, achieving this separation becomes more difficult as F 1 decreases. These challenges, as well as the poor results of our exponential fits, underscore our motivation for developing a fitting procedure that takes into account higher-order terms.
Fitting to the low-temperature, semiclassical solution We next consider fitting our numerical data to semiclassical form of the OTOC at low temperatures (see Section below) [6]: where U is the confluent hypergeometric function, we set ∆ = 1/4 (as expected from the Schwarzian action), and λ fit and N are fitting parameters associated with the Lyapunov exponent and the system size, respectively. This function provides a phenomenological model for capturing higher-order effects that occur after the initial exponential growth (i.e. saturation behavior). Nevertheless, we emphasize that the exact form of the function is only rigorously justified at low temperatures where the SYK model is described by the Schwarzian action. As before, we perform least squares regression within a window defined by F 0 ≤ F (t) ≤ F 1 . We then extrapolate the fitting parameter λ using a quadratic extrapolation function. In Fig. 4, we summarize the results of this approach. In general, we find better agreement with theoretical predictions than with the previous exponential fits, especially at high temperatures. Upon closer inspection, however, it is evident that the fitted values for λ do not extrapolate to the theoretical predictions, regardless of the choice for F 0 and F 1 . We conclude that this fitting procedure is not robust against finite-size and temperature corrections away from the low-temperature, semiclassical limit where it was derived.

Finite-size rescaling method
Having ruled out the possibility of fitting our data to a simple functional form, we now introduce a model-free method for extracting the Lyapunov exponent based on finite-size scaling. The only assumption we make is that the OTOCs approximately obey a rescaling symmetry of the form, This symmetry is expected to hold for any many-body chaotic model governed by ladder diagrams close to the semiclassical limit. Based on this symmetry, we devise the following numerical procedure to extract λ. First, we compute the timescale at which the OTOCs reach a fixed valueF (t N ) = F 0 , for system size N ; this requires interpolating our numerical data and solving for the intercept at F 0 . Second, we estimate λ fit (N ) via a numerical derivative, i.e. 1/λ fit (N ) = (t N − t N −1 )/(log N − log (N − 1)). Finally, we extrapolate the derivatives to N → ∞ using a polynomial extrapolation function, e.g. λ fit (N ) = λ 0 + λ 1 /N + λ 2 /N 2 . The extrapolations for various temperatures are shown in Fig. 5.
A few comments are in order. First, the extrapolation is performed on a subset of system sizes whose lowest 20 eigenstates lie within ∆E = 1/β of the ground state. This criterion is meant to rule out systems that are dominated by the discreteness of the energy spectrum, for which no effective (replica-diagonal) action exists. Furthermore, to avoid overfitting, we use a quadratic extrapolation function for temperatures corresponding to β ≥ 5.6 and a linear extrapolation for lower temperatures. The reported error bars on our final results correspond to the standard error of the fitting parameter, λ 0 .
Second, we note that the rescaling procedure depends on the value of F 0 . In the large N limit, the choice of this parameter is arbitrary, as the rescaling symmetry (4) is expected to hold for all values of F 0 . At finite sizes,  however, there are higher-order corrections that break the rescaling symmetry, particularly (i) at early times due to the microscopic cutoff, and (ii) at late times due to the crossover to power-law decay. We thus expect an intermediate choice of F 0 to provide the best approximation. Our results in the main text are based on F 0 = 0.25. In Fig. 3, we show that a different choice, F 0 = 0.16, provides consistent results. This demonstrates that our rescaling procedure is not overly sensitive to the exact value F 0 .

Unregularized OTOCs
We repeat the scaling procedure for the unregularized correlators and compare to the results with regularization, as shown in Fig. 6. For temperatures above βJ ≈ 10, regularization has little effect on λ fit (N ) (i.e. the estimate for λ at system size N ). However, for lower temperatures, the values for λ fit (N ) using the unregularized correlator show a much weaker temperature dependence and are much further from the theoretical predictions than in the regularized case. We attribute this discrepancy to the finite-size scaling of the magnitude of the OTOC rather than to a difference in the Lyapunov exponent (see section below). In particular, this implies that for sufficiently large system sizes (i.e. N (βJ) 3 ), we would expect the estimated Lyapunov exponent for the unregularized correlator to converge on the theoretical values. Our numerics are consistent with this prediction, though the limited range of system sizes makes it difficult to validate it directly.

LARGE N SOLUTIONS
In this section, we describe the large N , semiclassical solutions for the dynamics of the SYK model. These results were derived previously via either a diagrammatic approach or from the saddle-point of a disorder averaged effective action [7][8][9][10].

Schwinger-Dyson equations
The large N solution for the average Green's functions is given in imaginary time by the self-consistent Schwinger-Dyson equations: where ω are the Fourier components with respect to τ . The real-time version of these equations is obtained by setting τ = it. At low temperatures (βJ 1) and long timescales (τ J, tJ 1), the derivative term iω can be neglected, leading to an emergent conformal symmetry. The solution in this limit is where b ≈ 0.531 [9]. More generally, the Green's functions can be computed at all temperatures by solving (5) through an iterative numerical approach [9]. This procedure converges quickly in imaginary time and yields the results shown in Fig. 2(a) of the main text. In real time, the numerical analysis is more subtle, and we found that the most stable approach is the implementation proposed in [11]. We present these results in Fig. 2(b) of the main text and in Fig. 10 below. We also rely on the real-time correlators to compute the Lyapunov exponent and magnitude of OTOCs, as described in the following section. For both real and imaginary time, we benchmarked the numerical solutions by comparing to (6) at low temperatures.

Kernel equation
The leading order behavior for OTOCs is computed via a set of diagrams known as ladder diagrams. We begin by defining and make a growth ansatz of the form where t 12 = t 2 − t 1 . The exponent is determined by solving the eigenvalue equation with eigenvalue one. Here K R is the retarded kernel given by where G R (t) is the retarded Green's function and G W (t) is known as the Wightman function.
Crucially, the Wightman function depends on the regularization. In the case of the regularized OTOC, we have G W (t) = G(t + iβ/2), which can be computed numerically from the Schwinger-Dyson equations. To determine the Lyapunov exponent, we then perform the following numerical procedure. First, we solve for the eigenvalues of K R for a given value of λ. This relies on the numerical results for G R (t) and G W (t) and the discretization of time into M steps. Second, we perform a binomial search to find λ corresponding to a maximum eigenvalue of one. Finally, we repeat the procedure with different values for M and extrapolate to estimate λ in the continuous limit. The numerical results for λ are shown in Fig. 1(b) of the main text; we verify that the results agree with the low-temperature limit, λ ≈ 2π/β.

Regularized vs. unregularized exponent
We now consider the Lyapunov exponent in the case of the unregularized OTOC. In principle, one can obtain the exponent by calculating the Wightman function with no imaginary-time separation, i.e. G W (t) = G(t), and repeating the numerical procedure outlined above. However, the numerical analysis is more challenging, as the kernel is no longer Hermitian; in particular, we found that the (now complex) eigenvalues are very sensitive to numerical errors that arise from discretization and the imprecision of the Green's functions.
Nevertheless, a simple proof demonstrates that the exponents are the same in the two cases. We start by defining the kernel ansatz in the regularized case as and in the unregularized case as where G (r) W (t) = G(t) and γ u,r (t 12 ) is the normalizable eigenvector with respect to inner product: The proof is then as follows. By definition, (12) can be written in terms of the G We next reparameterize the time variables as t 1 = t 1 + iβ/4, t 2 = t 2 − iβ/4, t 3 = t 3 + iβ/4, and t 4 = t 4 − iβ/4. Using the fact that the integrand goes to zero at t 3,4 → ±∞, we can deform the integration contour and get: In addition, using the fact that γ r (t) is normalizable, one can show that γ r (t + iβ/2) is also normalizable with respect to the inner product for the unregularized case. We thus recover the regularized equation (11) and identify the relation γ (r) (t) = γ (u) (t + iβ/2). In summary, changing the regularization has no effect on the growth exponent but only on the the eigenfunction f (t) (which controls the magnitude of OTOCs, as shown below). More generally, this argument applies to any degree of regularization, i.e. G W (t) = G(t + iη), and to any other many-body chaotic model whose OTOCs are described by ladder diagrams.

Magnitude of OTOCs
Recently, Gu and Kitaev derived an identity that relates the magnitude of the leading order growth term in OTOCs with other quantities computable from the Schwinger-Dyson equations and the kernel equation. We define 9 FIG. 7. The normalized magnitude, C1, of the leading-order term in the OTOC, C1/N e λt , as a function of temperature. The magnitude is calculated numerically (black) for the regularized OTOC using (17); at high temperatures C1 ≈ 1.4 and at low temperatures C1 ≈ 0.5βJ, in agreement with the semiclassical solution, (24) (blue). For the unregularized OTOC, the semiclassical solution predicts the low-temperature scaling C1 ∼ (βJ) 3 (red).
the normalized magnitude, as in the main text, as The identity is then given by where t B = k (λ) is the "branching time", k is the eigenvalue of the retarded kernel, (10), and (γ, γ) is given by We note that γ(t) and G W (t) are defined in the previous section, and both depend on the choice of regularization. For the regularized OTOC, we solve for C 1 as a function of temperature using (17) and the numerical solution of the kernel equation. These results are shown in Fig. 8. At high temperatures C 1 approaches 1.4, while at low temperatures C 1 ≈ 0.5βJ. This latter result is consistent with previous work and provides validation of our numerical methods.
For the unregularized OTOC, the magnitude can theoretically be calculated following the same approach; however, we confront the same difficulties regarding the kernel equation as in the previous section and thus leave this computation for a future work. Of course, in the limit βJ → 0, the magnitude must be the same regardless of regularization. Moreover, we will show in the following section that the low temperature scaling can be computed from the Schwarzian action, leading to the result C ∼ (βJ) 3 .  (5), while the Schwarzian action (19) captures finite N behavior at low temperatures. In the limit N β 1, both theories approach the conformal limit (6).
where f (τ ) is a function describing reparameterizations of time, i.e. τ → f (τ ). Prior work has established the relation between the coupling coefficient C Sch and parameters of the SYK model: C Sch ≈ 0.01N/J [9] [12]. This prefactor controls the size of fluctuations about the C Sch → ∞ saddle-point solution.
Recently, analytical methods have been developed to solve the full dynamics of the Schwarzian action, enabling the calculation of correlation functions at all orders in 1/C Sch . Furthermore, these developments have established a direct correspondence between the Schwarzian action and near AdS in 1+1 dimensions. In particular, the saddle-point solution is dual to classical gravity, while higher order 1/C Sch corrections correspond to gravitational fluctuations.

Two-point correlators
We consider the two-point function, G(z) = χ i (z)χ i (0) β , where z = it + τ is complexified time and the thermal average is computed at inverse temperature β. The exact result from the Schwarzian action is given by [10,13,14] where ρ(s) = s 2π 2 sinh(2πs) is the density of states, Γ(·) is the Gamma function, and the normalization factor is equal to The behavior of G(z) can be understood qualitatively in several regimes. At short times t C Sch , the integrals are well-approximated by the classical saddle point, leading to the aforementioned conformal solution (6). At late times t C Sch , however, the behavior is dominated by the low-energy edge of the spectrum, which respects the Wigner semicircle law: ρ(E) ∼ √ E or ρ(s) ∼ s 2 . This gives rise to a power law decay with an exponent independent of the operator dimension. In particular, we can identify two cases, depending on the temperature relative to C Sch : (a) High temperature, β C Sch : (b) Low temperature, β C Sch : While observing a clear separation between these regimes is challenging at finite-sizes, we expect our numerical results to be described quantitatively by the full functional form of G Sch (z). To this end, we compute the integrals in (20) numerically in Mathematica. The integrals converge quickly for the imaginary-time correlator, G Sch (τ ); for the real-time correlator, obtaining convergence requires us to introduce a small imaginary-time separation, i.e. G Sch (it) → G Sch (it + ). The numerical results in real and imaginary time are shown in Fig. 9 and 10, respectively, for the temperatures and timescales relevant to our study. We note that the retarded Green's function corresponds to the real part of G Sch . A key feature of this correlator is a non-monotonic decay with respect to temperature; this results from the non-trivial dependence of the phase on t and β in (23) and is in stark contrast with the prediction of the conformal solution or of the semiclassical solution of the SYK model (i.e. the Schwinger-Dyson equations).

Out-of-time-order correlators
Although previous studies have derived an exact integral expression for the OTOC analogous to (20), the integrals are significantly more complex and solving them numerically is beyond the scope of this study. Instead, we rely on the semiclassical limit of the OTOC, which is valid for C Sch → ∞, t 1/λ log C Sch . This expression was derived using the correspondence to quantum gravity (i.e. by summing over tree-level graviton diagrams) and is given by where z i = t i +iτ i are the complexified times at which the four operators are inserted, i.e. V (z 1 )W (z 3 )V (z 2 )W (z 4 ) , and U (·) is the confluent hypergeometric function. From (24), we can determine the behavior of different regularizations. The regularized OTOC corresponds to which leads to Crucially, this expression can be expanded in powers of e λt /C Sch . For example, setting ∆ = 1/4, leads to As C Sch ∼ N , this manifestly satisfies the ansatz presented in Eq. 5 of the main text. We further note that the magnitude of the leading order term agrees with the numerical results determined in the previous section.
To obtain different regularizations, we can reduce the imaginary-time separation between the operators. For example, setting corresponds to a symmetric separation with energy scale η. It is evident that any finite value for η leads to a welldefined OTOC with the same Lyapunov exponent as the fully regularized case. To represent the fully unregularized correlator, the naive expectation is to take the limit η → 0, which causes the denominator in (24) to vanish. For the SYK model, this UV divergence is clearly unphysical, and one should instead impose a microscopic cutoff of order J. The net effect is to enhance to growth term by a factor of (βJ) 2 and thus decrease the scrambling time by a factor of log(βJ) 2 . More precisely, the leading order term for the unregularized correlator is given by where C 1 is an order one prefactor. The exact numerical value of C 1 cannot be determined by these methods, as the microscopic cutoff corresponds to "smearing" the operators over an energy scale J rather than setting an exact value for η. To summarize, regularization changes the magnitude of the exponential growth but has no effect on the Lyapunov exponent. The intuition behind this conclusion can be understood from the dual gravitational theory, where the Lyapunov exponent corresponds to the coupling strength of the graviton interaction and the regularization corresponds to the initial energy of an incoming shock wave. While the energy of the initial state has no effect on the coupling strength, it determines the timescale at which nonlinear graviton effects become relevant, leading to the saturation of the correlator. Specifically, the unregularized correlator corresponds to a higher energy initial state, which reaches saturation at an earlier time. While this intuition applies directly to the SYK model at low temperatures, we expect the same qualitative effects to hold at all temperatures due to the form of the ladder diagrams; in the general case, the graviton interaction would be replaced by a 'scramblon' interaction whose strength is governed by λ(β).