Chaotic-Integrable Transition in the Sachdev-Ye-Kitaev Model

Quantum chaos is one of the distinctive features of the Sachdev-Ye-Kitaev (SYK) model, $N$ Majorana fermions in $0+1$ dimensions with infinite-range two-body interactions, which is attracting a lot of interest as a toy model for holography. Here we show analytically and numerically that a generalized SYK model with an additional one-body infinite-range random interaction, which is a relevant perturbation in the infrared, is still quantum chaotic and retains most of its holographic features for a fixed value of the perturbation and sufficiently high temperature. However a chaotic-integrable transition, characterized by the vanishing of the Lyapunov exponent and spectral correlations given by Poisson statistics, occurs at a temperature that depends on the strength of the perturbation. We speculate about the gravity dual of this transition.

Motivated by its potential applications in high-energy and condensed matter physics, and also because of its simplicity, research on fermionic models with infinite-range random interactions [1][2][3][4][5][6][7][8][9], now generally called Sachdev-Ye-Kitaev (SYK) models [10][11][12][13], has flourished in recent times [11,.Interesting research lines currently being investigated include not only applications in holography [10][11][12][13] but also in random matrix theory [25-27, 30, 32, 34], possible experimental realizations [19,35,36], and extensions involving nonrandom couplings [24,28], higher spatial dimensions [18,21,31,37,38], and several flavors [39].A natural question to ask [18,21,24,31,[37][38][39] is to what extent holographic properties are present in generalized SYK models.For instance, similar features are observed for nonrandom couplings [24] and in higher-dimensional realizations of the SYK [37,38] model.However, in some cases, the addition of more fermionic species can induce a transition to a Fermi liquid phase [31] or a metal-insulator transition [18,21], which, at least superficially, spoils a holographic interpretation.Here we study the stability of chaos and holographic features of a generalized SYK model consisting of N fermions in 0 + 1 dimension with infinite-range two-body random interaction perturbed by a one-body random term where χ i are Majorana fermions so {χ i , χ j } = δ i j .The couplings J i jkl and κ i j are Gaussian-distributed random variables with zero average and standard deviation √ N respectively [10,12].We study the model [12,13] by introducing replica fields, averaging over disorder and decoupling the replica fields by two Hubbard-Stratonovich transformations that allow the integration of the original fermionic variables.The resulting partition function is expressed in terms of the bilocal fields G(τ 1 , τ 2 ) and Σ(τ 1 , τ 2 ): Z = [DG][DΣ]e −NS eff , where The saddle point equations in imaginary time, which become exact in the large-N limit of interest, are: where ω n = (2π/β)(n + 1/2), G n ≡ G(iω n ) and Σ n ≡ Σ(iω n ).In the long time, strong coupling limit, where conformal symmetry holds, the solution of the Schwinger-Dyson (SD) equations is dominated by the one-body term and therefore [10,12] the zero-temperature entropy always vanishes.The low-temperature limit of the specific heat is directly related to the leading correction to the conformal Green's function.A simple power-counting argument in the SD equations suggests that it has contributions from both terms.Therefore we expect the specific heat still to be linear with a slope c that may depend on κ.We confirm these results by exact diagonalization of Eq. (1) with J = 1.For a given set of parameters, we have obtained at least 10 6 eigenvalues.We have computed, following [26,30], the entropy at zero-temperature s 0 and the specific heat by using standard thermodynamic relations and a finite-size scaling analysis.As was expected, we have found a vanishing s 0 for any κ and a linear specific heat with c ∝ N f (κ) with f ∼ 0.5/κ for small κ and a steady increase for larger κ.We note that all these features, including s 0 = 0 [40], are consistent with the existence of a gravity dual.We have confirmed these results by an explicit evaluation of the free energy from Eqs. (2) and (3).See Suppermental Material, Secs.A and B for more details.Next we employ level statistics [41] to investigate the effect of the one-body perturbation in the quantum chaotic features of the model.
Level statistics.-For that purpose, we compute, from the exact diagonalization of Eq. ( 1), the level spacing distribution P(s), the probability to find two consecutive eigenvalues  4) as a function of κ for different N's.For sufficiently large N we observe a crossing at κ c ≈ 66 which suggests the existence of a chaotic-integrable transition.Results for larger N would be necessary to confirm it.See main text for an explanation of the absence of crossing for small N.Both P(s) and r were computed by using ensemble and spectral average in a window comprising 10% of the eigenvalues around the center of the spectrum.Results are robust to changes in the percentage of eigenvalues provided that the spectrum edges are avoided.
that probes the system dynamics for times of the order of the Heisenberg time ∼ /∆ with ∆ the mean level spacing.For an insulator, or a generic integrable system, it is given by Poisson statistics, P P (s) = e −s [41], while for a quantum chaotic system it is given by WD statistics [44] which is well approximated by the Wigner surmise, P W (s) ≈ 32 π 2 s 2 exp −4s 2 /π , for systems with broken time reversal invariance [41].For a meaningful comparison with these predictions, unfolding the spectrum [41] is necessary so that ∆ = 1 by a local fitting of the numerical spectral density by a smooth function which is subsequently employed to rescale the spectrum.
Results for P(s) as a function of κ, depicted in Fig. 1, show a gradual crossover from Poisson to WD statistics as κ decreases, which is also observed in the tail of the distribution (see inset).Unlike the standard SYK model [25,26], the one-  5) for N = 30, J = κ = 1 and different β's.For small β, we observe the correlation hole [42,43] followed by a ramp typical of quantum chaotic systems.As β increases, that probes the tail of the spectrum, results are not conclusive.Lower: A finite-size scaling analysis, also with J = 1, of the average adjacent gap ratio r β Eq. ( 4), where, unlike the previous figure, the average is weighted by the function exp(−β(E i + 2E i+1 + E i+2 )/4) but the spectrum is not unfolded.Here we have excluded the ten smallest eigenvalues from the analysis because r for such eigenvalues at the spectral tail is anomalous high even in the large κ limit.For details, see the supplemental material.For β = 0.2, which probes the low energy part of the spectrum, we observe a crossing at κ c ≈ 25 that seems to indicate a chaotic-integrable transition.However the size dependence is too weak to confirm the existence of the transition.body term breaks time reversal invariance for all N.As was mentioned previously, for J = 0 the Hamiltonian is effectively noninteracting which suggests that Poisson statistics applies in the N → ∞ limit.In order to determine whether Poisson statistics is robust for J κ, and therefore a transition occurs at a finite κ, we carry out a finite-size scaling analysis employing as scaling variable the adjacent gap ratio [45][46][47], for an ordered spectrum The average adjacent gap ratio for a Poisson distribution is r P = 2 ln(2) − 1 ≈ 0.386 while for WD statistics is ≈ 0.599 [48].In Fig. 1, r is depicted as a function of κ for different N's.Only 10% of the total number of eigenvalues, located around the center of the spectrum, are employed in the calculation.We stick to values of N for which time reversal symmetry is broken even in the absence of one-body term in Eq. (1).As was expected, except for κ 10, r is very close to the WD result for any N.Only for κ ≥ 25, a crossover to Poisson statistics is observed.For N ≤ 22, we did not observe a crossing point in the plot, and moreover, r gradually approaches the WD prediction, which suggests that the system is quantum chaotic for any κ.However, for N ≥ 30, we observe a crossing at κ = κ c ≈ 60.This is an indication of a transition from chaos to integrability, though it would be necessary to explore larger N's to confirm it.A possible explanation for the different behavior for small N is that the lowest eigenvalues (the most infrared part of the spectrum) are strongly correlated.The reason is that the one-particle sector for κ → ∞, which controls the lowest energy properties, is known [30] to be described by a skew-orthogonal random matrix.The number of eigenvalues related to one-particle states decreases exponentially with N, which would explain why its contribution is only relevant for sufficiently small N.
We note the existence of the gravity dual is related to the properties of the model in the low-temperature, strong coupling limit described by the tail, not the bulk of the spectrum studied above.Moreover, we are also interested in the nature of level statistics for shorter timescales where random matrix theory predicts level rigidity [41].In order to investigate these issues, we compute the connected spectral form factor of the unfolded spectrum where Z(t, β) = Tre −βH−iHt and β > 0. For quantum chaotic systems, and also for the unperturbed SYK model [26,30], we expect a correlation hole [30,42,43,49,50] for intermediate times followed by a ramp, related to the level rigidity observed in quantum chaotic systems [41].We note that, by increasing β, we probe the tail of the spectrum.Results, depicted in Fig. 2, show that for κ = 1 the ramp is still observed for sufficiently small β.In order to clarify the situation for larger β, which probes the spectral correlations of the smallest eigenvalues, we again carry out a finite-size scaling analysis of the averaged adjacent gap ratio r [Eq.( 4)], but the average is weighted by β (see caption of Fig. 2) so that the low energy part of the spectrum is singled out.A crossing seems to be observed (see lower plot of Fig. 2) at κ = κ c ≈ 25.This is a signature of a chaos-integrable transition in the tail of the spectrum.However, results are not conclusive because the size dependence is weak for large κ.We only note that, in qualitative agreement with the results of next section, κ c for β = 0.2 is smaller than in Fig. 1, where effectively β ≈ 0. It is worth mentioning that similar chaoticintegrable transitions in level statistics have been previously studied in the context of nuclear physics [51] and quantum chaos [7,[52][53][54] in somehow related models such as complex fermions with infinite-range interactions and a random diagonal one-body term or interacting systems with short-range interactions [55].
In summary, the finite-size scaling analysis is not fully conclusive to detect the chaotic-integrable transition.In order to confirm it, we investigate next out-of-time-order fourpoint correlation functions where quantum chaotic features are characterized by a finite Lyapunov exponent [56][57][58].Out-of-time-order four-point correlation function.-Inthe semiclassical limit, the time evolution of certain out-of-timeorder correlation function experiences a period of exponential growth [58,59] around the Ehrenfest time t * ∼ λ −1 L log( /S 0 ), where S 0 is a typical action of the system and λ L is the classical Lyapunov exponent.By contrast, for nonchaotic systems, the growth of t * with is only power law [60].The application of these ideas in high-energy physics, where is traded by a parameter ∼ 1/N that controls small quantum gravity corrections, has led to the proposal that black holes are quantum chaotic [56] with a Lyapunov exponent that saturates a recently proposed universal bound λ L ≤ 2πk B T/ [57].We now study whether these chaotic features are present in Eq. ( 1).We compute λ L from the following [12,57,61] out-oftime-order correlator Tr ρ(β) where is inserted [57] along the thermal cycle to regularize the otherwise divergent operator.It is possible to show that F (t 1 , t 2 ) satisfies where G R (ω) = G(iω n → ω+i0 + ) is the retarded Green's function in real frequency and G lr (t) is the Wightman function obtained from G lr (ω) = 2ie −βω/2 1+e −βω Im(G R (ω)).In order to compute these two-point functions we follow the strategy employed in Refs.[12,31].We analytically continue iω n → ω + i0 + the saddle point equations (3) and solve them using the spectral representation of the retarded Green's function.Substituting the ansatz F (t 1 , t 2 ) = e λ L (t 1 +t 2 )/2 f (t 12 ), where t 12 = t 1 − t 2 , into Eq.( 7) and expressing it in the frequency domain, we obtain the following eigenvalue equation for f (ω), where ω = ω 1 − iλ L /2 and g lr (ω) = dte iωt G lr (t) 2 .Finally, we compute λ L by imposing the existence of a nondegenerate eigenvalue equal to one so that Eq. ( 9) is satisfied.Results, depicted in Fig. 3, show the system displays chaotic behavior, namely, the Lyapunov exponent λ L is finite, for all studied values of κ and sufficiently high temperature.However, even in the strong coupling limit, λ L never FIG. 3. Lyapunov exponent λ L for the model Eq. ( 1) with J = 1 as a function of the inverse temperature β = 1/T and κ.From top to bottom, κ = 0, 0.2, 0.5, 1, 2. A finite λ L , which is a signature of quantum chaos, is observed for a fixed κ and not too low temperature.For sufficiently low temperatures, and κ > 0, we identify a T * (κ) such that λ L = 0 for any T < T * (κ) which signals a chaotic-integrable transition.(Inset) Critical value of κ = κ c at which the transition takes place as a function of β.Dots result from fitting the numerical data of the main plot near the transition.The solid line is the analytical expression from Eq.( 15) valid in the large-q limit.Agreement with the numerical data is reasonable for κ/J 1.
approaches the bound λ L = 2πk B T/ .Indeed, for a given temperature, λ L decreases as κ increases and eventually vanishes for sufficiently strong κ or, for a fixed κ, for sufficiently low temperature.Therefore, quantum chaos is robust to the introduction of a relevant one-body perturbation but only if it is weak enough and the temperature is high enough.We now confirm these results analytically by studying the following model with q/2-body interactions, where κ i j and J i 1 ,i 2 ,...,i q are again Gaussian-distributed random variables with zero average and q-dependent variances κ 2 qN , 2 q−1 q (q−1)!J 2 N q−1 , and q 1, respectively.As before, we fix J and use κ and β as the only parameters.The key insight is that the retarded kernel in this model, given by simplifies considerably in the limit of q 1, allowing for an analytical solution of the eigenvalue problem in Eq. (7).First, we proceed in the same way as for the model (1) to find an effective action and the associated saddle point equations analogous to Eqs.( 2) and (3).In the limit q 1, the saddle point equations can be consistently expanded in terms of yielding a non-linear boundary value problem for g, with θ = τ/β ∈ (0, 1).The retarded kernel Eq.( 11) is thus obtained by analytical continuation of the resulting G(τ) and, as mentioned before, for q 1 it is given by a simpler expression, where g(τ) is the solution of Eq.( 13) and is given explicitly as a power series in κ/J 1 in the Supplemental Material, Sec. C. We again use the ansatz F (t 1 , t 2 ) = e λ L (t 1 +t 2 )/2 f (t 12 ) in order to rewrite the eigenvalue problem in Eq.( 7) as a Schrödinger equation for f (t 12 ).The eigenstates of the resulting equation are found perturbatively in κ/J 1, giving a correction O(κ 2 ) to the Lyapunov exponent.This correction is given in terms of an integral that, for low temperature βJ 1, is approximated by The transition occurs when λ L = 0, which leads to a βdependent critical κ = κ c .For instance, κ c (β=133) ∼ 0.2J and κ c (β=53) ∼ 0.5J, which is in good agreement (see Fig. 3) with numerical results for q = 4.We refer to Sec.C of the Supplemental Material for additional details.
Finally, we note these types of transitions are generic [55], so it would be interesting to identify their gravity dual.We speculate with the possibility that the gravity dual of the transition studied in this Letter is a Hawking-Page transition where the black hole and thermal gas phases correspond to the chaotic and integrable phase, respectively.In conclusion, we have found that the SYK model perturbed by a random one-body term is still chaotic in the limit of sufficiently high temperature or weak perturbation.However, for a given strength of the perturbation, the system undergoes a chaotic-to-integrable transition for sufficiently low temperatures which may have a gravity dual interpretation.
Note added: Close to completion of this work, we became aware of three papers [62][63][64] that study a somehow similar generalized SYK model though the focus of these papers is rather different.Ref. [64] studies quantum quenches while the other two investigate a two fermion species generalization of the model in which a transition occurs by tuning the number of fermions.
We thank D. Anninos, S. Banerjee, and P. Sabella-Garnier for illuminating discussions.Part of the computation in this Letter has been done using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo. A. M. G. acknowledges partial financial support from a QuantEmX grant from ICAM and the Gordon and Betty Moore Foundation through Grant GBMF5305.The work of M. T. was partially supported by Grants-in-Aid No.In this appendix we present explicit results for the low temperature limit of the entropy and the specific heat that, for space limitations, could not be included in the main text.The zero temperature limit of the entropy, s 0 , was obtained by exact diagonalisation of the Hamiltonian Eq.( 1) for different N and κ.We note that for any finite N the entropy will vanish in the T → 0 limit [65].Therefore, to justify a zero entropy in this limit it is necessary to extrapolate the finite N numerical results to the N → ∞ limit.However, this does not seem strictly necessary as the N dependence, depicted in Fig. 4, is very weak, especially for κ ≥ 1.A simple extrapolation of the curves for larger N to the T = 0 limit leads to a zero-temperature entropy which is smaller than 10 −3 for all κ's.These results are fully consistent with the N → ∞ results, which have been obtained by fitting In order to obtain log Z N we solve Eq. ( 3) as described previously in Refs.[12,38].We use a Fast Fourier transform to switch between frequency and time domains and solve iteratively until convergence.For ω n = 2πT (n+1/2) we take −N ω /2 < n < N ω /2−1, with N ω = 2 24 and N t = 4N ω points in the frequency and time domains.
We now move to the study of the low temperature limit of the specific heat C(T ) per Majorana obtained by exact diagonalisation.For that purpose we employ the following thermodynamic expression: where Z is the partition function, . . .stands for ensemble average and k labels the eigenvalues for a given disorder realisation with average Ē.
We carry out quenched averages, namely, the specific heat is computed separately for each disorder realisation.The final specific heat is the arithmetic average over all disorder realisations.We then fit the low temperature limit by a low order polynomial in temperature.The coefficient of the linear term is the specific heat coefficient, c, which in units of N and with the coupling constant set to one, is π/(6κ) for κ → ∞ and ≈ 0.4 for the unperturbed two-body SYK model (κ = 0) [12,26].As shown in Fig. 5, for large κ, c is indeed very close to π/(6κ).Only for κ 1 we observe a moderate increase of c.  16) and the exact diagonalisation of Eq.( 1) for N = 34 and different κ's.The specific heat is clearly linear in the low temperature limit with a slope that it is close to the κ = 0 prediction c = π/6 for any κ ≤ 1.This is a further confirmation that the ground state and lowest energy excitations of the Hamiltonian Eq.( 1) are mostly controlled by the random-mass term.We cannot study arbitrarily small κ because it would require to reach very low temperatures that compromise the numerical accuracy of the results.We have checked that the obtained specific heat coefficient c is robust to changes in the fitting interval.We have observed that for N ≤ 30 the value of c does not have a monotonic dependence on N which makes difficult to carry out a finite-size scaling analysis.For that reason, unless otherwise stated, the fitting is restricted to the largest N = 34 that can be reached numerically.Finally, in Fig. 6, we compare the specific heat coefficient c obtained from the large-N fitting of log Z N , as explained above with the exact diagonalisation result given in Fig. 5. Deviations are consistent with 1/N corrections only taken into account in the latter.FIG. 6. Specific heat coefficient c as a function of κ.The exact diagonalisation result, shown in blue crosses, is the slope of the fitting curve in Fig. 5.The black dots are obtained from fitting −β F N = log Z N obtained from the numerical solution of the large-N saddle point equations Eq.( 3).The dashed line corresponds to c = π/(6κ), the analytical value when the Hamiltonian only contains the random-mass term [12].Differences between the two results are consistent with 1/N ∼ 0.03 corrections only retained in the exact diagonalisation result.

Appendix B: Analytical Low Temperature Thermodynamic properties
We carry out an analytical calculation of the low temperature thermodynamic properties of the generalised one+twobody SYK model in the limit of large N. Since the exact solution for the one-body Hamiltonian is known analytically, it will be convenient to develop the perturbative approach around this exact ground state.For this purpose, it will be convenient to consider a rescaled model given by where the standard deviation of the random couplings J i jkl and κ i j are given by √ 6κ N 3/2 and 1

√
N respectively [66].The disorder averaged effective action associated with the rescaled Hamiltonian is given by where is the Euclidean many-body correlator.
We will proceed in two steps, first we find a perturbative solution of the rescaled SD equations and then we compute the low temperature limit of the free energy by substitution of the obtained Green's function in the action.
Perturbative solution of the SD equation.-TheSD equations derived from Eq.( 18) can be written in Fourier space as where we have defined σ(ω) = dτ e iω n τ G β (τ) 3 .For κ = 0, the system is noninteracting and Eq.( 19) simplifies to an algebraic quadratic equation in G, which is solved exactly by The subindex '0' in G 0 (ω n ) indicates it is the solution to the κ = 0 noninteracting limit.For κ > 0, Eq.( 19) is an integral equation, and cannot be solved exactly.We proceed with a perturbative solution in the limit κ 1.Let Inserting in Eq.( 19) and expanding in κ, we find where σ 0 (ω) = dτ e iωτ G 0 (τ) 3 acts as a source for the first order correction.The prefactor can be written exactly, However σ 0 (ω n ) requires more effort since we first need G β 0 (τ).Note that the only temperature dependence of G 0 (ω n ) is through the Matsubara frequencies, and thus disappear upon analytic continuation iω n → + i0 + .This is expected since κ = 0 is essentially a noninteracting system, and the spectral density should not be temperature dependent.However, for any κ > 0 the system is interacting, and we do expect non-trivial temperature dependence in G(ω n ).As we will see below, this comes exactly from the non-linearity induced by the σ 0 contribution.
Finite temperature.-Weanalytically continue iω n → +i0 + Eq.( 20) to get the zeroth order (in κ) retarded Green's function which corresponds to the noninteracting limit Note that, as discussed above, temperature dependence has disappeared.The zeroth order spectral function is where I A (x) is the characteristic function of the set A. Note this is the well-known Wigner semicircle distribution for the spectral density of a random matrix.It should come at no surprise since for κ = 0 the Hamiltonian of the system is a sparse random matrix.The spectral function allow us to analytically continue G 0 (ω n ) to the whole complex plane as and provide a useful way of getting the different Green's functions.To get the finite temperature Matsubara Green's function, we set z = iω n and sum over the frequencies G β (τ) = n∈Z e iω n τ G(iω n ), as usual in the Matsubara formalism: In principle we have everything we need to compute σ 0 (ω n ) and thus the first order correction g(ω n ).However, in the lack of a closed expression, we proceed with a low temperature expansion of Eq.( 27).More, precisely we expand the integrand of Eq.( 27), where ρ(λ) = ρ 0 (λ) is given in Eq.( 25), for large β.We then perform the integration order by order giving where we define θ = τ/β ∈ (0, 1).We now expand G β 0 (τ) 3 in 1/β and perform the Fourier transform to get σ 0 (ω n ) = β 0 dτ e iω n τ G β 0 (τ) 3 .Inserting this into Eq.(22), together with Eq.( 23), we get an expression for g(ω n ) that should be added to the zeroth order solution Eq.( 20) (see Eq.( 21)).Analytically continuing this expression iω n → + i0 + and expanding in low frequencies gives an approximation for the analytic continuation of Eq.( 21) where we have conveniently separated the real and imaginary parts.The real part of the retarded Green's function is odd, while the imaginary part which gives the spectral density is even.Note that the effect of interactions is mainly to renormalise the coefficients of the free κ = 0 model, with both zero and finite temperature contributions.This reflects the fact that the interactions are irrelevant compared to the hopping term.At zero temperature β = ∞, corrections only appear at order 2 .This means that the low frequency behaviour of the ground state of the system is unchanged.
Low temperature expansion of the free energy.-Tocompute the free energy density of the model, we have to insert the finite temperature saddle point solution we obtained in the sections above into the effective action Eq.(18).As mentioned earlier, we expect the low temperature expansion to be given by where E 0 is the ground state energy density, s 0 the zero temperature entropy density and c the specific heat coefficient.Our aim in this section is to compute these coefficients for small κ.One can check that the integral terms in Eq.( 2) do not contribute to these lowest order coefficients, and the leading contributions come from the Tr log term.By definition, it is given by We can write this sum over Matsubara frequencies as the residue of the complex integral γ dz 2πi log G(z)n(z)e z0 + where n(z) = (1 + e βz ) −1 and the contour γ englobes the poles at z = iω n .The integrand also has a branch cut along the real axis.Since the integral decay at infinity, we can deform γ to wrap around the branch cut in the real axis.This leads to where G R ( ) = G(iω n → + i0 + ).For κ = 0, G R ( ) is simply given by Eq.( 24).Thus, giving where in the last equality we used the Sommerfeld expansion for the Fermi-Dirac distribution n( ) = ).The log term is exponentially decaying, and does not contribute in Eq. (30).For κ = 0 we restore the units by taking the coupling κ i j to have a standard deviation parametrized by a scale K: K/ √ N. Therefore, we conclude that for κ = 0 we have E 0 = (1 − π) KN, s 0 = 0 and c = π 6 N K .These values are consistent with the previous results in the literature for the SYK model at q = 2 [12].Note that the term proportional to β −1 in the low temperature expansion only depends on derivatives of the argument evaluated at = 0. Thus for κ = 0 the specific heat coefficient only depends on the low frequency properties of the retarded Green's function.
For κ > 0, we only have access to the IR low frequency result in Eq.( 29).Naively inserting this in Eq.( 32) leads to UV divergences due to unboundedness of the integrals.Note this was not needed for κ = 0 since the compact support of the spectral function comes naturally from the exact solution.But UV divergences are expected since the perturbative IR solution does not hold for large frequencies.We therefore introduce an UV cutoff λ to regularise the integrals.We will show that this leads to cutoff dependent ground state energy, but cutoff independent entropy and specific heat.This is expected given that, as discussed above, the specific heat is determined by the IR behaviour of G R .Although there is some freedom in the choice of λ, it is largely constrained by non-perturbative properties of the spectral function ρ( ) = 2Im G R , which must satisfy the sum rule d 2π ρ( ) = 1 and positivity ρ( ) ≥ 0. These constraints give us a consistent way to fix the cutoff: we choose λ such that λ −λ d 2π ρ( ) = 1 [67].This ensures that | | < λ saturates the spectral weight, and therefore for | | ≥ λ the perturbative solution breaks down.It is easy to check that this choice automatically satisfies positivity of ρ.With this discussion in mind, we take where we define a( , κ, β) = tan −1 Re G R Im G R with G R given in Eq. (29).As in Eq.( 34), the integral over (λ, ∞) gives an exponentially decaying term that does not contribute to the thermodynamic coefficients.The remaining integral over (−λ, λ) can be integrated numerically using the cutoff estimated from the sum-rule for ρ. [68].To extract the specific heat coefficient, we subtract the zero temperature result and multiply by a factor β 2 .According to Eq.( 30), the result should asymptote to c/2 as β → ∞.These are shown in Fig. 7 and Fig. 8 respectively.FIG. 7. Numerical integration of β 2 (F/N − E 0 ), which asymptotes to c/2 at low temperatures.
One can observe in Fig. 7 an order 10 −3 correction in c/2 for small κ.While this correction is consistent with the exact diagonalisation results in Appendix A, it can also be an artefact of perturbation theory.We thus study the dependence of c/2 in κ as we go higher orders in perturbation theory for G R ( ).This is shown in Fig. 9.
The specific heat coefficient c/2 increases with κ.However, the higher order we go in perturbation theory for small frequencies , the smaller is the increase for a given κ.This suggests that a fully non-perturbative calculation of c(κ) could lead to a κ independent specific heat coefficient at least in the limit κ 1.
Indeed this can also be understood analytically.First we need to identify which term gives this contribution.The integral over the π/2 factor in Eq. (35) gives only a cutoff depen- From Eq.( 29), we can check that d d a| =0 = − 1 2 for any κ and β.Thus this term gives identical to that for κ = 0. Since a is now temperature dependent, possible corrections to the specific heat coefficient can come from the integral of a over (−λ, 0) or from the higher order odd derivatives.However a simple series expansion of a reveals that the leading order temperature dependence in a comes at order O( 9 ), which is beyond the scope of the perturbative result in Eq. (29).One can check that, had we gone only to order O( 4 ) in G R , the first order temperature dependence would have been at order O( 5 ).Going to order O( 6 ) pushes the leading order temperature dependence of a to order O( 7), and finally going to order O( 8 ) pushes it to O( 9 ).This analytical argument is fully consistent with the evaluation of c from Eq.( 32) depicted in Fig. 9.
In Appendix A, the low temperature thermodynamic coefficients were studied numerically by exact diagonalisation.Corrections to the κ = 0 specific heat coefficient c = π/6 and entropy density s 0 = 0 for κ 1 were found to be of order O(10 −3 ) or lower.The analytic calculations from this Appendix confirm these results.They also corroborate the claim that the ground state of our model is dominated by the κ = 0 limit of Eq.( 1).
Appendix C: Analytical calculation of the Lyapunov exponent with q/2-body interaction and q 1 In this appendix we aim to give analytical support to the numerical results depicted in Fig. 3 where it was shown that, for a fixed value of the coupling constant κ > 0, the Lyapunov exponent λ L always vanishes for sufficiently low temperatures.We have managed to obtain analytical results not for the Hamiltonian Eq.( 1) but for a closely related model where the two-body term is replaced by a q/2-body, with q 1, interaction keeping the one-body random perturbation of Eq.( 1).The motivation for this choice is that, without the one-body perturbation, the Lyapunov exponent can be computed analytically at any coupling in the large q limit, providing an explicit setup for studying the saturation of the chaos bound at low temperatures [12].
By following closely the method of Ref. [12], we show below that the perturbative expansion in κ/J 1 captures the relevant physics.For any fixed κ/J 1, we identify a range of temperatures for which the Lyapunov exponent is non-zero, though never saturates the bound on chaos.We also show that it vanishes for sufficiently low temperatures.This is in full agreement with the numerical results for the twobody model which corroborates the existence of a chaotic-tointegrable transition in this type of generalised SYK models.
The model and the perturbative solution.-Asdiscussed in the main body of the Letter, we work with following Hamiltonian which is the q/2-interaction generalization of the original model Eq.(17).As before J i 1 i 2 ...i q and κ i j are Gaussian distributed random variables with zero average and variances 2 q−1 q (q−1)!J 2 N q−1 and κ 2 qN respectively.Note that, following [12], we have rescaled the coupling constants: κ 2 → κ 2 /q and J 2 → J 2 2 q−1 /q.The large q limit is then defined by taking q 1 and keeping the rescaled couplings fixed [12].We follow the same replica procedure described in the main body of the Letter to get the following effective action As in [12], in the limit q 1 we can consistently expand G as Inserting the above in the saddle point Eq.( 19) and expanding in q, in Euclidean time we can simplify it to where θ = τ/β ∈ [0, 1).Together with the finite temperature boundary conditions g(0) = g(1) = 0, this equation defines a non-linear boundary value problem for g.For κ = 0, the solution is given by Note that for ν = 0 we have βJ = 0 while for ν = 1, βJ = ∞.Thus ν ∈ [0, 1] parametrises the flow of βJ.We linearise Eq.( 40) around g (0) (the κ = 0 solution).More explicitly, we substitute g(θ) = g (0) (θ) + κ J 2 g (1) (θ) + O((βκ) 4 ) into Eq.(39) and then into Eq.(40) to get the equation satisfied by g (1) : where we changed the θ-coordinate to and the corresponding boundary conditions are g πν 2 = g (1) − πν 2 = 0.The solution of this boundary-value problem is given by where we have defined, Lyapunov Exponent.-Asdiscussed in the main manuscript, the out-of-time -order four point correlator is generated by the convolution with the real-time retarded kernel Eq. (7).At large q, this simplifies considerably, K R (t 1 , t 2 , t 3 , t 4 ) = θ(t 13 )θ(t 24 ) 2J 2 e g(τ=it 34 +β/2) + q −1 κ 2 .(46) Thus, at large q the second term is sub leading.By noting that ∂ t θ(t) = δ(t), we can convert Eq. (7) The equation above has the form of a one-dimensional Schrödinger equation for the eigenfunction f with eigenvalues  iy) .For κ = 0, V (0) (y) = 2 sech 2 y is the well studied Pöschl-Teller potential.In particular, it has a bound state with energy E (0) λ = −1 and normalised eigenstate f (0) (y) = 1 √ 2 sech y.This implies that for this eigenstate we have λ L = 2π β ν and therefore since F(t 1 , t 2 ) = e λ L (t 1 +t 2 )/2 f (t 1 −t 2 ) this bound state correspond to an exponential growth for the out-of -time-order four-point function.
We are interested in studying what happens to this bound state when we add the κ/J 1 correction to the potential.Or in a more suggestive notation, when we consider V(y) = V (0) (y)+ κ J 2 V (1) (y) where V (1) (y) = e g (0) (x→iy) g (1) (x → iy) with g (1) (x) given by Eq. (43).By standard quantum mechanical perturbation theory, the correction E λ = E (0) λ + κ J 2 E (1)  λ is given by E (1)  λ = f (0) | 2 cosh 2 (y) g where 2δE(ν) ≡ B(ν) + 19  18 − log 2. Therefore the correction to the bound-state energy is given by Letting we obtain the correction to the Lyapunov exponent, In Fig. 10 we plot Lyapunov exponent for different values of κ with J = 1.This is to be compared with the exact numerical the main text we removed 10 eigenvalues from the lowest end of the spectrum for each sample.
Similar dependence of r i on the eigenstate index i is also observed for the other values of N. In order to extract the features of the model that should survive for large N, we have removed ten lowest eigenvalues in obtaining r β for the lower part of Fig. 2 in the main Letter.The sample average of the gap ratio r i = min(δ i , δ i+1 )/ max(δ i , δ i+1 ), in which δ i = E i − E i−1 is the difference between the neighboring energy eigenvalues in the ordered spectrum, plotted against i.

FIG. 1 .
FIG.1.Upper: P(s) for N = 34 and different κ's with κ in units of J = 1.We clearly observe a crossover from Wigner-Dyson (WD) to Poisson statistics as κ increases.Lower: Finite-size scaling analysis of the averaged adjacent gap ratio r Eq. (4) as a function of κ for different N's.For sufficiently large N we observe a crossing at κ c ≈ 66 which suggests the existence of a chaotic-integrable transition.Results for larger N would be necessary to confirm it.See main text for an explanation of the absence of crossing for small N.Both P(s) and r were computed by using ensemble and spectral average in a window comprising 10% of the eigenvalues around the center of the spectrum.Results are robust to changes in the percentage of eigenvalues provided that the spectrum edges are avoided.

FIG. 2 .
FIG.2.Upper: Connected spectral form factor g c (t) for the unfolded spectrum from Eq. (5) for N = 30, J = κ = 1 and different β's.For small β, we observe the correlation hole[42,43] followed by a ramp typical of quantum chaotic systems.As β increases, that probes the tail of the spectrum, results are not conclusive.Lower: A finite-size scaling analysis, also with J = 1, of the average adjacent gap ratio r β Eq. (4), where, unlike the previous figure, the average is weighted by the function exp(−β(E i + 2E i+1 + E i+2 )/4) but the spectrum is not unfolded.Here we have excluded the ten smallest eigenvalues from the analysis because r for such eigenvalues at the spectral tail is anomalous high even in the large κ limit.For details, see the supplemental material.For β = 0.2, which probes the low energy part of the spectrum, we observe a crossing at κ c ≈ 25 that seems to indicate a chaotic-integrable transition.However the size dependence is too weak to confirm the existence of the transition.
JP26870284 and No. JP17K17822 from JSPS of Japan. A. R. B. is funded through a research programme of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO).B. L. is supported by a CAPES/COT grant No. 11469/13-17.Appendix A: Numerical Thermodynamic properties

18 FIG. 4 .FIG. 5 .
FIG.4.Entropy S as a function of temperature T from the exact diagonalisation of Eq.(1) for different N's and κ where k B stands for the Boltzmann constant.The number of eigenvalues employed is mentioned in the main text.The weak dependence on N is consistent with a vanishing zero temperature entropy.

FIG. 10 .
FIG. 10.Lyapunov exponent λ L , Eq.(51), as a function of the inverse temperature βJ, in units of 1/k B , for different values of κ and J = 1.Note that βλ L 2π ∼ |β − β * | vanishes linearly near the transition, with a slope of approximately −2κ FIG. 11.The sample average of the gap ratio r i = min(δ i , δ i+1 )/ max(δ i , δ i+1 ), in which δ i = E i − E i−1 is the difference between the neighboring energy eigenvalues in the ordered spectrum, plotted against i.