Thouless time for mass-deformed SYK

We study the onset of RMT dynamics in the mass-deformed SYK model (i.e. an SYK model deformed by a quadratic random interaction) in terms of the strength of the quadratic deformation. We use as chaos probes both the connected unfolded Spectral Form Factor (SFF) as well as the Gaussian-filtered SFF, which has been recently introduced in the literature. We show that they detect the chaotic/integrable transition of the mass-deformed SYK model at different values of the mass deformation: the Gaussian-filtered SFF sees the transition for large values of the mass deformation; the connected unfolded SFF sees the transition at small values. The latter is in qualitative agreement with the transition as seen by the OTOCs. We argue that the chaotic/integrable deformation affect the energy levels inhomogeneously: for small values of the mass deformation only the low-lying states are modified while for large values of the mass deformation also the states in the bulk of the spectrum move to the integrable behavior.


Introduction and Summary
The common feature of classical chaotic systems is their high sensitivity to the initial conditions. Concretely, they show an exponential dependence of the dynamical evolution to the initial conditions (for example the initial position) of the form where λ L is the Lyapunov exponent which measures the strength of chaos in the system under investigation. 1 The extension of the concept of classical chaos to the quantum mechanical framework is not straightforward: from equation (1.1) it is clear that the classical notion of chaos relies on the notion of trajectory and so its extension to the quantum world is problematic. Traditionally, quantum chaos has been associated to the well-known Bohigas-Giannoni-Schmit (BGS) conjecture [1], which can be stated as follows: spectra of quantum systems, whose classical analogues are chaotic, show the same fluctuation properties as predicted by random matrix theory (RMT). From the BGS conjecture we see that the quantum chaotic properties of a system are encoded in the statistical properties of its spectrum, and so they can be analyzed using statistical methods. In particular, quantum chaotic systems should show the typical features of RMT, i.e. level repulsion and spectral rigidity. 1 The BGS conjecture have received in the years a large number of numerical and analytical confirmations and it is nowadays accepted as true, even though a proof is still lacking.
On the other hand, another diagnostics of quantum chaos has become very popular in recent days for the study of many-body 2 quantum systems: starting from an old work by Larkin and Ovchinnikov [6] and treating equation (1.1) semi-classically, Kitaev has proposed the out-oftime-ordering correlators (OTOCs) where W and V are 2 generic operators and β is the inverse temperature, as an alternative way to characterize quantum chaos. In particular, if C W V (β, t) displays a behavior of the form where N is the number of degrees of freedom, the system is considered to be chaotic. The definition of quantum chaos based on (1.2) has some very nice features, that have been explored heavily in recent years: it provides a quantum definition of the Lyapunov exponent, which has been later shown [7] to satisfy the universal bound saturated by systems which are dual to Einstein gravity -the so-called maximally chaotic systems. 3 Moreover, (1.2) provides an interesting time scale -the scrambling time -which 1 For an introduction to the BGS conjecture, the field of quantum chaos and RMT techniques, the reader can look at, for example, [2][3][4]. 2 Similar studies have been performed also for single-body quantum systems [5] but the results are not conclusive yet. 3 However some maximally chaotic systems not directly dual to Einstein gravity have been recently discussed [8,9].
can be viewed as the time at which the exponential behavior in (1.3) becomes dominant and quantum effects of chaos appear.
Given the two definitions of chaos we just reviewed, it is natural to ask how (and if) they are related to each other. Of course, finding such a connection would be very important, in order to construct a unique theory of quantum chaos and to find a unique definition of quantum chaotic systems. 4 However, from the way in which they are defined, it is clear that finding a connection between the two definitions is very hard. From a technical point of view, they use completely different tools and languages: the definition of chaos based on the BGS conjecture uses a statistical approach and statistical tools, while the definition of chaos in terms of the OTOCs is formulated in field theoretical (or quantum mechanical) terms. From a conceptual point of view, the two definitions should apply at time scales which are in principle very different from each other: the definition based on the OTOCs is semi-classical and so it is valid at earlier times than the definition based on the BGS conjecture which applies at very late times, when the system has lost completely the memory of the initial dynamics and it can be described by some random hamiltonians which just preserve the symmetry of the original dynamics.
To address this kind of questions, the SYK model [10][11][12][13][14][15] and its generalizations and variants 5 can be useful toy models: Indeed, one of the main attractive features of SYK is that the OTOC can be computed (exactly in the large N and large q limit and numerically in the large N limit with finite q) and shows maximal chaos. Moreover, starting from the works [55][56][57], it has been shown that its spectrum clearly shows RMT features. In particular, [55] has studied the so-called Spectral Form Factor (SFF) for the SYK model, which is nothing but the analytically-continued thermal partition function (we denote with β the inverse temperature) g(t, β) ∝ | Tr(e −(β+it)H SY K )| 2 , (1.5) and has shown that for large enough values of the time variable t it shows the same behavior of RMT. One of the appealing feature of the SFF is that it provides a quantum mechanical flavor to the usual RMT quantities (like the level repulsion and the spectral rigidity) and so it could be in principle closer to the language used in the OTOCs description.
Indeed, some studies addressing the connection between the RMT properties of SYK and the OTOCs appeared in recent days: very recently, [58] has studied the onset of RMT behavior in the (Gaussian filtered) SYK SFF 6 and found that, in agreement with the results of [59], the time-scale at which the RMT behavior appears -which we will call Thouless time -is not the 4 Indeed, a priori there could be systems which are chaotic according to one definition and not chaotic according to the other. 5 For example, see the U (1) symmetry and flavor generalizations [16][17][18][19][20][21], supersymmetric generalizations [22][23][24][25][26][27], SYK(-like) models without disorder [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43], gravity duals [44][45][46][47][48][49][50][51] and the Schwarzian theory [52][53][54]. 6 As we will review, the Gaussian filtered SFF is simply the SFF in which the contribution coming from the 3 scrambling time. Moreover, in [60] it has been argued that, while the RMT dynamics is able to capture the relevant physics at very late times, it fails to reproduce the early time features of chaos and that, for example, RMT would scramble almost immediately -an unwanted feature which would violate the chaos bound (1.4). All together, these results suggest that the connection between the early-time definition of chaos, based on the OTOCs diagnostics, and the late-time definition of chaos, based on RMT analysis, could be very hard (or impossible) to find -even in a very particular model like SYK.
In this paper, motivated by the negative results just mentioned, we will try to test a weaker possibility: even though it is by now clear that the scrambling time and the Thouless time are not the same time-scale -and even though it is clear that the RMT description is unable to capture aspects of the early-time manifestations of chaos -it is still possible that the Thouless time and the scrambling time are functionally related, in a perhaps highly non-trivial way, and that systems which are more chaotic according to the OTOCs criterion develop the RMT behavior at earlier times -i.e. their Thouless time is smaller than the Thouless time of systems which have smaller Lyapunov exponents.
To this purpose, we will focus our attention on a specific variant of the SYK model -which we will call "mass-deformed" SYK model -introduced and studied in [61]. It is characterized by an hamiltonian that, beyond the quartic interaction term of the standard SYK model, contains an additional quadratic piece, a random mass term (1.6) As we will review following [61], this model shows a transition from a chaotic behavior, when the constant κ is small enough and so the dynamics is dominated by the quartic chaotic piece, to an integrable behavior, when the constant κ is large and the dynamics is dominated by the quadratic, integrable, term.
The transition from chaotic to integrable behavior has been observed, in [61], both by a RMT analysis and by studying numerically (and analytically in the usual large q limit) the large N behavior of the OTOCs. However, the two diagnostics of chaos do not see the transition in the same way: the RMT analysis shows the transition from chaos to integrable behavior for values of κ ∼ 15 or larger, while the OTOCs see the transition for values of κ which are much lower, around κ ∼ 1. We will argue that this discrepancy is not just due to the fact that the OTOCs are computed in the large N limit and the RMT analysis is performed at finite N . Instead we will provide arguments that this discrepancy is a consequence of the way in which states in the bulk of the spectrum are exponentially magnified. In this way, it has been shown in [58] that the onset of RMT behavior can be detected to earlier times than in the standard SFF. the chaotic/integrable transition affects the spectrum of the model: indeed, we will see that the integrable transition starts from the low-lying states, i.e. the states closer to the ground state and only after, i.e. for larger values of κ, it affects the states in the bulk of the spectrum. In agreement with this pattern, we will conclude that the OTOCs are mostly affected by the tail of the spectrum, i.e. by the states which are closer to the ground state, while they are pretty insensitive to the states in the bulk of the spectrum, and for this reason they see the chaotic/integrable transition for lower values of κ. On the other hand, the RMT analysis is mostly affected by the states residing in the bulk of the spectrum and so it sees the transition at much higher values of κ when also the states in the bulk show the transition to the integrable behavior. We will provide confirmations of this picture both by means of RMT observables depending only on the energy spectrum (level spacing distribution andr statistics) and by RMT observables depending also on the eigenvectors, like the Inverse Participation Ratios (IPR). All these results point towards a clear separation of the chaotic/integrable transition between the low-lying modes and the states in the bulk, as can be seen in Figure 7b which shows how the IPR displays a different behavior between the low-lying states and the states in the bulk of the spectrum.
With the aim of finding a connection between the behavior of the OTOCs and the RMT analysis, and to find a quantity which can provide an estimation of the Thouless time in agreement with the scrambling time, we will look for RMT observables which could be able to capture the same physics of the OTOCs and that, in particular, are able to see the chaotic/integrable transition for lower values of κ. To this end, we will test both the Gaussian filtered SFF, as introduced by [58], and the connected unfolded SFF, which has been introduced in [62] and that is simply the SFF computed using the unfolded spectrum and removing the disconnected piece. The Gaussian filtered SFF detects the chaotic/integrable transition at large values of κ, in agreement with the RMT analysis of [61]. This result is not surprising: by construction the Gaussian filtered SFF magnifies the bulk of the spectrum with respect to the tail and hence it sees the transition at large values of κ. 7 On the other hand, we will see that the Thouless time computed using the connected unfolded SFF, especially for low enough values of the temperature, shows a pattern which is qualitatively similar to the pattern seen by the OTOCs. In particular, we will see that the chaotic/integrable transition arises for low values of κ and that, at κ ∼ 1 the chaotic features mostly disappear. We will then use the connected unfolded SFF to estimate the Thouless time for different values of κ and N and we will see that the chaotic/integrable transition arises around k ∼ 1 and that it becomes sharper and sharper when N increases. These features, and the qualitative agreement between the OTOCs and the connected unfolded SFF, suggests that the connected unfolded SFF could be an important observable to study the relation between the early-time chaos diagnostics and the BGS conjecture, in the mass-deformed SYK model as well as in other models. Of course, our study has several technical and conceptual limitations, that we plan to address elsewhere and that we will discuss in the final section of this paper.
The paper is organized as follow. In Section 2, we will review the main features of the mass-deformed SYK model, as introduced in [61].
In Section 3, we will review the notion of the Thouless time as the time-scale at which the RMT behavior appears. We will mention about the differences between this time-scale and the scrambling time and we will introduce our main diagnostic tools: the connected unfolded SFF, as introduced in [62], and the Gaussian-filtered SFF, introduced very recently in [58].
In Section 4, we will present our numerical results for the computation of the connected unfolded SFF. We will discuss that it seems to capture a chaotic/integrable transition at small values of κ ∼ 1 (roughly compatible with the transition measured by the OTOCs) as well as that it suggests that the transition could be temperature-dependent. This second observation leads us to conclude that the chaotic/integrable transition does not affect the spectrum homogeneously, but it mostly affects the states in the tail of the spectrum.
In Section 5, we will test this possibility by studying three diagnostics of chaos: the level spacing distribution, the Inverse Partecipation Ratios (IPR) and the IR diversity. The latter two are especially interesting since they do not require any unfolding procedure and, more importantly, because they involve also the energy eigenvectors, and not just the energy spectrum. We will see that all these observables point to a clear separation of the chaotic/integrable transitions for the states in the tail of the spectrum and the states in the bulk, with the latter requiring a higher value of κ to move to the integrable dynamics.
In Section 6, we will compute the Thouless time as measured by the connected unfolded SFF. We will confirm that the Thouless time exhibits a clear transition for small values of κ. Transition that seems to be sharper and sharper when N increases.
In Section 7, we will discuss the behavior of the Gaussian filtered SFF for different values of N and κ. We will see, in agreement with the results of [58], that its behavior is, at the level of precision we can achieve, almost independent on N . Moreover, we will see that it detects clearly a transition to the integrable dynamics for large values of κ, approximately κ ∼ 15.
In the discussion, we will present our conclusions and prospects of future work. In Appendix A, we will support the evidences of Section 5 by computing ther-parameter statistics. We will also discuss some puzzles related with the unfolding procedure that we hope to clarify in the future. 6 We will analyze the following model, introduced and studied in [61], which generalizes the SYK hamiltonian adding a quadratic piece (a so-called random mass term) (2.1) As usual, the χ i are N Majorana fermions in 1 dimension satisfying {χ i , χ j } = δ ij and the coupling constants J ijkl are random variables extracted from a Gaussian distribution of mean value 0 and variance 6 N 3 . The random mass term contains the antisymmetric tensor K ij , whose entries are randomly extracted from a Gaussian of mean value 0 and variance 1 N . Finally, we have the parameter κ which controls the chaotic properties of the model: for κ = 0 the model is equivalent to the standard SYK model and hence it is maximally chaotic. As we mentioned in the Introduction, turning on κ the model undergoes a chaotic/integrable transition which, however, is detected differently by the OTOCs and the RMT analysis: it appears around κ ∼ 1 when we consider the OTOCs and around κ ∼ 15 when we consider RMT tools, [61].
Since our goal is to find RMT diagnostics which can be, in principle, sensitive to the scrambling physics -and hence to the type of chaos that the OTOCs diagnose -we will look for RMT quantities which are mostly affected by the tail of the spectrum. To this purpose, we will propose the connected unfolded Spectral Form Factor (SFF) as an interesting RMT quantity to study. We will show that it matches more closely the behavior of the OTOCs, even if we stress that a very precise comparison between the OTOCs results and the results we will obtain in this paper would require to extrapolate our analysis to the large N limit, a problem that we would like to address in the future. Instead, our focus in this paper will be mostly on the behavior of the connected unfolded SFF as a function of κ. We will focus most of our computations on the case of N = 30, but we will also comment some features and caveats of the N = 28 and N = 32 cases. 8 8 As we will see, the case of N = 30 has a better behavior at small values of κ with respect to the cases of N = 28 and N = 32. The reason for this, which will be discussed extensively in Appendix A, is that the quadratic deformation in (2.1) forces the RMT ensemble to fall in the GUE class. When N = 30 also the pure SYK model, with κ = 0, belongs to GUE, while the cases of N = 28 and N = 32 at κ = 0 belong to GSE and GOE, respectively. In this way, by turning on κ we do not break any symmetry and we do not change the RMT ensemble, obtaining a more homogeneous behavior of the Thouless time for small κ.

The Thouless time t Th
In the characterization of the chaotic properties of a system via the OTOCs, the strength of chaos is usually quantified by the Lyapunov exponents: systems with higher values of the Lyapunov exponents are more chaotic than systems with lower (or vanishing) values of the Lyapunov exponents. Moreover, there is a notion of maximally chaotic systems: they are systems with Lyapunov exponents saturating the chaos bound with β the inverse temperature. These quantitative features of chaos are not usually discussed in the RMT approach: in this framework one is interested in studying how much the spectrum of the model under investigation resembles the spectra occuring in RMT. Hence, quantities that are usually studied in this approach are the level repulsion and the spectral rigidity, for example; whereas the concept of maximal chaos and Lyapunov exponent (or the associated scrambling time) are not considered in this case. Given the implications that the scrambling time and maximal chaos have in holography, it would be very interesting to introduce these concepts also in the RMT approach to chaos.
Given these motivations, we conjecture that some remnants of the early-time chaos (as diagnosed by OTOCs and Lyapunov exponents) could be analyzed in the RMT regime by a precise characterization of the so-called Thouless time. The Thouless time is defined as the time after which the SFF of the system under study behaves like the SFF of RMT and we conjecture that systems with higher Lyapunov exponents (which are more chaotic and scramble earlier) should show RMT behaviors at earlier times with respect to systems which are less chaotic. Hence, our goal is to compute numerically the Thouless time for the model (2.1) at different values of the variable κ, with the aim of testing whether the chaotic/integrable transition as seen by this observable is compatible with the transition seen by the OTOC.
Let us emphasize that, very recently, [58] has argued that the time in which the RMT behavior appears is not the scrambling time, in agreement with the argument of [59]. However, we are suggesting the weaker possibility that, even though the time we will measure is not the scrambling time, it could be functionally related to the scrambling time -in a perhaps highly non-trivial way -and so it could be a useful indirect probe of the scrambling time and of the scrambling physics.
A way to evaluate the Thouless time has been discussed, for a different model, in [62]. The computation goes as follows: one computes the connected piece of the SFF using the unfolded spectrum. Concretely, by connected unfolded SFF we will mean the following quantity: where we denoted by {Ẽ 1 ,Ẽ 2 ,Ẽ 2 N/2 } the unfolded energy levels (we will denote the original energy levels with {E 1 , E 2 , E 2 N/2 }) and the symbol . . . denotes the average over many ensemble realizations. 9 The time at which the connected unfolded SFF develops the ramp typical of RMT is the Thouless time.
Let us discuss, via an example, the importance of considering the connected unfolded SFF. In Figure 1 we report the full SFFs, computed for the case of N = 30 and with increasing values of κ, before and after the unfolding procedure. We recognize that, without doing the unfolding, the time at which the ramp starts to appear is decreasing as a function of κ and this is not the behavior we expect. The models with higher values of κ are less chaotic (as expected) because the ramp are more horizontal, but it would be hard to relate this feature to the Thouless time. After performing the unfolding, we see that the time at which the ramp appears is closer to be an increasing function of κ (as expected) and also the ramps show the same slope but, still, at low values of κ it is not possible to find a monotonically increasing trend.
The necessity of considering the connected piece of the SFF is easy to understand: as already discussed in [55], the full SFF is dominated by the disconnected piece at early times, which has a slow fall-off (it goes like 1/t 3 in SYK) and "hides" the ramp behavior of the connected SFF. The typical way to cure this problem simply consists in removing the disconnected piece of the SFF and consider the connected piece and this is the strategy we will follow in this paper.
On the other hand, it has been suggested in [58] that, in many-body quantum systems like SYK, removing the disconnected piece could be not enough: for this kind of systems, also the connected SFF shows a fall-off which goes like 1/t 3 and, hence, it still hides part of the ramp behavior. For this reasons, a different observable has been proposed [58, eq. (19)], which can be viewed as a "Gaussian filtered" SFF with α a free parameter and without unfolding the spectrum. We can intuitively think about this observable as the SFF computed for a spectrum in which the density of the states ρ(E) has been rescaled to e −α E 2 ρ(E) and it is immediate to see that such a rescaling increases the importance of the states of the bulk of the spectrum (we are assuming that the spectrum is centered at E = 0) with respect to the states in the tail. They show that such an observable is more suitable for many-body chaotic systems like SYK, in the sense that it is able to detect the onset of RMT behavior at earlier times than the connected SFF. Of course, this observable by construction is almost insensitive to changes of the spectrum in the tail but it has the nice feature that it is not necessary to unfold the spectrum. 10 In Section 7, we will compare the results, for the values of κ we will probe, for both the connected unfolded SFF and for the Gaussian filtered SFF and we will give an explanation of their different behaviors as functions of κ.

Numerical results on the connected unfolded SFF
Following the suggestion outlined in [58, App. A] we used a fifth order polynomial fitting function to unfold each ensemble realization separately.
We computed the connected unfolded SFFs for N = 28, 30, 32 and for the cases with κ = 0, 0.1, 0.2, 0.35, 0.5, 0.75, 1, 2, 5, 10, 15, 25 as well as the cases with higher values of κ (40 ≤ κ ≤ 90). We have worked at β = 0, 5 × 10 −5 , 10 −4 , 10 −3 . The results for κ ≤ 25 are summarized in Figure 2 and 3, while those for κ ≥ 40 are displayed in Figure 4. For β = 0 we immediately see that the oscillatory behaviors, typical of the SFFs at very high temperatures, make very hard to identify the Thouless times. For this reason, in the following we will mostly focus our attention on the higher values of β.     Moving to the case β = 5 × 10 −5 , 10 −4 , some more interesting behaviors appear. We recognize that the expected behavior of the Thouless time clearly shows up. For κ sufficiently small (say, κ ≤ 1) the time at which the RMT behavior appears is clearly an increasing function of κ: in this sense, we can say that the recipe based on the connected unfolded SFF is working as expected and that the time at which the RMT behavior becomes manifest in the connected unfolded SFF is sensitive to the value of κ. From this point of view, we see that the connected unfolded SFF "cures" the problems we pointed out in Figure 1.
Another feature which is worth to remark is that, especially for N = 30, 32, the connected unfolded SFF is almost unchanged when passing from κ = 1 to κ = 10 and that, for these values of κ, the hints of chaos in the graph are very low: this suggests that the transition from chaos to integrability happens for κ ∼ 1 and that, after the transition, the connected unfolded SFF is insensitive to the value of κ. This intuition is confirmed by Figure 4, in which we plot the connected unfolded SFF for large values of κ. We see clearly that, for these values of κ, we cannot find any remaining hints of chaos in the graph and that everything is unchanged when κ gets increased.
From the results for β = 10 −3 , an interesting point which is important to observe is that the chaotic features at κ ≥ 1, that at β = 10 −4 were still vaguely visible, are now completely washed out: at these values of β and κ the unfolded connected SFF is already showing the integrable behavior. On the other hand, we notice that, for N = 30 and N = 32 the ramps of the connected unfolded SFF do not have exactly the same slope. Indeed the various slopes are slightly different from each other, a behavior that we already noticed in a much more prominent way in the left Figure 1 before the unfolding procedure: we believe that this behavior is signaling that our unfolding procedure is not entirely correct in the tail of the spectrum and that as a result the slopes show some differences.
All together, these preliminary observations suggest that the connected unfolded SFF is able to see the transition from the chaotic to integrable regime in the mass-deformed SYK model: when the relevant parameter κ gets increased the Thouless time also increases accordingly. Moreover, the transition from chaos to integrability seems to be temperature-dependent and at β = 5 × 10 −5 , 10 −4 happens around κ ∼ 1, while it happens for smaller values of κ at β = 10 −3 . Let us emphasize that the value at which the chaotic/integrable transition appears, κ ∼ 1, is compatible with the transition which in [61] was observed by studying the OTOCs. This is different from the transition observed by the RMT analysis which instead appeared at higher values of κ, say κ ∼ 15 or even bigger. However, we stress that our results have several limitations, that we already mentioned in the Introduction and that we will discuss more in Section 8.
In the next Section, we will confirm and will provide an explanation to these features.

Chaoticity depends on energy levels
In this Section we will argue that the chaotic/integrable transition is not affecting the energy levels homogeneously but it start to affect the energy levels close to the ground state -the tail of the spectrum -and then, i.e. for larger values of κ, it moves to the bulk of the spectrum. This pattern is suggested by the results of Section 4, in which we have described that the connected unfolded SFF is seeing a transition from chaotic to integrable regime for values of κ ∼ 1 and that the precise value of κ at which this transition appears is temperature-dependent. Let us recall how the value of β affects the SFF (we will consider, for simplicity, the expression for the full, not-unfolded, SFF). Given a certain energy spectrum {E 1 , E 2 , . . . , E N }, the SFF can be written as From (5.1) it is clear that the inverse temperature β plays the role of a cut-off, telling us how many energy levels effectively contribute to the SFF: indeed, when the products β (E i + E j ) become large, the contributions to the SFF are exponentially suppressed. Combining this observation with the results of Section 4 -which showed that the chaotic hints in the connected SFF for κ ∼ 1 disappear when passing from β = 10 −4 to β = 10 −3 -it becomes natural to guess that the transition from chaos to integrability starts to affect the energy levels close to the ground state and that later (i.e., for higher values of κ) move to the bulk of the spectrum. This would explain the behavior, as a function of the temperature, of the connected unfolded SFF for κ ∼ 1: at β = 10 −4 the SFF is still probing a certain amount of the states in the bulk of the spectrum, which are in the chaotic regime for κ ∼ 1; on the other hand, for β = 10 −3 the SFF is dominated by states which are closer to the ground state where, at κ ∼ 1 the transition to the integrable regime would have already taken place.
To confirm the behavior just described we will compute three different diagnostics of quantum chaos in RMT: the level spacing distribution, the Inverse Participation Ratios (IPR) and the IR diversity (and the structural entropy). We will compute them, at various values of κ, for different portions of the spectrum, obtained by keeping only a finite number of levels, L, ordered from the ground state. In this way, we will see that the chaotic/integrable transition arises at lower values of κ for small values of L, while appears at larger values of κ when L is large.
We want to stress that, contrary to the level spacing distribution, IPR and IR do not require any unfolding procedure. The fact that they also show a separation in the chaotic/integrable transition between the tail and the bulk suggests that this separation is not an artifact due to the unfolding. Moreover they could be more easily related to the physics of the OTOCs, since their definition involves also the energy eigenvectors and not just the eigenvalues.

Level spacing distribution
We studied the level spacing distribution 11 for the unfolded upper block of the N = 30 hamiltonians averaged over 312 ensembles, including only the first L = 100, 500, 1000, 5000, 10000 ordered energy levels, for the values of κ = 0, 0.5, 1, 5. Note that here we unfold the whole spectrum first, and then keep the first L eigenvalues in the unfolded spectrum. This is what the factor e −β (Ẽ i +Ẽ j ) in the connected unfolded SFF (5.1) effectively does. The results are presented in Figure 5.
We clearly see that, for κ = 0, the spectrum shows level repulsion even for energy levels 11 We caution the reader that the SFF and the level spacing distribution are actually probing two, a priori, different manifestations of quantum chaos: the level spacing distribution is a short-range observable, while the SFF is mostly affected by the global structure of the spectrum (it is a long-range observable). However, in typical systems both the long-range and the short-range observables show a similar behavior, chaotic or integrable.  This result confirms our hypothesis that the non-chaotic features start to appear close to the ground state and that, in passing from β = 10 −4 to β = 10 −3 the SFF probes regions of the spectrum where at k = 1, 5 a transition between chaotic and integrable behaviors is happening. Indeed, after the unfolding, we remind that the SFF effectively probes around 1000 energy levels at β = 10 −3 , and around 10000 energy levels at β = 10 −4 , and we see from Figure 5 that these are the values of L at which the transition is taking place for κ ∼ 1.
Another conclusion we reach is that, given the qualitative agreement between the transition from chaos to integrable behavior as seen by the OTOCs and by the connected unfolded SFF, also the behavior of the OTOCs is mostly controlled by the tail of the spectrum, while it is pretty insensitive to the features of the spectrum in the bulk.
In Appendix A, we also study a dimensionless quantity, ther-parameter, which characterizes the level spacing statistics. The averagedr-parameter also supports the distinct transition at κ ∼ 1 and 10 for the tails and the bulk of the spectrum, respectively. However, we have to say that such separation is less evident (but still visible) for ther-parameter computed without performing the unfolding (Figure 21).

Inverse Participation Ratios (IPR)
To provide further evidences to the picture we described in Section 5.1, we compute other chaos diagnostics, inverse participation ratio (IPR) and IR entropy, as defined in [64][65][66][67][68]. In particular, [68] analyzed the IR diversity, which is the normalized exponential of the IR entropy, of the SYK model. These quantities measure the randomness, or the delocalization, of the energy eigenvectors.
Let us consider a normalized energy eigenstate |E n (n = 1, 2, · · · , D): where {|e m } is a set of given basis vectors and D is the dimension of Hilbert space. Then, the inverse participation ratios (IPR) ξ n is defined by In addition, the IR entropy S IR n is given by and it is convenient to evaluate the IR diversity [68]: Note that the IR diversity Ω n and the IPR ξ n correspond to the first Rényi entropy S n and (the exponential of) the second Rényi entropy S (2) n , respectively, defined by These observables measure how widely an energy eigenstate spreads over a given basis. Therefore, they depend on the basis but nevertheless they have played a successful role in characterizing quantum chaos. Of course, these quantities are useful if a "natural" basis is at disposal in the system under investigation. For example, if one choose the energy eigenstates as a basis, the energy eigenstate would be completely localized with respect to that basis, a result which is not very satisfactory. On the other hand, by picking up the site basis in a lattice system, one can study how the system is spatially localized.
In this Section, we focus on the IPR of the mass-deformed SYK model, and we will present a qualitatively similar result for the IR diversity in Section 5.3. The IPRs of eigenvectors of a random matrix drawn from GUE, GOE and GSE form a Gaussian distribution around the mean value 0.500, 0.333 and 0.500 with deviation ≈ 0.00394, 0.00423 and 0.00385, respectively (See Figure 6). Here, we evaluated a normalized IPR, where the normalization factor is given by the dimension of the eigenvector D = 2 14 . Of particular significance are the values of the standard deviations: they are very small, which means that all the IPRs are very close to the average values.
For the IPR of the mass-deformed SYK model, one can naturally choose the "spin-chain" basis where our hamiltonian corresponds to all-to-all random interactions among N/2 sites. Then, the IPR (and, the IR diversity) will measure the "localization" with respect to this basis. We focus on the N = 30 mass-deformed SYK, which clearly exhibits the distinct chaotic/integrable transition of the low-lying energy subsectors as we have seen in the level spacing distribution of Section 4 (also, see ther-parameter analysis of Appendix A). Moreover, we analyze the even parity subsector, and therefore we normalize the IPRs of the energy eigenvectors by D = 2 14 .
We compute the IPR of N = 30 mass-deformed SYK model with various κ's. In particular, we present the distribution of the IPR for κ = 0, 1, 10, 100 with respect to the normalized energy in Figure 7a. For small κ 1, the IPRs are sharply concentrated around ξ ≈ 0.5 over the whole spectrum. As expected for N = 30 SYK model, this distribution is close to the that of the GUE. As κ is increased (e.g., 0 < κ < 10), the IPRs of the low-lying spectrum are decreasing and dispersing while the bulk of spectrum keeps staying around ξ ≈ 0.5 (See Figure 8). Once again, of particular importance are the values of the variances more than the mean values. Indeed, it is from the values of the variances that we can really measure a departure from the RMT behavior.
This implies that the energy eigenstates of the low energy spectrum become "more localized" with respect to the basis at κ ∼ 1 and move away from the delocalization regime in 1 < κ < 10. On the other hand, the bulk of the spectrum still belongs to delocalization regime of the GU E even at κ = 10. For κ > 10, the bulk of spectrum begins to deviate from the delocalization regime.
The distinct transitions of the IPRs of the tails and the bulk of spectrum can also be seen in the mean value of the IPR versus κ in Figure 7b. For full spectrum, the transition of the IPR takes place at κ ∼ 10. On the other hand, the low-lying spectrum (L = 50, 100) cross the transition at κ ∼ 1. However, we stress that from the mean value of the IPR one cannot immediately conclude that a chaotic/integrable transition is taking place, since it could be also a transition from GUE to GOE: what really discriminates between these two cases are the variances of the distributions, plotted in Figure 8. 12 The behavior we described is consistent with the fact that the energy eigenvectors of the tails of the spectrum are typically more sensitive to a perturbation than the bulk of the spectrum. 13 This difference in the delocalization/localization transitions of the low-lying and the bulk of spectrum is also consistent with the previous results based on the spectral observable, the 12 In Section 5.3 we will see that also the structural entropy allows to distinguish between these two extreme cases. 13 We would like to thank Subhro Bhattacharjee for pointing out this. level spacing distribution, of Section 5.1 (as well as ther-parameter in Appendix A). The IPR analysis of the distinct transitions is based on the eigenvectors, which provides independent confirmation of our hypothesis on the chaotic/integrable transition. Note that correlation functions, in particular the OTOCs, also depend not only on the spectrum but also on the eigenvectors. In particular, since the low-lying spectrum have a significant influence on the OTOCs, the Lyapunov exponent should be greatly affected by the chaotic/integrable transition of the IPRs of the low-lying spectrum rather than the bulk one. Hence, we believe that this result provides a further supporting evidence that the connected unfolded SFF, for sufficiently large values of β, indeed sees the same physics of the OTOCs.

IR diversity and structural entropy
Now we present similar results for the IR diversity, as defined in (5.5), for both RMT (Figure 9) and the even parity sector of N = 30 mass-deformed SYK model ( Figure 10). Like the IPR in Section 5.2, we observed the delocalization/localization transition of the tails of the spectrum of N = 30 SYK model at κ ∼ 1 while it takes place at κ ∼ 10 for whole spectrum (see Figure 10 and Figure 11).
Although the IPR and the IR entropy (IR diversity) are independent observables, they seemingly show a qualitatively similar behavior. To study more quantitatively the difference between the IPR and IR entropy we compute another observable: the so-called structural entropy [69] defined by the difference between the first and the second Rényi entropy: The structural entropy is non-negative and bounded by the normalized IPR or the so-called spatial filling factor [69], i.e., Note that the structural entropy 14 is 0 if and only if the corresponding eigenstate is distributed uniformly over a subset of the given basis (i.e., all non-zero |ψ mn | 2 are equal). This can easily be seen by the Jensen's inequality with the definitions of IPR in (5.3) and IR entropy in (5.5) [69]. The structural entropy characterize the shape of the distribution of the probability |ψ mn | 2 and has been applied in quantum chemistry and disordered system etc. (see [70] and reference therein).
We also present the structural entropy of the GUE, GOE and GSE in Figure 12a and N = 30 SYK model in Figure 12b.
From the structural entropy, one can distinguish the integrable regime from the GOE class. Recall that the mean value of the IPR and the IR diversity of the integrable regime (i.e., large κ=0. κ=1. κ=10. κ=50. κ=100. κ) of the SYK model is close to those of the GOE although the shape of the distributions are different and hence, simply from the mean value of the IPRs, one cannot distinguish the integrable behavior from the GOE behavior. But, the averaged structural entropy of the integrable regime (e.g., κ = 1000) is 0.482 (See Figure 12c), which is far from the structural entropy ≈ 0.369 of the GOE class (See Figure 12a). The structural entropy also show the different chaotic/integrable transitions of the tails and the bulk of spectrum (See Figure 13).    Having obtained a qualitative understanding on the behavior of the Thouless time as computed by the connected unfolded SFF in Section 4, and having confirmed these preliminary observations by the arguments in Section 5, we now want to characterize in a more quantitative way the Thouless time as a function of κ. We restrict in this Section the analysis to the values κ ≤ 15 and omit κ ≥ 25 since, as we already observed, these are the values in which the connected unfolded SFF sees the transition from chaos to integrable behavior around κ ∼ 1.
As we stressed this is also the range in which the chaotic to integrable transition is seen by studying the OTOCs, as shown by [61, Figure 3].
The results in Figure 2    We have also estimated the error by the maximal/minimal values of t where the connected, unfolded SFF is valued 1.1 × (minima), following a similar criterion to the one used in [55]. The  results are shown in Figure 14. As we can see, the determination of the Thouless time looks accurate for low values of κ, except the case with N = 30 and β = 5 × 10 −5 where we could not determine t TH for κ ≤ 0.75 since the local minimum is hidden by the large oscillation. For large values of κ unfortunately the error bars are large, because, for κ = 0.75, 1, the connected unfolded SFF around the dip shows an almost flat behavior and so it is hard to determine the true minimum. It is possible that, by taking more samples, these error bars could be reduced a bit but the flatness of the connected unfolded SFF in that area prevents the error bars to 26 become small even for very large ensemble realizations. On the other hand, it is possible that increasing N and going to the large N limit, the flatness could be reduced.
Interestingly we see that, especially at β = 10 −4 , the transition from the chaotic to integrable behavior becomes sharper and sharper when N gets increased. This is an expected feature which the determination of the Thouless time using the connected unfolded SFF seems to be able to capture.

The Gaussian filtered SFF
In this Section, we want to compute the Gaussian filtered SFF, as defined in eq. (3.3), for the mass deformed SYK model, to have a comparison with the results we obtained for the connected unfolded SFF.
We remind that the Gaussian filtered SFF has been introduced in [58] as a variation of the full not-unfolded SFF, in which the contributions coming from the tails of the spectrum are exponentially suppressed with respect to the contributions coming from the bulk of the spectrum. It has been shown in [58] that this procedure allows to see the onset of RMT behavior at earlier times than the standard SFF. It is not clear that the time at which the RMT behavior shows up in the Gaussian filtered SFF is the true time at which the RMT physics appear in the SYK dynamics, but at least this time provides an upper bound to the true time at which the RMT dynamics shows up.
As we already mentioned, the Gaussian filtered SFF, Y (α, t), depends on an additional parameter α, which controls the dimension of the energy window which gets magnified by the Gaussian filtered SFF. The authors of [58] have found an optimum value of α for the SYK model corresponding to α = 2.9 and, in the study of the Gaussian filtered SFF for the mass deformed SYK model, we have taken the same value of α.
In Figure 15 we have reported the Gaussian filtered SFF for the cases of N = 28, 30 and 32, averaged over 603, 312 and 150 ensemble realizations, respectively, at various values of κ. It is immediate to recognize that the Gaussian filtered SFF is not seeing the transition from chaos to integrable dynamics at the low values of κ: indeed, from Figure 15, we see that the plot is pretty insensitive to the value of κ when we are in the regime κ ≤ 1.
On the other hand, we see that the transition from chaos to integrable dynamics is welldescribed by the Gaussian filtered SFF at higher values of κ: starting from κ = 2, we see that the graph displays, after the initial Gaussian fall-off, an intermediate range in which the behavior is not very clear. In total, we see that the RMT behavior, still present at κ = 2, shows up later, a signal that the chaotic features of the model are becoming weaker. The RMT behavior continues to be reduced when we further increase κ and we finally see the transition to the full integrable dynamics around κ ∼ 15. We recall that this value of κ for the transition was already observed in [61] from the RMT analysis, which therefore roughly agrees with the results of the Gaussian filtered SFF. By contrast, we have seen in the previous Sections that the connected unfolded SFF is seeing a transition at lower values of κ, κ ∼ 1, and that this value is in rough agreement with the results from the OTOCs.
Another interesting feature to observe, which agrees with the observations of [58], is that the behavior of the Gaussian filtered SFF is almost unchanged when N increases.
We can understand this discrepancy between the connected unfolded SFF and the Gaussian filtered SFF as a consequence of the way in which the integrable behavior enters in the system as a function of κ, that we have described in Section 5: we have seen that the integrable behavior affects first the states closed to the ground states and that later (for higher values of κ) affects the states in the bulk. Since the Gaussian filtered SFF is, by construction, mostly controlled by the states in the bulk of the spectrum, it sees the integrable transition at high values of κ, while it is pretty insensitive to the changes in the states close to the ground states, which start to appear at low values of κ. The connected unfolded SFF does exactly the opposite: it is mostly controlled by the states close to the ground state and so it sees the transition at lower values of κ. This transition is also the transition that seems to control the behavior of the OTOCs. On the other hand, it does not show any dependence on κ when this is large.
Summarizing, we see that the connected unfolded SFF and the Gaussian filtered SFF -at least for the particular case of the mass-deformed SYK model and as functions of κ -seem to be complementary: they are dominated by different portions of the spectrum and so they show different features of the transition from chaotic to integrable behavior.

Discussion
In this paper, we have studied the onset of RMT dynamics in a mass-deformed version of the SYK model [61], where the usual SYK Hamiltonian with quartic interaction is perturbed by a quadratic term (a random mass-term). The strength of the quadratic deformation is controlled by a coupling constant κ. As observed already in [61], when κ = 0 the model is equivalent to the standard SYK model and so it is maximaly chaotic and shows RMT behavior. On the other hand, when κ is very large, the model is dominated by the quadratic term and hence it is not chaotic. It has been further shown in [61] that this chaotic/integrable transition happens at certain finite value of κ. However, the precise value of κ at which the transition takes place is dependent on the probe one is using to test chaos: using the OTOCs the transition becomes evident for small values of κ (κ ∼ 1), while the RMT analysis performed in [61] (like the level spacing distribution using large portions of the spectrum) detects the chaotic/integrable transition for larger values of κ, around κ ∼ 15.
With the aim of clarifying this point, we studied the chaotic/integrable transition as seen by other two probes of quantum chaos: the connected unfolded SFF and the newly introduced [58] Gaussian-filtered SFF. These two observables have the property to probe different regions of the spectrum: the connected unfolded SFF, at sufficiently low values of the temperature, probe the tail of the spectrum -i.e. the states closer to the ground state -while the Gaussian-filtered SFF by construction is focused in the bulk of the spectrum. We have checked that these two observables also see the transition to the integrable regime at different values of the coupling constant κ: the connected unfolded SFF sees a sharp transition of the Thouless time at small values of κ, which qualitatively agrees with the behavior seen by the OTOCs. On the other hand, the Gaussian-filtered SFF sees the transition to the integrable regime at values of κ in qualitative agreement with the RMT analysis of [61].
This observation suggested us that the chaotic/integrable transition does not take place homogeneously on the entire spectrum, but that states closer to the ground state (in the tail) move to the integrable behavior for lower values of κ and that the connected unfolded SFF (and the OTOCs as well) is mostly controlled by these states. On the other hand, the states in the bulk of the spectrum are less influenced by the quadratic perturbation and so they move to the integrable regime at larger values of κ, and these are the values at which the Gaussianfiltered SFF sees the transtion to the integrable regime. To confirm this picture, we studied the transition to the integrable regime for portions of the spectrum including only L levels starting from the ground states. We performed this confirming study both by analysing the level repulsion (and ther-statistics) and by studying the IPR distribution and the IR diversity (which depend also on the energy eigenstates) and found agreement between all these methods.
Given our results, we are tempted to conjecture that the connected unfolded SFF might be an interesting observable that deserves further studies and that it could be important in order to find a contact point between the early time definition of chaos (based on the OTOCs) and the late time definition of chaos. For example, an intriguing question that one would like to address is the following, already formulated in [41]: is there a way to characterize maximally chaotic systems by looking at their spectrum? It is possible that the answer to this question could be related to the study of the connected unfolded SFF.
Even though the connected unfolded SFF has some intriguing features that we described, we have to say that there are subtleties, already discussed in [58], related to unfolding and that they should be clarified in order to make more robust the computations performed with the connected unfolded SFF. In the unfolding procedure, we have adopted the strategy that we first fit the number of states n(E) = θ(E − E i ) of the spectrum of each ensemble realization by a fifth order polynomial in E, which we call n av (E), to define the unfolded eigenvalues E i = n av (E i ). Though the purpose of unfolding is to remove the global structure of the level distribution and to keep only the local fluctuations, hence n av (E) should almost coincide with n(E) everywhere, in our current method we observe that the deviation is larger at the edges. Indeed we also observed that if we perform the fitting on the edges and bulk separately, we obtain some different results. Hence, our results depend somehow on the way in which the unfolding is realized; even though the checks we performed, some of them are independent on the unfolding, suggest that the general picture should be robust. However, since there is not an a-priori correct choice of fitting Ansatz and hence unfolding is not a well-defined procedure, so far we cannot fully trust our results. In the random matrix model there is a canonical way to define n av (E): we should use the "average density" ρ av (E) ≡ i 1 E−E i and define n av (E) = E −∞ dE ρ av (E ). Since the observables in the (mass-deformed) SYK model are also defined as ensemble averages, it could be another plausible choice to define n av in the same way. 15 Of course, it would be very nice if one could obtain an analytical good approximation of the large N distribution ρ av (E). For this purpose, the results of [71], [72] might be very useful. 15 We thank Antonio García-García and Fidel Schaposnik Massolo for useful discussions on this point.
Once a robust unfolding procedure will be established, it will be very interesting to extrapolate the large N limit of the connected unfolded SFF results. To this purpose, the methods of [73] could be very useful.

Acknowledgement
We would like to thank Subhro Bhattacharjee, Antonio García-García, Guy Gur-Ari, Masanori Hanada, Sungjay Lee, Loganayagam R. and Fidel Schaposnik Massolo for extensive discussions and/or correspondence. We thank Korea Institute for Advanced Study for providing computing resources (KIAS Center for Advanced Computation Abacus System) for this work.

A Ther-parameter statistics
In this Appendix, we compute the distribution of the ratio of the nearest-neighborhood spacing, the so-calledr-parameter statistics [74,75]. We will compute it for the full spectrum of the massdeformed SYK model (at various values of κ) in Appendix A.1 and for subsectors of the massdeformed SYK spectrum in Appendix A.2 (again, at different values of κ). In agreement with the results of Section 5, we will see that the chaotic/integrable transition affect the different portions of the spectrum inhomogeneously. However, we will notice that the results of Appendix A.2 seem to be dependent on the unfolding, in the sense that the separation between the transitions as seen by the low-lying modes and the states in the bulk of the spectrum is more evident when we consider the unfolded spectrum than when we consider the not unfolded spectrum (even though they are still visible also in this case). This is an unexpected feature, since it is known that the results of ther-parameter statistics should be independent on the unfolding (though this statement is usually applied to ther-parameter statistics computed over the full spectrum). We do not know whether this behavior is an artifact due to the fact that our unfolding procedure is not optimal and, also for this reason, we decided to study the IPR and IR diagnostics (in Sections 5.2 and 5.3) to be more confident about the robustness of our findings.

A.1 Ther-parameter of full spectrum
Ther-parameter statistics of SYK(-like) models has been discussed in [37,68,76]. Given the nearest-neighbor energy spacing s n ≡ E n+1 − E n of an ordered energy spectrum {E n }, thẽ r-parameter is defined byr n ≡ min(s n , s n−1 ) max(s n , s n−1 ) = min r n , 1 r n , where we defined r n ≡ sn s n−1 . Note that, by definition, the r-parameter lies in [0, 1]. 16 The Wigner-like surmise for 3 × 3 matrices gives the distribution of the ratio of consecutive level spacing: where β = 1, 2 or 4 for GOE, GUE or GSE, respectively, and Z β is a normalization constant. From the distribution of r's, one can evaluate the mean value ofr-parameter by The mean values ofr for Poisson, GOE, GUE and GSE are summarized in Table 1 [75].
The random matrix classification of SYK models has been studied in [55-57, 77, 78], which is summarized in Table 2. The quartic interaction in our hamiltonian (2.1) belongs to GOE, GUE or GSE class depending on N (See Table 2). On the other hand, the quadratic interaction  plays a role of (random) mass term and can be analytically solved. Hence, it does not show the random matrix universality. Therefore, when the quadratic interaction dominates the quartic one in (2.1) (i.e., for large κ), one can expect that the level spacing exhibits Poisson statistics. On the other hand, for κ = 0, it will shows the corresponding random matrix behavior, and there would be a transition between the random matrix class and the Poisson statistics as κ increases. However, since the quadratic interaction breaks the symmetry of the quartic one for N ≡ 0, 4 (mod 8), the perturbation by the quadratic interaction could lead to additional transition to other matrix class, for sufficiently small 17 κ. This transition between different classes can also be captured by the averagedr-parameters. Hence, we now discussr-parameters as a function of κ.
As already remarked, ther-parameter analysis would not require the unfolding procedure. However, since we have unfolded the spectrum for the connected SFF, we also consider the unfolded spectrum in addition to the spectrum without unfolding. Furthermore, we discuss the difference between ther-parameters of the unfolded spectrum and the spectrum without unfolding.
First, let us consider N ≡ 2, 6 (mod8) where level statistics is expected to be GUE for κ = 0. In this case, the quadratic interaction does not break the symmetry of the quartic one. Indeed, we observe only the chaotic/integrable transition from GUE to Poisson 18 in Figure 16 which becomes sharper and clearer as N increases. Note that the averagedr-parameter for small N (e.g., N = 10 and N = 14) is far from the Poisson statistics. This is a consequence of the fact that the level spacing for the small N does not have enough statistics.
For N ≡ 0 (mod8), the quartic interaction belongs to the GOE class. However, the quadratic one breaks the symmetry of the GOE, which leads to a transition from GOE to GUE at small κ ∼ 10 −2 (see Figure 17). It would be interesting to study the intermediate transition from GOE to GUE in larger N . For instance, one can check whether in large N the intermediate transition  Figure 17: The mean value of ther-parameter for N = 16, 24 ≡ 0 (mod8) of (a) the spectrum without unfolding (b) the unfolded spectrum. The average is performed over M = 2000, 100 ensembles, respectively. We observe transitions from GOE to GUE and from GUE to Poisson.
takes place at non-zero value of κ or not. But, the limitation of numerical computation makes it difficult to perform the large N interpolation of this transition. As κ increases further, one can observe the second transition from GUE to Poisson at κ ∼ 10 due to the dominating quadratic interaction (see Figure 17).
Note that ther-parameters of the spectrum without unfolding and the unfolded spectrum for N ≡ 0, 2, 6 (mod8) do not show crucial difference (See Figure 16 and 17). However, one  can see a huge discrepancy in the mean value of the unfoldedr-parameter for N ≡ 4 (mod8) in Figure 18. Recall that At κ = 0 the quartic interaction is supposed to belong to the GSE class in which we have r ≈ 0.67617. However, ther-parameter of the unfolded spectrum at κ = 0 is 0.65429, 0.450978 and 0.0576467 for N = 12, 20 and 28, respectively (See Figure 19b), which are far from the GSE class. On the other hand, ther-parameter of spectrum without unfolding at κ = 0 is close to the expected value. i.e. 0.67102, 0.670404 and 0.666468 for N = 12, 20, respectively (See Figure 18a).
As soon as the small quadratic interaction appears, it breaks the symmetry of the quartic Hamiltonian, which makes ther-parameter close to zero for both cases (see Figure 18 and 19). As κ increases, ther-parameter increases until it reaches ther-parameter of the GUE class. Then, for example, ther-parameter for N = 28 exhibits the GUE universality from κ ∼ 10 −2 to κ ∼ 10. Moreover, as κ is increased further, we can also observe the transition from GUE to Poisson around κ ∼ 10. Note that ther-parameter of the spectrum without unfolding and the unfolded spectrum exhibit a qualitatively similar behavior for non-zero κ.
Although ther-parameters for small κ < 0.1 show distinct behavior depending on N , they all approach to the GUEr-parameter. In addition, as κ increases further, the transition from GUE to Poisson appears at κ ∼ 10 for every sufficiently large N (See Figure 20). In this chaotic/integrable transition, one cannot find any qualitative difference between the spectrum without unfolding and the unfolded spectrum.  Figure 19: The mean value of ther-parameter for N = 12, 20, 28 ≡ 4 (mod8) of (a) the spectrum without unfolding (b) the unfolded spectrum. The average is performed over M = 4000, 500, 100 ensembles, respectively. We observe a transition from GUE to Poisson.

A.2 Ther-parameter of subsectors of the spectrum
In the previous Section, we focus our attention on ther-parameter for whole spectrum where the chaotic/integrable transition seems to occur around κ ∼ 10. This result is consistent with the Gaussian filtered SFF, introduced in [58] and studied for this particular model in Section 7. However, as we discussed in Section 4, and as we confirmed with the analysis of Section 5, the low-lying spectrum seems to be responsible for the chaotic/integrable transition in the connected unfolded SFF. Hence, we now study ther-parameter of the low-lying spectrum without unfolding as well as the low-lying unfolded spectrum. 19 For the case of N = 30, we observe that the low-lying levels (L = 50, 100) begin to transit to the value of Poisson at smaller value of κ < 10 than the whole spectrum (See Figure 21). Moreover, as we include more levels, the value of κ at which the transition occurs becomes larger. And, for L = 5000, 10000 low-lying levels, the chaotic/integrable transition is close to that of whole spectrum (κ ∼ 10). Note that during the transition, ther-parameter of the lowlying spectrum is affected by the unfolding procedure in contrast to the bulk of the spectrum. Namely, the transition of the low-lying unfolded spectrum seems to be sharper and takes place at lower value of κ ∼ 1 than that of the low-lying spectrum without unfolding.
In Figure 22, we also evaluate the transition ofr-parameter of the same portions of the spectrum for N = 26, 28 as in N = 30 case. For example, we take L = 25, 50, 250, 500, 2500 levels for N = 28 case. Unlike N = 30 case, the low-lying levels for N = 26, 28 do not have enough statistics to show the transition to Poisson. Also, the unfolding procedure lower thẽ r-parameter of the low-lying levels during the transition while the bulk of spectrum is not affected. Furthermore, the effect of the unfolding on the low-lying spectrum becomes more drastic for larger value of N (See Figure 21 and Figure 22).
As we already mentioned in the beginning of this Appendix, this difference in the behavior of the low-lying modes with respects to the bulk of the spectrum, which becomes sharper when N gets large and by using the unfolded spectrum, could be a hint of some interesting large N behaviors of the model. It is possible that this behavior is somehow magnified artificially by the unfolding procedure but given the results of Sections 5.2 and 5.3 we are optimistic about the robustness of our findings.