Krylov complexity from integrability to chaos

We apply a notion of quantum complexity, called"Krylov complexity", to study the evolution of systems from integrability to chaos. For this purpose we investigate the integrable XXZ spin chain, enriched with an integrability breaking deformation that allows one to interpolate between integrable and chaotic behavior. K-complexity can act as a probe of the integrable or chaotic nature of the underlying system via its late-time saturation value that is suppressed in the integrable phase and increases as the system is driven to the chaotic phase. We furthermore ascribe the (under-)saturation of the late-time bound to the amount of disorder present in the Lanczos sequence, by mapping the complexity evolution to an auxiliary off-diagonal Anderson hopping model. We compare the late-time saturation of K-complexity in the chaotic phase with that of random matrix ensembles and find that the chaotic system indeed approaches the RMT behavior in the appropriate symmetry class. We investigate the dependence of the results on the two key ingredients of K-complexity: the dynamics of the Hamiltonian and the character of the operator whose time dependence is followed.


Introduction
The notion of "complexity" is playing an increasingly important role in a number of physical contexts [1], from computational condensed matter all the way to holographic spacetime [2]. As suggested by the colloquial meaning of 'complexity', such a quantity should capture the notion of how 'complicated' a physical system is. Quantum mechanically such a notion could refer to states, say with respect to some chosen 'simple' reference state, or operators, or perhaps some combination thereof. In fact, a particularly natural notion of complexity is associated with the time evolution generated by the Hamiltonian itself. A mathematically precise definition of the complexity of time evolution under a given Hamiltonian is given by Krylov complexity [3][4][5] or 'K-complexity' for short. To date, several aspects of Krylov complexity have been studied in various setups and systems, for example [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].
Unitary evolution under a quantum Hamiltonian sends an initial operator O 0 to its Heisenberg-evolved time-dependent version e iHt O 0 e −iHt , exploring thus over time the space spanned by successive commutators of the form [H, [H, · · · [H, [O 0 ]]]. In this way the initial operator O 0 explores a larger and larger subspace of the Hilbert space of operators, and a natural notion of complexity should quantify how quickly this spread occurs, and furthermore how big a subspace of the Hilbert space of operators is eventually explored. As described in section 2 below, K-complexity captures exactly this notion of spread mathematically and allows us to quantitatively distinguish different physical systems by the efficiency of this spread. This is done by transforming the intuitive idea of exploring higher and higher commutators, as above, into an orthogonal basis of the Hilbert space of operators, and studying the Heisenberg dynamics with respect to this basis.
In addition to the behavior of K-complexity for given individual quantum systems, such as the SYK model [3,5,6], 2D CFTs [18,19], and more general symmetry-based Hamiltonian systems [20,21], it is interesting and important to categorize the possible Krylov phenomenologies according to more universal criteria. One of the most interesting of these is clearly the behavior of K-complexity in the class of chaotic quantum systems as opposed to that of integrable ones, initiated in [7] for systems away from the thermodynamic limit.
Quantum integrable systems, such as the strongly interacting XXZ chain [25] or the quadratic SYK model, are less efficient at exploring Krylov space, as evidenced for example by their reaching a lower saturation value of K-complexity at late times. By mapping the dynamics of operator spreading to an off-diagonal Anderson-like hopping problem on the Krylov chain this under-saturation is linked to the (partial) localization of the wave function on the Krylov chain [7]. Maximally chaotic systems feature a late-time complexity saturation value which is exponential in the number of degrees of freedom [5]; interacting integrable systems saturate at quantitatively lower values as compared to chaotic models due to localization effects in Krylov space [7], and free systems typically depict complexity saturation values at late times which are linear, or polynomial, in the number of degrees of freedom [5]. Related work, in the thermodynamic limit, includes [10] for systems featuring many-body localization, and [11] for the transverse-field Ising model with and without an integrability breaking term. We will focus on K-complexity at long time scales for finite systems away from the thermodynamic limit [4,5,7].
In this paper we explore K-complexity in a class of quantum systems that show integrable to chaotic phase transitions as a function of certain control parameters. We also characterise, for the purpose of comparison, the behavior of K-complexity in random matrix theory, both with and without time reversal symmetry. Interestingly we find agreement between the late-time behavior of a quantum chaotic Hamiltonian with that in the random matrix ensemble of the right symmetry class.
In the remainder of this paper we will introduce and review background material regarding K-complexity (Section 2), and summarize what is known about the evolution as a function of time (see Figure 1). Section 3 will introduce the main working horse of this study, the XXZ spin chain as well as two different integrability breaking deformations. In section 4 we will present numerical results performed on the integrable XXZ chain and its chaotic deformations exploring the behavior of K-complexity through the integrabilitychaos transition with particular regard to the late-time saturation value. Section 5 es- Figure 1: This table summarizes the general behaviour expected from K-complexity, particularly for finite chaotic systems with S degrees of freedom (Λ is the bandwidth of the system). It is based on [3,4] and [5].
tablishes the analogous results in pure random matrix theory (RMT), and categorizes the K-complexity behavior of chaotic systems with and without time reversal symmetry. We shall find agreement between the chaotic spin chain and the appropriate RMT universality class at sufficiently late time. We end with a discussion of our results in Section 6.

Review of K-complexity and its late-time saturation value
Krylov complexity is a measure of operator complexification as it evolves in time. It is defined by constructing an orthonormal basis starting with the operator itself and constructing orthogonal directions by iteratively commuting it with the Hamiltonian. It was introduced in [3] as a probe of quantum chaos in the thermodynamic limit, and in [4] it was suggested as a measure of operator complexity at all time scales for finite systems with operators satisfying the Eigenstate Thermalization Hypothesis (ETH) [26][27][28][29][30]. In [5] K-complexity was computed for complex SYK 4 systems and it was shown numerically that its time-dependent profile fits the one expected from quantum computation as well as from holography [2,31]. K-complexity was computed in [7] for the XXZ model which is a strongly interacting many-body integrable system, and was shown to saturate at late times at values below those found for SYK 4 which is a maximally chaotic system [32][33][34]. This paper aims to bridge the gap between the integrable and the chaotic by introducing a Hamiltonian which interpolates between the two, and studying K-complexity for a fixed type of local operator.
We now briefly review the definition of K-complexity. Given a Hamiltonian H, an operator O and an inner product (A|B) = 1 D Tr(A † B) where D is the Hilbert space dimension, the Krylov basis is defined by an iterative orthonormalization procedure known as the Lanczos algorithm: Compute A n If A n = 0: STOP Otherwise: define norm b n = A n and normalized operator O n = A n /b n .
Here, A = (A|A) is the norm of an operator A and it is to be understood that b 0 = 0 and O −1 = 0. In this way we construct a complete ordered orthonormal basis, the Krylov chain, adapted to the operator's time-evolution. The value of n at which the algorithm terminates is the Krylov space dimension, denoted by K, and in [5] it was shown that it satisfies K ≤ D 2 − D + 1. This bound is saturated in all the cases studied in this paper. The orthonormalization coefficients b n are called the Lanczos coefficients.
The time-evolving operator can now be expanded in the Krylov basis where φ n (t) can be thought of as the wavefunction over the Krylov basis, which satisfies, via the Heisenberg equation, a Schrödinger-like equation with boundary conditions φ −1 (t) = 0 and φ n (t = 0) = δ 0n . From Unitarity, since the initial operator is normalized at the first step of the Lanczos algorithm, the wavefunction φ n (t) is normalized at all times: n=0 |φ n (t)| 2 = 1. K-complexity is defined as the time-dependent average position over the Krylov chain 3) The behavior of C K (t) at different time scales for finite (chaotic) systems with S degrees of freedom is summarized in Figure 1. It is associated with the behavior of the Lanczos coefficients at different n scales via the wavefunction (shown schematically in the same figure). The Krylov elements {O n } K−1 n=0 can be thought of as a basis of sites |O n ) in a chain of length K. The action of the Liouvillian L ≡ [H, ] relates different sites on the Krylov chain i.e. L|O n−1 ) = b n |O n ) + b n−1 |O n−2 ) and its matrix is tridiagonal. We shall denote the eigenvalues of the Liouvillian by ω i and its eigenvectors by |ω i ), ( In this formulation, the time-evolution of the operator is given by From (2.1), the wavefunction φ n (t) is the projection (O n |O(t)). Hence the time-dependent transition amplitude is and the long-time average of |φ n (t)| 2 is given by Note that the late-time average of the transition amplitude is normalized, From (2.3) and (2.8) the late-time saturation value of K-complexity is which will be the main object of study of this paper. Before moving on to a more concrete study of K-complexity, we need to clarify the role of the connected and disconnected contributions to the various quantities just introduced.

Effect of operator's trace on late-time saturation value of K-complexity
In this section we discuss the effect of the operator's trace on the late-time saturation value of K-complexity. It can be directly linked to the influence of the disconnected part of the two-point function on its late-time plateau and, just like when studying the latter, in order to probe universal effects due to chaotic behavior, one may work with operators with a zero one-point function or subtract it explicitly if it is initially non-zero. For a hermitian normalized operator written in the energy basis as in (2.5), the two-point function is given by From which Q 00 is obtained by setting n = 0 in (2.8): where in the last step we assumed the absence of degeneracies or rational relations in the energy spectrum. In Appendix A it is shown that using the "ETH" ansatz for RMT [35], i.e.
where O is O(1) and the matrix r ab is drawn from a Gaussian ensemble with zero mean and unit variance, the scaling of Q 00 is as follows (2.14) In the case of the two-point function, to better probe the spectral correlations in the system, one usually studies the connected part which amounts to using the traceless operator:Õ Using the traceless version of the operator (2.15), it is shown in Appendix A that together with the ansatz (2.13), the connected version of Q 00 behaves as (2.16) From (2.10), if Q 00 is significantly large compared to Q 0n for n > 0 (taking (2.9) into account), the saturation value of K-complexity will be pulled down to smaller values. From (2.14) it is clear that Q 00 can be as large as permitted by normalization for an operator with non-zero trace, its value being controlled by the one-point function, but this does not reflect any universal behaviour of the autocorrelation function. In order to study universal features, one must work with operators with a zero one-point function, since in that case the two-point function is directly equal to its connected part, and ETH predicts that the latter plateaus at 1 D , while the corresponding transition probability plateaus at Q (c) 00 ∼ 1 D 2 . In chaotic systems, this is consistent with the observation that Q (c) 0n approaches O 1 D 2 ∼ 1 K for all n = 0, ..., K − 1. This means that for chaotic systems the late-time transition probability is more uniform as a function of n, compatible with the normalization K−1 n=0 Q 0n = 1 and implying a K-complexity long-time average which approaches ∼ K 2 as seen for example in complex SYK 4 [5]. For the effect of the 1-point function on the saturation value of K-complexity in the complex SYK 4 model see Appendix A.1.

XXZ and its integrability breaking
The Heisenberg XXZ spin chain is an integrable model which exhibits Poisson level-spacing statistics. The model consists of nearest-neighbor spin interactions where S α i = 1/2σ α i and σ α i are the Pauli matrices with α = x, y, z. In a series of papers it was shown that even the addition of a local operator such as to the XXZ Hamiltonian can break its integrability [36][37][38][39][40][41] and spectral statistics will show chaotic behaviour. Another type of integrability breaking term [42] we will consider is the next-to-nearest-neighbour operator We will demonstrate the transition from integrability to chaos by studying the distribution of the ratios of consecutive level spacings [43,44], and show that increasing the strength of the integrability breaking term from zero will result in a transition in the spectral behaviour from integrable to chaotic. For an ordered set of energy eigenvalues {E i } D i=1 , consecutive level spacings are defined as s i = E i+1 − E i and consecutive ratios are defined as the set r i = s i /s i−1 . The distribution P (r) was computed in [43] for the random matrix ensembles GOE, GUE and GSE. It is useful to define the quantityr i = min r i , 1 r i with distribution P (r) = 2P (r)θ(r − 1) whose mean r can be used as an indicator to distinguish an integrable system from a chaotic one. For a Poissonian distribution of level-spacings P (s) = e −s , the distribution of r is given by P (r) = (1 + r) −2 [44] and r = 2 ln 2 − 1 ≈ 0.38629. For the Wigner ensembles (GOE, GUE and GSE) distinguished by their Dyson index (β = 1, 2 and 4 respectively) it was shown in [43] that a very good approximation for practical purposes is P

Choice of sector and local operator
The XXZ Hamiltonian commutes with the operator representing the total spin in the zdirection and is invariant under reflection with respect to the edge of the chain, represented by the parity operator P [45]. To avoid degeneracies in the Hamiltonian spectrum we will work in a sector with fixed total spin and parity. To study K-complexity we will use open boundary conditions and focus on a local operator O which respects these two symmetries and keeps the computation within the chosen sector: where i is chosen to be near the center of the chain. This operator has non-zero trace and we remove its trace according to (2.15) before performing the Lanczos algorithm. Appendix A.2 discusses the effect of the one-point function on the saturation value of K-complexity in pure XXZ. When adding the integrability breaking term (3.2), we keep within the sector by using an odd-length chain and situating the impurity at the middle of the chain: The integrability breaking term (3.3) commutes both with M and with P .

r-statistics for XXZ with integrability breaking terms
We now present results for the r-statistics of the following interpolating Hamiltonians: and where H XXZ is given in (3.1), H d in (3.6) and H N N N is given in (3.3). We work with various values of J zz and set J = 1 in all cases. Some of the Hamiltonians and operators used in the numerical computations were constructed using the QuSpin package [46]. The Lanczos algorithm and K-complexity computations were performed using the codes we developed in [5] and [7]. Figures 2a and 2b show ther statistics for the Hamiltonians (3.7) and (3.8) respectively. We plot the distributions ofr for various values of the coefficient of the integrability breaking terms, as well as the mean values r . We compare the results for both P (r) and r with the analytical results for Poisson and GOE mentioned in Section 3. We see that increasing the strength of the integrability breaking term makes the system transition from displaying integrable statistics to displaying chaotic statistics. Note that after the transition, increasing the value of the coefficient of the integrability breaking term even further makes the system less chaotic, as can be seen in Fig

K-complexity and integrability-chaos transition
In [7] it was shown that the saturation value of K-complexity is sensitive to the integrability/chaos of a model, by comparing results for complex SYK 4 systems with results for XXZ systems of similar Krylov space dimensions. It was argued that the time evolution on the Krylov chain given by Equation (2.2) can be mapped to an Anderson problem with off-diagonal disorder. Higher disorder would imply some amount of localization for the Liouvillian eigenvectors and hence a smaller saturation value of K-complexity, while less disorder would imply less localization and higher saturation values of K-complexity. In this section we study the Lanczos coefficient statistics and saturation value of K-complexity for the interpolating Hamiltonians given by (3.7) and (3.8), with an operator of the type (3.5).   Figures 4b and 5b, this is consistent with the phenomenology described in [7] namely that the saturation values of K-complexity will increase with decreasing disorder in the Lanczos coefficients.
By increasing the value of the coefficient of the integrability breaking term we interpolate from a fully integrable model (XXZ) to a chaotic model, as can be seen through the r transition in Figure 2. In Figure 3 we plot the distribution of the log of ratios of consecutive Lanczos coefficients log(b n /b n+1 ). The mean of this distribution is ≈ 0 and the standard deviation generally decreases with the strength of integrability breaking, indicating less disorder in the Lanczos sequence. Indeed, we find consistently that the saturation value of K-complexity is affected by the strength of the integrability breaking term, and generally increases with the value of the integrability breaking coefficient, as can be seen in Figures 4 and 5. The latetime saturation value of K-complexity as a fraction of the Krylov space dimension can be read off from the vertical lines in the figures, where the x-axis was scaled according to the corresponding Krylov space dimension. Another interesting aspect is the timedependent profile of K-complexity at various time scales and for different integrability breaking strength, for which results are presented in Figure 6. Again we find a consistent relationship between the strength of the integrability breaking term and the value of Kcomplexity.

RMT results and dependence on universality class
This Section gathers results on the saturation value of K-complexity for different operators in Random Matrix Theory, to be used as a reference to compare with the results obtained  where DH is a flat measure, the standard deviation σ sets the energy units (and was set to 1 in the numerics), and H is a complex hermitian or a real symmetric matrix depending on whether we work with the Gaussian Unitary Ensemble (GUE) or with the Gaussian Orthogonal Ensemble (GOE), respectively 1 [47].

Influence of the structure of the seed operator
A detailed numerical study reveals that the behavior of K-complexity, and in particular its late-time saturation value, is not only controlled by the statistics of the Hamiltonian spectrum, but also influenced by the structure of the operator under consideration. As an extreme illustration of this, Appendix B shows analytically that an operator that is constant in the energy basis, which is a very atypical observable in any system, features a late-time K-complexity saturation value of ∼ K 2 regardless of the spectrum of the underlying Hamiltonian. In contrast, a typical operator in RMT should satisfy the RMT operator Ansatz (see e.g. [35]) for its matrix elements in the energy basis: where all {r ab } are independent random numbers 2 drawn from a normal distribution with zero mean and unit variance; they are either real or complex depending on the universality class at hand. The one-point function term O in (5.2) will not be important for the current analysis because, as explained in Appendix A, we shall work with traceless operators. Operators satisfying the Ansatz (5.2) can be constructed as sparse operators in the basis in which the Hamiltonian is drawn from the Gaussian ensemble, or as random matrices with independent entries. In both cases, the change-of-basis matrix that brings the operator to the energy basis is a random unitary drawn from the Haar measure and for sufficiently large D they both agree with the structure (5.2). Results on the late-time behavior of K-complexity for both operator choices in the different universality classes can be found in Figure 7, which suggests that the saturation value of K-complexity is sensitive to the universality class to which the Hamiltonian belongs as well as to the choice of operator and, in particular, to whether the operator breaks time-reversal or not. Note that, in general, the complexity saturation values are below K 2 .

Deviations from the RMT Ansatz: ETH operators
In [5] we studied complex SYK 4 , which is a chaotic system with richer features than just RMT and displayed a complexity saturation value close to K 2 . Two features may be regarded as responsible for that behavior: the Wigner-Dyson statistics satisfied by the Hamiltonian spectrum, and the fact that the operators studied satisfy the eigenstate thermalization hypothesis (ETH) [26,29,30,35,48], which is an extension of the RMT Ansatz (5.2) that accounts for a smoothly varying density of states: where r ab are independent (up to hermiticity), identically distributed normal random variables with zero mean and unit variance, and E ≡ (E a + E b )/2 and ω ≡ E a − E b are, respectively, the average energy and the energy difference between the corresponding levels. O(E) and S(E) are the microcanonical one-point function and entropy, respectively, and the function f O (E, ω) gives the Fourier transform of the connected two-point function, sometimes denoted spectral function [3]. Disregarding the E-dependence, the highfrequency tails of this function are known to be bounded from above by an exponential profile: where E T is the Thouless energy, which is itself constrained by a system-dependent upper bound, and controls the regime of applicability of RMT. For the sake of the current analysis, Operator choices: a sparse operator in the basis in which the Hamiltonian is drawn, a random real operator drawn from a Gaussian distribution in the same basis, and a random complex operator again drawn in the same basis. We observe that the random operators saturate at higher values as compared to the sparse operator, and within them, the one that breaks time-reversal (i.e. the complex one) has the highest complexity saturation value. Right: Same choices of operator, but for a Hamiltonian drawn from GUE. In this case time reversal is anyway broken by the Hamiltonian itself, which is why the two random operators have quantitatively very similar features, both having a complexity saturation value slightly higher than that of the sparse operator.
we generated operators following the Ansatz (5.3) where the E-dependence was taken to be constant and the ω-dependence was chosen to saturate the bound (5.4) with an adjustable Thouless energy. The (rescaled) off-diagonal elements in the energy basis {r ab } were chosen to be either real or complex. The saturation value of K-complexity as a function of the Thouless energy for the different choices of Hamiltonian and operator are depicted in Figure  8. The different choices of Hamiltonian and operator can be classified according to how they comport regarding time reversal. If we define the time reversal transformation T as an anti-unitary transformation that acts as complex conjugation T * = K in the basis in which the Hamiltonian is drawn from the Gaussian ensemble, we can make the following identifications: • GOE + real r ab : This situation matches that of a time-reversal preserving operator in a system with a Hamiltonian that preserves T , as they both are real in the computational basis 3 , and therefore the operator will still be real in the energy basis.
• GOE + complex r ab : This case describes the situation in which the Hamiltonian is T - Figure 8: K-complexity saturation value as a function of the Thouless energy for different ETH operators (either real or complex in the energy basis) with Hamiltonians drawn from two different RMT ensembles (GOE and GUE). The horizontal lines mark the asymptotic value for E T → ∞, corresponding to the case of observables satisfying the pure RMT operator Ansatz. In order to mod out system-dependent scaling of the energy spectrum due to the choice of normalization of the Hamiltonian, the Thouless energy was normalized by the total bandwith Λ, allowing for potential comparisons with other systems.
invariant but the operator is not. The operator matrix elements in the computational basis will be complex and, since the change-of-basis matrix for going to the energy basis is a real orthogonal matrix, it will also have complex entries in the energy basis.
• GUE + complex r ab : Since the Hamiltonian already breaks time reversal, the matrix of eigenvectors expressed in coordinates over the computational basis will be a random unitary, and hence in general the operator will have complex entries in the energy basis regardless of whether it was real or complex in the computational basis (i.e. regardless of whether it is invariant under T or not, respectively.) • GUE + real r ab : Along the lines of the previous point, we shall conclude that this configuration is just an atypical case, not particularly physically meaningful for discussions regarding time reversal. We have nevertheless still kept the results for this case in Figure 8 for the sake of completeness of the analysis.
Disregarding for the discussion the seemingly unphysical GUE+real configuration, Figure 8 illustrates the fact that in systems where time-reversal is broken, either by the Hamiltonian or by the operator, depict a systematically higher K-complexity saturation value. We have also observed a continuous dependence of the saturation value on the Thouless energy that can pump up the former from the lower limiting value attained when the observable satisfies the pure RMT operator Ansatz throughout the spectrum.

ETH in the deformed XXZ
As the integrability-breaking defects studied in Section 3 are made stronger, the spectrum of the Hamiltonian of the deformed XXZ chain transitions from Poissonian statistics to Wigner-Dyson statistics. At the same time, it is possible to see that the seed operator under consideration transitions from a non-ETH regime when d is small to having and ETH structure when d attains the value that makes the spectrum of the Hamiltonian chaotic. This phenomenon was already studied in works like [49][50][51].
Here we present results on ETH checks for two extreme values of d = 0, 0.94 for the system and operator that were analyzed in Figure 4b. Figure 9 displays the result. In the integrable regime, the operator does not fulfill the ETH Ansatz because the fluctuations are not Gaussian. In the chaotic regime ( d = 0.94) the operator is seen to agree with the ETH Ansatz displaying a Thouless energy normalized by the spectral bandwidth of roughly E T Λ ∼ 0.05. In this chaotic regime, we have that the spectrum of the Hamiltonian is chaotic and that the operator fulfills the ETH Ansatz with a certain Thouless energy; since these are precisely the only two ingredients defining the systems studied in Subsection 5.2 and depicted in Figure 8, one can compare the K-complexity saturation values. The universality class at hand for our deformed XXZ is "GOE+real", and we note from Figure  8 that indeed, for a Thouless energy satisfying E T Λ ∼ 0.05 one expects a K-complexity saturation value around 0.4K, consistent with what we found in Figure 4b when d = 0.94.

Discussion
We have explored the behavior of K-complexity of a strongly coupled integrable model with an integrability breaking deformation both at the integrable point and in the chaotic phase. The purpose of this study is to delineate what kind of behavior this notion of complexity has in a chaotic system as opposed to a strongly coupled integrable one. We have found further strong evidence for the picture proposed in [5,7] namely that an exponentially large Kcomplexity saturation value at late times is a generic feature of a quantum chaotic system, while integrable systems, even strongly coupled ones which do not show exact degeneracies of energy levels, have quantitatively lower saturation values. We studied K-complexity and its late-time saturation value for XXZ systems with two types of integrability breaking terms and found that increasing the value of the coefficient of the integrability breaking term, results in an increasing value for the late-time saturation of K-complexity. We further compared the late-time saturation value in the chaotic regime to the results for RMT in the corresponding universality class, finding reasonable agreement. Along the way, we noted that non-zero operator one-point functions can influence the late-time behavior of K-complexity in a similar fashion to how the disconnected piece of the two-point function governs the late-time regime of the correlator if it is not subtracted.
We will end this discussion with a number of open questions and further avenues of research surrounding K-complexity. Firstly, it would clearly be of great interest to develop Figure 9: ETH checks for the system and operator studied in Figure 4b. On the left we show the integrable regime d = 0 while on the right we consider the chaotic regime with d = 0.94. Since we are studying operators with a zero one-point function, the ETH check can focus on just the off-diagonal elements of the operator in the energy basis. We note, in agreement with previous works (such as [49][50][51]), that the difference between the integrable and the chaotic phase is subtle: in both cases it is possible to extract a smooth envelope f (ω ab ) for the off-diagonal matrix elements of the operator as a function of the energy difference ω ab ≡ E a − E b , whose high-frequency tail can be fitted by an exponential form f (ω) = Ce −ω/E T yielding in both cases a very similar Thouless energy E T /Λ ∼ 0.055 (the r-value of the fit being around 0.992), where Λ denotes the spectral bandwidth. The difference between the integrable and the chaotic case is that in the former the fluctuations of the off-diagonal elements are not quite Gaussian, and hence they cannot be claimed to fulfill ETH, whereas they become more Gaussian in the latter, where the integrabilitybreaking defect is stronger. a more analytical understanding of the non-trivial phenomena we have uncovered in this paper, perhaps by attacking it from the angle of 'Krylov localization' as in [7], that is to analytically establish localization of the relevant part of the wave-function on the Krylov chain. A further interesting avenue to explore concerns time-dependent Hamiltonians, and how one might develop a viable generalization of Krylov complexity in such circumstances 4 . Finally, since a particularly interesting application of K-complexity is in the context of holographic duality, it is desirable to incorporate K-complexity into the holographic dictionary. In this context it is intriguing to remark that [52] recently proposed a discrete holographic dual of the aperiodic XXZ chain. It may be of interest to generalize our results to this case, to allow a more direct comparison to a simple discrete holographic code.
with unit variance (and hence the elements r ab are also of order D 0 ). Note that the operator (A.1) with the matrix elements given by (A.2) is normalized 5 according to the operator inner product The autocorrelation function is given by Due to normalization of the operator (A.3), the two-point function (A.4) starts at 1, i.e. φ 0 (0) = 1. The Ansatz (A.2) has some implications on the late-time behavior of φ 0 (t), which we can study by performing a long-time average: We can now use that, for ω = 0: With this, assuming no exact degeneracies in the energy spectrum, the long-time average of (A.4) eliminates the contribution of the off-diagonal matrix elements and yields: where in the second equality we have used the Ansatz (A.2). We thus conclude that the long-time average of the two-point function is dominated by the square of the one-point function. This fact is also in qualitative agreement with large-N factorization, i.e. in the thermodynamic limit two-point function becomes disconnected at late times. In order to probe spectral correlations we can choose to subtract explicitly the onepoint function squared from the auto-correlation function (A.4), which defines the so-called connected two-point function: where, to alleviate notational crowding, we have implicitly assumed that O is hermitian, and we have made use of the fact that the one-point function is time-independent, O(t) = O . Again, making use of the expression of the operator O in the energy basis, we can write (A.8) as: where we have used that O = 1 D Tr[O] is the infinite-temperature one-point function of the operator O. As defined in (A.8), the connected two-point function is not normalized so that its value at t = 0 is exactly one, but this is not important because we can still prove that φ (A.10) as can be seen plugging in the Ansatz (A.2). We note that the leading term in (A.10) is the first term of the last expression, consisting of a sum of D 2 numbers of order one, divided by a D 2 factor, and hence φ Plugging the Ansatz (A.2) in (A.11) we again find that the terms involving the order-one quantity O cancel out, yielding: The quantity inside the braces is of order D 0 . We thus conclude that the connected twopoint function has a long-time average of order 1 D , and that this is deduced from the ETH-like Ansatz (A.2). To argue that φ (c) 0 (t) actually plateaus at 1 D , one should prove that its long time variance is (exponentially) suppressed 6 , so that the function remains close to its long-time average at late times. We shall do that later, when studying the long time average of the square of the two-point function. But before that, we can note that there was a simpler way to derive the previous results, by defining a new operator O obtained by subtracting the one-point function from the initial operator O: (A.13) 6 We shall refer to quantities of order 1 D or smaller as exponentially suppressed because the Hilbert space dimension is typically exponential in the number of degrees of freedom S of the system, i.e. D ∼ e S .
Using the operator Ansatz for O given in (A.2), we note that the matrix elements of O in the energy basis are given by: where all r ab are of order one and the matrix R ≡ ( r ab ) is related to the matrix R ≡ (r ab ) through: In particular: from where it is apparent that all r ab are of order one and that they follow the exact constraint Tr[ R] = D a=1 r aa = 0. This operator redefinition is useful because we immediately note that the connected two-point function of O is identically equal to the full two-point function of O: (A.17) And thus, using the expression of O in the energy basis (A.14) it is straightforward to see that: from where it is immediate that φ  which for n = 0 takes the form (2.12). From the Ansatz (A.2) we find that: We thus find that the long-time-average of the square of the autocorrelation function behaves like The time averaged transition probability Q 00 takes an order-one value controlled by the onepoint function. In order to see an exponentially suppressed plateau, we again need to work with the connected two-point function φ (c) 0 (t), and the associated probability P (c) 00 (t) := φ (c) 0 (t) 2 (note that the two-point function is always real provided that the operator is hermitian, which we assume), whose long-time average we shall denote Q (c) 00 . As we showed in (A.10), φ (c) 0 (t = 0) ∼ 1, and therefore P And hence P  In a previous work on cSYK 4 [5], we studied the K-complexity of hopping operators, h ij = c † i c j + h.c. These operators have a zero one-point function, and hence their two-point function is connected. However, non-universal effects due to a non-zero one-point function can be probed if we consider, for example, an on-site number operator n i = c † i c i . Indeed, since in cSYK 4 we work in fixed occupation sectors, the one-point functions of the on-site number operators are constrained by the relation: where N is the total number operator and L is the number of sites (or rather, the number of complex fermions). Hence, (A.24) together with the fact that n i ≥ 0 implies that at least one on-site number operator needs to have a non-zero expectation value whenever N > 0. In fact, the chaotic character of cSYK 4 seems to distribute equally the expectation value accross all the n i , and given a fixed occupation sector we can thus estimate n i ∼ N L for all i = 1, ..., L, i.e. the one-point function equals the filling ratio.
This non-zero one-point function controls the averaged transition probability Q 00 , which becomes of order one in system size and, as discussed above, has the effect of lowering the K-complexity saturation value in (2.10). Figure 10 illustrates this claim.
In order to avoid this, we can seed the Lanczos algorithm with the subtracted version of the operator, O, which by construction has a zero one-point function. Following the usual arguments [3,53], we find that the Lanczos coefficients b n of O are in one-to-one correspondence with the moments µ n of the connected two-point function of O, φ (c) 0 (t), and that the K-complexity long-time average is given by: 7 Actually, in order to show that P    Vertical lines mark the estimated K-complexity saturation value. For this system size, K = 15751 ∼ 10 4 , and we observe that Q 00 ∼ 1 for the full number operator, signaling that indeed the one-point function dominates the late-time regime of the transition probability. The rest of the Q 0n with n > 1 for the full number operator seem to be rather uniformly distributed, but with a lower value due to the constraint K−1 n=0 Q 0n = 1. This eventually enforces a complexity saturation value of ∼ 0.2K, much below the naively expected ∼ K 2 , which the hopping operator does display. Conversely, both for the hopping operator and for the subtracted number operator, which have a zero one-point function by construction, all the transition probabilities are rather uniformly distributed around 1 K , yielding a K-complexity saturation value much closer to K 2 .
where Q 00 = Q (c) 00 . Therefore, as argued above, for an ETH operator this last quantity will be exponentially suppressed, hence not competing with the other uniformly distributed Q 0n with n > 0 and allowing for a K-complexity saturation value closer to K/2. This is illustrated in Figure 10.

A.2 Role of the one-point function in XXZ
This may raise some concern regarding previous work in XXZ [7], as the operators used in that case were on-site Pauli sigma matrices, whose one-point function in fixed-magnetization Hilbert space sectors are also constrained by the value of the total magnetization in the given sector, through: Figure 11: Long-time averaged transition probabilities for an instance of XXZ studied in [7], this time also considering a version of the operator where the non-vanishing one-point function has been subtracted. In contrast to the cSYK 4 case, this time the subtraction of the one-point function doesn't alter drastically the K-complexity saturation value, since in this case undersaturation is due to the monotonously decaying profile of Q 0n that we associated to Anderson localization on the Krylov chain in [7], together with the fact that even the connected part of the two-point function is itself not exponentially suppressed at late times due to the integrable nature of the system.
where, for XXZ, N denotes the number of chain sites. Since XXZ is integrable, we don't necessarily assume that all σ z n are similar, but it is anyway clear from (A.26) that in general they need not be zero. This might make one think that the XXZ calculations should be re-made taking as an input the subtracted version of on-site Pauli matrices. Figure 11 illustrates that, in this case, subtracting the one-point function doesn't alter qualitatively the results because the connected part of the two-point function in this integrable system is already not exponentially suppressed, and hence the late-time value of the two-point function doesn't change drastically if one subtracts the disconnected part from it. such that K−1 i=0 |O i | 2 = 1. We will call such an operator a 'flat' operator. We recall from [7] that for an odd number of elements in the Krylov basis, which is the case when no degeneracies are present and the operator has a non-zero projection over all Liouvillian frequencies (in such a case K = D 2 − D + 1 which is an odd number), the Liouvillian eigenvector at the middle of its spectrum has zero eigenvalue the first element is given by where we used (B.1) directly. The next element is where in the second equality we used the fact that the Krylov elements are normalized, hence 1 = i =middle |(O 1 |ω i )| 2 + |(O 1 |ω middle )| 2 and since from (B.5), (O 1 |ω middle ) = 0 we deduce that i =middle |(O 1 |ω i )| 2 = 1.
The rest of Q 0n can be computed in a similar manner, and we conclude that for n ≥ 1 One can check that this result is normalized where we used K = D 2 − D + 1 and from normalization of |ω middle ) we know that 1 D + K−1 2 n=1 X 2 n D = 1. The value of K-complexity can then be estimated as follows: where C K middle is the K-complexity of the eigenvector |ω middle ) 8 and by definition 0 ≤ C K middle ≤ K. Hence it is found that for a flat operator which for large enough D indicates that independently of the spectrum or Lanczos coefficients data. We show numerically that this is indeed the case by studying flat operators with Hamiltonians of GOE statistics and Poissonian statistics in Figure 12.