Adiabatic eigenstate deformations as a sensitive probe for quantum chaos

In the past decades, it was recognized that quantum chaos, which is essential for the emergence of statistical mechanics and thermodynamics, manifests itself in the effective description of the eigenstates of chaotic Hamiltonians through random matrix ensembles and the eigenstate thermalization hypothesis. Standard measures of chaos in quantum many-body systems are level statistics and the spectral form factor. In this work, we show that the norm of the adiabatic gauge potential, the generator of adiabatic deformations between eigenstates, serves as a much more sensitive measure of quantum chaos. We are able to detect transitions from non-ergodic to ergodic behavior at perturbation strengths orders of magnitude smaller than those required for standard measures. Using this alternative probe in two generic classes of spin chains, we show that the chaotic threshold decreases exponentially with system size and that one can immediately detect integrability-breaking (chaotic) perturbations by analyzing infinitesimal perturbations even at the integrable point. In some cases, small integrability-breaking is shown to lead to anomalously slow relaxation of the system, exponentially long in system size.


I. INTRODUCTION
Finding signatures of chaos in the quantum world has been a long-standing puzzle [1][2][3]. In the last few years exciting progress has been made on characterizing the effects of chaos on dynamical properties of quantum manybody systems, see Fig. 1 [4][5][6][7][8][9][10][11]. Classical chaos is usually described in terms of an exponential sensitivity to initial conditions [12]. However, the quantum world precludes any definition of chaos in terms of physical trajectories due to the Heisenberg uncertainty principle. Alternatively, chaos can be defined in terms of the absence of integrability. Classical Liouville-Arnold integrability is formulated in terms of independent Poisson-commuting integrals of motion. Again, although there have been many attempts to characterize quantum integrability in a similar way, no such unique definition exists [13][14][15][16].
In the last two decades, Random Matrix Theory (RMT) [17][18][19] has shown outstanding success in the understanding of quantum chaos. Following the work of Wigner [20,21], Bohigas, Giannoni, and Schmit [22] conjectured that the energy-level statistics of all quantum systems whose classical analogues are chaotic should show level repulsion and belong to one of three universal classes depending upon their symmetry: the Gaussian orthogonal ensemble, the Gaussian unitary ensemble, or the Gaussian symplectic ensemble. On the other hand, according to the Berry-Tabor conjecture [23], integrable systems have uncorrelated energy levels and exhibit Poissonian level spacing statistics. These ideas were later extended to generic quantum systems and tested numerically under the general framework of the eigenstate * dsels@g.harvard.edu In systems without spatial locality this could lead to fast scrambling, allowing one to identify a Lyapunov exponent. Systems with spatial locality are further characterized by an additional, so-called, butterfly velocity. While this ballistic propagation ends at times O(L), diffusive dynamics continues up to the Thouless time O(L 2 ). All local dynamics has now come to a stop, nonetheless operators keep spreading over operator space, becoming increasingly more complex. This process continues for exponentially long times, only stopping at the Heisenberg time exp(S(L)).
thermalization hypothesis (ETH) [24][25][26][27][28][29]. By now, the emergence of the random matrix behavior of quantum eigenstates is an accepted definition of quantum chaos. Numerically, two additional steps are required before one can accurately compare the statistical properties (e.g. through level statistics or the spectral form factor [30,31]) of a particular quantum system with the predictions of RMT: (1) remove any symmetries; and (2) rescale the spectrum, setting the local mean level spacing to unity (also called unfolding the spectrum). Firstly, if symmetries are not removed, energy levels in different symmetry sectors don't have any correlations, so that spectra of chaotic systems can show Poissonian distributions [32,33]. However, finding all symmetries of a manybody Hamiltonian is computationally hard without any physical intuition, since this effectively involves searching for all possible (local) operators that commute with the Hamiltonian. Secondly, there are various methods to unfold the spectrum, and it is known that statistics, especially ones measuring long-range correlations, can be sensitive to the adopted unfolding procedure [34]. Moreover, the procedure can also exhibit finite-size effects. In light of these issues, it is advisable to rather use the ratio of two consecutive level spacings [35,36].
Here we propose an alternative tool to detect chaos in quantum systems, based on the rate of deformations of eigenstates under infinitesimal perturbations. Mathematically, the distance between nearby eigenstates (also known as the Fubini-Study metric [37][38][39][40] ) can be expressed as the Frobenius norm of the so-called adiabatic gauge potential (AGP) [37,[41][42][43], which is exactly the operator that generates such deformations. From ETH, it is straightforward to show that this norm should scale exponentially with the system size in chaotic systems [37]. In this sense, quantum chaos manifests itself through an exponential sensitivity of the eigenstates to infinitesimal perturbations, which can be viewed as an analogue to classical chaos, reflected in the exponential sensitivity of trajectories to such perturbations. Moreover, unlike standard probes of RMT such as the spectral form factor and level statistics, which only depend on the eigenvalues of the Hamiltonian, the AGP norm is sensitive to both the level spacings and the specific kind of adiabatic deformation (perturbation).
We find that the norm of the AGP shows a remarkably different, and extremely sensitive, scaling with system size for integrable and chaotic systems: polynomial versus exponential (Fig. 2). In our method, we do not need to remove any symmetries before computing the AGP norm. We show that one can detect chaos through the sharp crossover between the polynomial and exponential scaling of the norm. The sensitivity of this norm to chaotic perturbations is orders of magnitude greater than that of the aforementioned methods. Using this approach, we find several, previously-unexpected, results for a particular but fairly generic integrable XXZ spin chain with additional small perturbations: i) The strength of the integrability-breaking perturbation scales exponentially down with the system size, much faster than in previous estimates [44,45]; ii) integrabilitybreaking perturbations immediately lead to an exponential scaling of the norm of the AGP, showing that chaotic perturbations can be already detected in the integrable regimes; and iii) in the presence of small integrabilitybreaking terms, the system can exhibit exponentially slow relaxation dynamics, which is similar to the slow dynamics observed in some classical nearly-integrable systems like the Fermi-Pasta-Ulam-Tsingou (FPUT) chain [46][47][48]. We also find that such relaxation dynamics are very different for observables conjugate (see Eq. (2) be-low) to integrable and chaotic directions (perturbations) of the Hamiltonian. We find similar results for an Ising model, where the integrability is broken by introducing a longitudinal field. 8  The connection with relaxation is not surprising, since one representation of the AGP is in terms of the long-time evolution of a local operator conjugate to the coupling. Hence, our results relate to recent studies of information propagation through operator growth in quantum manybody systems [49][50][51], where chaotic and integrable systems are again expected to exhibit qualitatively different behavior (e.g. in operator entanglement [52,53] and Lanczos coefficients [54]). Whereas most of the previous works mainly focused on short-time effects, here we effectively focus on dynamics and operator growth at times that are exponentially long in the system size ( Fig. 1).

II. ADIABATIC GAUGE POTENTIAL
Before proceeding, let us define the adiabatic gauge potential (AGP) and discuss some of its key properties. Given a Hamiltonian H(λ) depending on a parameter λ, the adiabatic evolution of its eigenstates as we vary this parameter is generated by the AGP as (in units with = 1): Using the Hellmann-Feynman theorem, it is easy to see that the matrix elements of the AGP between such eigenstates are given by where ω mn = E m (λ) − E n (λ), ∂ λ H is the operator conjugate to the coupling λ and we have made the dependence on λ implicit. The diagonal elements of A λ can be chosen arbitrarily due to the gauge freedom in defining the phases of eigenstates. A convenient choice consists of setting all diagonal elements equal to zero. For simplicity we will assume there are no degeneracies in the spectrum, but as will be clear shortly, this assumption is not necessary and does not affect any of the results below. We define the L 2 (Frobenius) norm of this operator as: where D is the dimension of the Hilbert space. It follows from ETH that this expression should scale exponentially with the system size in chaotic systems ||A λ || 2 ∼ exp [S], where S is the entropy of the system [37]. Within ETH, the off-diagonal matrix elements of local operators, including ∂ λ H, scale as m|∂ λ H|n ∝ exp[−S/2] [25,28] while the minimum energy gap between states, ω mn , scales as exp [−S]. The scaling of individual matrix elements was already explored in the literature to study the crossover between chaotic and nonergodic behavior, e.g. in the context of disordered systems [55,56]. As we will demonstrate, the exponential scaling of the norm of the AGP can be used to detect the emergence of chaotic behavior in the system with tremendous (exponential) precision.
However, Eq. (2) is not particularly convenient: the norm of the exact AGP can be dominated by the smallest energy difference between eigenstates, and as such it is highly unstable and difficult to analyze, especially close to the ergodicity transition. Accidental degeneracies in the spectrum that are lifted by ∂ λ H also cause the norm to formally be infinite. To resolve this issue, it is convenient to instead define a 'regularized' AGP as follows: where µ is a small energy cutoff. This has a clear physical intuition: instead of considering transitions (matrix elements) between individual eigenstates, we now only consider transitions between energy shells with width µ. For eigenstates with |ω mn | µ, this reproduces the exact AGP, whereas in the limit |ω mn | µ, the AGP no longer diverges but reduces to a constant. Alternatively, within the operator growth representation (see Eq. (8) below), µ −1 has the interpretation of a cutoff time. Numerically, this regularization has the immediate advantage that it gets rid of any problem with (near-)divergences. Note that µ does not need to be system-size independent for this. Interestingly, as long as µ ∝ exp[−S], the norm of the AGP within chaotic systems should also remain proportional to exp[S]. We can use this flexibility in defining µ to our advantage, choosing it to be parametrically larger than the level spacing to eliminate any effect of accidental degeneracies, but still exponentially small to minimize the deviation from the exact AGP. We find that choosing µ(L) ∝ L exp[−S(L)], where L is the system size, is the most convenient choice (see Appendix A).
Within the ETH ansatz, all relevant off-diagonal elements are parametrized by a function f λ (ω), as introduced by M. Srednicki [25], according to with σ nm a random variable with zero mean and unit variance. Combined with Eqs. (3) and (4), this yields for chaotic quantum systems, with |f λ (ω)| 2 the average overĒ of |f λ (Ē, ω)| 2 (see Appendix C). One can thus interpret the norm of the regularized AGP as an integral transform of |f λ (ω)| 2 . Alternatively, it is convenient to rewrite the regularized AGP as a time integral [57][58][59]: where sgn(t) is the sign function and is the operator conjugate to the coupling λ in the Heisenberg representation. The exponential factor exp[−µ|t|] can be seen as a particular choice of a filter function in the context of quasi-adiabatic continuation [60][61][62]. Moreover, it's to be noted that Eq. (8) remains valid for classical systems [57,58] and therefore the scaling of the AGP norm can be used to detect classical chaos, which we leave for future work. From the representation of Eq. (8), it is clear that the inverse of the parameter µ plays the role of a cutoff time, limiting the growth of (∂ λ H)(t) in the operator space. Note that this time is much longer than the ones generally studied in literature (e.g, the time scale characterizing the ballistic propagation of information where v LR is the Lieb-Robinson velocity and L is the system size) [49][50][51]. One of the outcomes of our work is that an exponential sensitivity to detecting the onset of chaos requires access to exponentially long time scales (Fig. 1).

III. NUMERICAL RESULTS
We can now compare with results for the AGP in integrable/non-ergodic models. Specifically, we move to the analysis of the norm of the regularized AGP for a specific integrable XXZ model with open boundary conditions [63][64][65][66][67], whose Hamiltonian is given below: We will now consider the effects of various integrabilitybreaking terms. Although the thermodynamics of the above model can be solved exactly using the Bethe ansatz [63][64][65][66][67], we still don't have access to matrix elements of general local operators n|∂ λ H|m and the exact AGP remains out of reach even in the integrable limit. Consequently, there are also no results on the scaling of the AGP with increasing system size.
For reference, we also analyze an Ising model in the presence of a longitudinal field whose Hamiltonian is given below: where open boundary conditions are chosen for the chaotic Ising model. This model has a trivially-integrable limit at zero longitudinal field h z = 0, which maps to a system of free fermions [68]. In this non-interacting (free) limit, the AGP can be computed analytically [37,69] (see Appendix B). In the presence of the longitudinal field, this model shows a Wigner-Dyson type distribution of the energy level spacings, which is particularly pronounced at the parameters h x = ( √ 5 + 5)/8 and h z = ( √ 5 + 1)/4 [70]. We will use these values when computing the AGP in the chaotic regime.
In Fig. 2, we show the AGP norm scaled by the system size ||A λ || 2 /L [71] for the interacting XXZ model and the Ising model both at the chaotic and non-interacting points. This clearly shows the remarkably different scalings with system size L for chaotic, integrable and free models. For chaotic models, the scaled AGP norm shows the exponential scaling expected from ETH. For the free model, the scaled norm is system-size independent up to exponentially small corrections away from the critical point (see Appendix B). For the integrable XXZ model, the scaled AGP norm shows a nontrivial polynomial scaling: ||A λ || 2 /L ∝ L β . We find that the exponent β is nonuniversal and depends on the choice of the anisotropy ∆ (see Appendix D). We have chosen λ = h x for both the integrable and non-integrable Ising models and λ = ∆ for the XXZ model.
While the exponential scaling of the AGP norm in the chaotic regime and the constant AGP norm in the free model are expected, the polynomial scaling of this norm of the XXZ integrable model is very interesting and leads to non-trivial conclusions. Recently, LeBlond et al. [72] have shown that the matrix elements of local operators in this integrable model are not sparse (as compared to the matrix elements of non-interacting integrable models). The latter implies that Eq. (7) for the AGP norm still applies, where |f λ (ω)| 2 can also be found from the Fourier transform of the symmetric correlation function (see Appendix C). Since we have chosen µ to be exponentially small in the system size and ||A λ || 2 is polynomially (not exponentially) large, the function f λ (ω) must vanish as ω → 0. This behavior is to be contrasted with chaotic systems where at small ω this function saturates at a constant value, in agreement with the Random Matrix Theory [28].

IV. INTEGRABILITY BREAKING
Having established the scaling of the AGP norm in three different regimes, we will move to the analysis of integrability breaking by small perturbations and focus on a more generic XXZ model. As an integrability-breaking term, we choose a magnetic field coupled to a single spin in the middle of the chain, acting as a single-site defect, where (L + 1)/2 stands for the smallest integer greater than or equal to (L + 1)/2. Then we analyze the AGP for the total Hamiltonian as a function of the integrability-breaking parameter d . Interestingly, in Ref. [73] it was argued based on the same model that even a single site defect is sufficient to induce chaos in the thermodynamic limit.
In Appendix E, we analyze an extensive integrability-breaking perturbation by considering H = H XXZ + ∆ 2 V with V = i σ z i+2 σ z i and find the results to be very similar. A challenging question is how quickly chaos emerges when a non-ergodic, or integrable system, is subjected to an integrability-breaking perturbation. In classical systems with few degrees of freedom, it is known from KAM theory that integrable systems are stable against small perturbations [74][75][76]. It is widely believed that chaos is generally induced by infinitesimal perturbations in the thermodynamic limit [44,45,77,78], with the potential exception of many-body localization [79,80], although the precise scaling of the critical perturbation strength with the system size remains an open question. A standard limitation of numerical approaches (using e.g., level statistics or spectral form factor) addressing this question is the small system sizes amenable to simulations, where it is possible to reliably extract the data. We will show here how one can get more bang for one's buck by studying the AGP norm.
In Fig. 3 a) we show the scaling of the norm of the AGP as a function of the system size for different perturbation strengths d . We choose the zero magnetization subspace of the XXZ chain with number of spins up N ↑ = L/2 , where L/2 stands for the largest integer less than or equal to L/2, and for the direction of the AGP we choose λ = ∆, i.e. as in Fig. 2. For the cutoff, we choose µ = LD −1 0 , where D 0 is the dimension of zero magnetization sector. From the figure, we clearly see a sharp crossover in the scaling of the norm of the AGP as a function of system size from the integrable power law behavior to the chaotic exponential behavior. The straight lines are obtained by a least squares fit, with the slope extracted for the largest d and then used for other perturbations. After the best fitting parameters were found, the critical system sizes were obtained for a particular defect energy at which the integrable (polyno-mial) and chaotic (exponential) curves intersect. These values are shown in the inset of Fig. 3 a), showing a clear exponential scaling of the critical perturbation strength with the system size.
Another very peculiar observation is that for small defect energies the slope extracted from the exponential fit, β = 1.28, is almost two times larger than the slope predicted by ETH, β = log(2) ≈ 0.69. Combining this result with Eq. (7), which we highlight works both in integrable and nonintegrable regimes, we conclude that at small ω the function |f λ (ω)| 2 should scale exponentially with the system size. This implies that the system must have exponentially long relaxation times, which are known to exist in classical chaotic systems like the FPUT chain [46][47][48]. Although we cannot rule out the eventual relaxation to the ETH value for system sizes greater than studied, our results here suggest that, while an exponentially small perturbation is sufficient to induce chaos in the system, it takes an exponentially long time for the system to thermalize. In Appendix E, we show that a similar behavior persists if we break the integrability by a small extensive perturbation, here chosen as the second nearest-neighbor Ising interactions. We found the same slope of β ≈ 1.28, ruling out that this scaling is induced by the ultra-local nature of the perturbation in Fig. 3 a). As the defect energy is increased further to large values (in particular, d = 1), we find that the slope of AGP norm's exponential growth reduces again to the ETH value of β ≈ log(2) (see Appendix F).
Consistent results are obtained for the Ising model (11), where one can consider breaking the integrability of the transverse field Ising model (h z = 0) by introducing a small non-zero h z -field, while probing an integrable direction.
The results are shown in Fig. 3 b). Like in the XXZ case, we observe a sharp crossover from the the unperturbed scaling of the AGP norm (cf. Fig. 2) to exponential scaling with an exponent that exceeds the ETH expectation, once again having implications on the long time relaxation of the system. Indeed, from expression (7) one finds that, for sufficiently small µ, For L > L * , where the AGP norm has exponential scaling, the norm becomes where C roughly is the value of the unperturbed AGP norm at L * . Recall that we observed a scaling of the critical perturbation strength like d ∼ e −αL , such that one finds where η = β/α, and κ = β − log(2), and we have neglected all polynomial factors in system size. For the XXZ model, the exponents are η ≈ 1.6 and κ ≈ 0.85 log(2) (cf. caption of Fig. 3). Because |f λ (ω)| 2 is the Fourier transform of the two-point correlation function of ∂ λ H (cf. (C2)), as ω → 0 it is proportional to the relaxation time of the system. Combining these considerations, we see that for the XXZ model we have with both κ and η of O (1). Similarly, for the Ising model, τ ∼ h η z e κL where η ≈ 1.8 and κ ≈ 1.28 log(2) (cf. caption of Fig. 3). We see that the relaxation time increases exponentially with the system size. For large system sizes it can saturate at some L-independent value, which should diverge as d → 0. This would reflect in crossover of the scaling of the AGP norm to the ETH result: ||A λ || 2 ∝ exp[S(L)] = exp[log(2)L]. While this scenario seems likely, we do not see any signatures for such a crossover within our numerics and thus cannot rule out more exotic scenarios for the behavior of the relaxation time with the system size. Moreover, at intermediate system sizes accessible to numerics we see an extremely stable exponential scaling of the AGP norm (and hence of the relaxation time), with the exponent β independent on the strength of the integrability-breaking perturbation as long as it is sufficiently small. To contrast the scaling of the AGP norm with more traditional approaches in Fig. 4, we show the mean ratio of energy level statistics as a function of defect energy for system size L = 16. Given subsequent energy level spacings s n = E n+1 − E n , this ratio is defined as For non-ergodic systems and Poissonian level statistics, r ≈ 0.386, whereas for chaotic systems and Wigner-Dyson statistics r ≈ 0.536. In this model, the average ratio r shows the crossover from non-ergodic to ergodic behavior at * d ∼ 0.1 [81]. This crossover value of d has a very weak dependence on the system size. In comparison, for the same system size L = 16 the AGP norm shows a clear crossover to chaos for a much smaller * d ∼ 10 −3 (see Fig. 3 a) ). For larger system sizes, the gap between the chaos thresholds extracted by these two methods becomes even larger. Moreover, we also estimated the critical perturbation strength using the spectral form factor for the same system size L = 16. Since it, in general, doesn't self-average [82,83], we added disorder to zzcoupling in the Hamiltonian (Eq. (13)) which reduces the sensitivity of this probe to detect chaos. We find * d ∼ 0.1, a value where the level statistics is roughly half way between Poisson and Wigner-Dyson (cf. Fig. 4). Such an estimate is consistent with the one extracted from a different model [84].
We believe that the reason that the AGP norm is so much more sensitive is because it effectively detects the change in the differential of the norm with the system size. The absolute value of the AGP norm at the threshold is still much closer to the integrable value than to the chaotic one. Such a differential is much harder to detect using other measures, e.g. the level spacing ratio, because this crossover is much smoother, and it is harder to define a sharp threshold. In Fig. 5 we show similar results, now choosing to deform the Hamiltonian in the direction of the integrabilitybreaking operator itself, i.e. λ = d for the XXZ chain and λ = h x for the Ising chain. We choose to work in the full Hilbert space with dimension D = 2 L . We find that the AGP norm shows exponential scaling even when d = 0, i.e. when the Hamiltonian is integrable. We find a good fit to the exponential scaling ||A λ || 2 ∼ e βL , with now β ≈ log(2), consistent with ETH. To our surprise, we find that the integrability-breaking direction in the coupling space can be detected already at the integrable point. Again we confirm that the results remain the same if we use an extensive integrability-breaking term instead (cf. Appendix F). Note the difference with the data presented in Fig. 2, where the perturbation itself did not break integrability and polynomial scaling was observed, and the calculation for Fig. 5, where the perturbation does break integrability and exponential scaling is observed: the AGP probes both the Hamiltonian and the perturbation.

V. CONCLUSIONS
We found that the properly-regularized norm of the adiabatic gauge potential, the generator of adiabatic deformations, can serve as an extremely sensitive probe of chaotic behavior. Within chaotic systems, this norm scales exponentially with system size, whereas it scales polynomially in interacting integrable systems and is approximately system-size independent in non-interacting systems for adiabatic deformations preserving integrability. For adiabatic deformations breaking integrability, exponential scaling is generally observed.
Using the present method to investigate the effects of an integrability-breaking perturbation on an interacting and integrable XXZ chain, we found that perturbations that are exponentially small in system size suffice to induce chaotic behavior. We also found that such a small integrability-breaking term leads to anomalously slow dynamics along the integrable directions, with the relaxation time scaling exponentially with system size. Such integrability-breaking perturbations can also be detected at the integrable point, where no anomalous dynamics occur and the relaxation is consistent with the standard diffusive behavior expected from ETH.
This motivates the use of the adiabatic gauge potential, which is connected with both deformations of eigenstates and operator dynamics, as a sensitive probe into either chaotic or integrable behavior of quantum many-body systems.
Unless stated otherwise, in all calculations we have chosen a cutoff µ = LD −1 , where D is the dimension of the Hilbert space. The prefactor L has been chosen to remove the logarithmic correction coming from the zerofrequency contribution of |f (ω = 0)| 2 = L in chaotic models (see Appendix C). This can also be motivated by plotting the AGP norm and comparing it w.r.t. different choices of cutoff. We first study this norm close to chaotic-integrable transition point and then later describe its effect deep in the chaotic regime.
When we are close to the chaotic-integrable transition point and the cutoff is too small (e.g. µ = L −1/2 D −1 ), then we find that the AGP norm is too sensitive to the exponentially close eigenstates, showing a non-smooth exponential scaling, which makes it hard to draw any conclusions (see Fig. 6 a) ). On the other hand, if the cutoff is too large (e.g. µ = L 2 D −1 ), then the AGP norm, albeit smooth, is no longer sensitive to the small strength of integrability-breaking perturbation(see Fig. 6 b)). In Fig. 6 c) with µ = LD −1 , we find that the rescaled AGP norm shows an exponential scaling that is both appropriately smooth and exponentially sensitive to integrabilitybreaking perturbations.
Deep in the chaotic (ergodic) phase, we find that the numerically-obtained scaling for the norm of the AGP is almost the same for the different choices of cutoff scaling we studied.

Appendix B: Derivation of AGP for the free model
As shown in Refs. [37,69], the AGP for changing the transverse field h x in a free Ising model with periodic boundary conditions is given by where the operators O l are given by the following Pauli string operators The norm of the AGP follows as where Tr [O l O p ] = 2 L+1 L was used since all strings of Pauli matrices are trace-orthogonal. The above expression was used to compute the AGP norm for the free model in Fig. 2 in the main text.
To obtain the scaling with system size, we can use the analytical expressions of α l for large enough system sizes [37], i.e. α l = h −l−1 x in the paramagnetic phase where h 2 x > 1. Using this, we find that Recall that the correlation length in the transverse field Ising model ∼ 1/ log h x .

Appendix C: Derivation of the AGP norm for chaotic and integrable systems
As mentioned in the main text [28], following ETH the matrix element of a local operator ∂ λ H can be written as: where ω = E m − E n ,Ē = (E n + E m )/2, σ nm a random variable with zero mean and unit variance, and f λ (ω,Ē) is related to the non-equal time correlation function of ∂ λ H. Its expression at infinite temperature (for details, see Ref. [28]) is given by: where {...} stands for the anti-commutator and connected correlation function The norm of the AGP can now be expressed in terms of the f λ (ω,Ē) function, where |f λ (ω)| 2 = n |f λ (ω,Ē)| 2 /D, the density of states is given by Ω(E n + ω) = e S(En+ω) and at infinite temperature Ω(E n + ω)e −S(Ē) = 1. As explained in the main text, the above expression is true for both interacting integrable and chaotic models. For sufficiently small µ one can assume f λ (ω) to become constant, hence it can be pulled out of the integral and the AGP norm thus becomes: Deep in the ergodic (chaotic) phase, we have |f λ (µ)| 2 = L [28] due to diffusive dynamics at long times and therefore, ||A chaotic || 2 ∼ e S(L) = e L log(2) for µ ∝ Le −S(L) . On the other hand, we find that close to the chaoticintegrable transition point ||A|| 2 ∼ e βL with β ≈ 2 log(2). More specifically, we found that in this regime slopes β/(2 log(2)) = 0.92 for XXZ model and β/(2 log(2)) = 1.14 for Ising model. Note that it follows directly from eq. (C3), and x 2 /(x 2 + 1) 2 ≤ 1/4, that Consequently, for any local perturbation the norm of the regularized AGP -where we set µ ∼ L exp(−S(L))can't grow faster than exp(2S(L)). Not only does it appear that this bound gets saturated when probing integrable direction ∂ λ H in models in which the integrability is weakly broken, it further implies that those observables ∂ λ H take exponentially long to relax. Indeed, the above scaling can only be achieved by effectively having |f λ (µ)| 2 ∼ e S(L) . Yet, the total spectral weight, dω |f λ (ω)| 2 , is only polynomially large in the system size, implying that the corresponding spectral weight must be localized in a region ∆ω ∼ e −S(L) . Combined with expression (C2), the latter implies ∂ λ H(t) takes exponentially long to relax to equilibrium.
For interacting integrable models we found ||A|| 2 ∼ L β , where the exponent β is non-universal. Since the norm is not exponential in system size, the function |f λ (µ)| 2 ∼ e −S(L) . This means that the function should vanish in the zero frequency limit, which (c.f. Eq. (C2)) implies oscillatory dynamics of the observable ∂ λ H(t).

Appendix D: Effects of the anisotropy XXZ model.
In this Appendix, we will again consider the XXZ Hamiltonian (Eq. (10)): where ∆ is the anisotropy, and we take ∆ = λ as the adiabatic deformation, but now at different values of ∆. We find that the slope of the AGP norm depends nontrivially on ∆ (Fig. 7).

Appendix E: NNN interactions in the XXZ chain
In the main text, we studied the effect of strictly local integrability-breaking operator (whose support is a single site). Looking into the effects of the locality, we here study an extensive integrability-breaking operator. We add a next-nearest-neighbor (NNN) interaction to the XXZ chain, with Hamiltonian given as: The above model is chaotic for large enough ∆ 2 [32]. We choose λ = ∆ (Fig. 8) and λ = ∆ 2 (Fig. 9). In the limit ∆ 2 → 0, when the above Hamiltonian (eqn. E1) is integrable, the former (latter) is the integrabilitypreserving (breaking) direction. As shown in Figs. 8 and 9, results are similar as for the strictly local perturbation studied in the main text. This implies our results are robust to the nature of the adiabatic deformation. Rescaled AGP norm ||A λ || 2 /L with λ = ∆ of the XXZ chain at ∆ = 1.1 shows a sharp crossover from polynomial to exponential scaling with system size, even for very small perturbation strengths ∆2. As ∆2 decreases, the system size where this crossover happens increases. Straight lines are the exponential fits with |A λ || 2 /L ∼ e 1.28L . Inset: The integrability-breaking defect energy scales exponentially with system size, i.e. ∆ * 2 ∼ e −0.9L . This is calculated for the symmetry sector with zero magnetization.

Appendix F: Universal slope of the AGP norm
Here we study the AGP norm in the XXZ chain in the limit when the magnitude of the integrability-breaking perturbation (either the local defect energy d or the NNN interaction strength ∆ 2 ) is of the same magnitude as the ∆/J energy scale. In this limit, we find that the AGP has an exponential scaling with system size characterized by an almost universal slope β ≈ log 2, which is close to the one predicted by ETH. Details about the model and its parameters are given in the caption of Fig.  10.