Krylov localization and suppression of complexity

Quantum complexity, suitably defined, has been suggested as an important probe of late-time dynamics of black holes, particularly in the context of AdS/CFT. A notion of quantum complexity can be effectively captured by quantifying the spread of an operator in Krylov space as a consequence of time evolution. Complexity is expected to behave differently in chaotic many-body systems, as compared to integrable ones. In this paper we investigate Krylov complexity for the case of interacting integrable models at finite size and find that complexity saturation is suppressed as compared to chaotic systems. We associate this behavior with a novel localization phenomenon on the Krylov chain by mapping the theory of complexity growth and spread to an Anderson localization hopping model with off-diagonal disorder, and find that localization is enhanced in the integrable case due to a stronger disorder in the hopping amplitudes, inducing an effective suppression of Krylov complexity. We demonstrate this behavior for an interacting integrable model, the XXZ spin chain, and show that the same behavior results from a phenomenological model that we define: this model captures the essential features of our analysis and is able to reproduce the behaviors we observe for chaotic and integrable systems via an adjustable disorder parameter.


Introduction
The notion of complexity is playing an increasingly important role in our understanding of quantum systems, from the concrete to the abstract.It is important to quantify the intuitive idea of how hard it is to simulate a given quantum system, either on a classical computer for practical purposes, or in fact on a quantum computer or a simulator in a putative laboratory setting [1].This set of questions has motivated the introduction of gate complexity [1], aiming to quantify how complicated a quantum circuit is needed to simulate a quantum system of interest.However, gate complexity appears to be hard to uniquely generalize from discrete settings to continuum ones.This is required for investigating complexity in the context of the states in the continuous Hilbert space of the world of quantum field theory.These are of additional interest if one wishes to import complexity into the domain of gravitational physics 1 .In this setting, complexity, appropriately defined, is being entertained as a fundamental quantity for our understanding of quantum gravityfor the most part via the AdS/CFT correspondence.Attempts are to concretely implement this by associating to it geometrical quantities such as the volume of the black-hole interior [6,7], or momentum [8,9], or indeed one of a multitude of other options [10].A challenge for this set of ideas is that while it is possible to quantitatively calculate the behavior of the bulk complexity, once defined, the mathematical definition of the corresponding fieldtheory boundary complexity to be equated to these geometrical quantities is not well enough defined.In particular, all the definitions suggested in [1,11] for the quantum boundary system require ad hoc add ins, be they arbitrary tolerance parameters for gate complexity or penalty functions for the Nielsen 'complexity geometry' definition.It therefore behooves us to explore the different options available, with particular regard to their possible extension to the continuum field theory setting.
At a very basic level, an appropriate notion of complexity should enable us to think of chaotic strongly correlated many-body systems, or chaotic strongly coupled quantum field theories, as being more complicated or complex than their free, or more generically interacting integrable cousins.In other words we expect chaotic theories to time evolve simple states and operators efficiently into more complex states and/or operators, while on the other hand integrable non-chaotic theories are expected not to increase complexity as fast.In this paper we will concretize this scenario for a particular notion of complexity, namely Krylov complexity, to which we will turn after a few more general remarks.
The notion of complexity has received attention in non-integrable theories [6,[12][13][14][15]] and much effort has been focused on exploring their chaotic properties which are supposed to parallel those of black holes, which are examples of maximally scrambling systems of finite entropy S in the microcanonical sense.The comparison to black holes leads one to expect an early period of exponential growth, which gives way at times of order t s ∼ log S to a linear growth phase which in turn lasts until t ∼ exp S. At this time complexity reaches a plateau where it assumes values proportional to the size of the Hilbert space -that is exponential in S. At very long time scales of order Poincaré time, e e S , we expect complexity to exhibit Poincaré recurrences (when applicable).At early times, the exponential rate of complexity growth has been linked to quantum versions of Lyapunov exponents [12,16] and their bounds 2 .
In order to put our arguments on a more concrete footing, in this paper we focus our attention on the notion of Krylov complexity [12], C K (t), (also referred to as 'K-complexity'), which depends both on the system Hamiltonian H and an operator of interest O. Intuitively C K (t) measures the spread of the operator O in operator space by characterizing its growth under Heisenberg evolution with respect to the Hamiltonian H. Specifically it does so by quantifying the average position of O(t) in a judiciously chosen basis, obtained by orthonormalizing operators obtained from successive commutations with the Hamiltonian, [H, . . ., [H, O]].The output of this orthogonalization procedure is an ordered orthonormal basis and a set of normalization coefficients, the Lanczos coefficients b n .As we shall explain further below, the spread of the operator over this basis can be understood in terms of a quantum-mechanical hopping model of a particle on an auxiliary space, the so-called Krylov chain, with hopping amplitudes determined by b n , and with K-complexity playing the role of the position expectation value on that chain.It should be noted that K-complexity is bounded from above by the Krylov space dimension, which is always smaller than the dimension of the Hilbert space of operators, through its construction.It also does not depend on a tolerance parameter as gate complexity does [1], or on a penalty metric as geometric complexity does [11].Other recent work related to the study of K-complexity includes [17][18][19][20][21][22][23].In particular, works like [17,19,22,24] have applied it to the realm of holographic conformal field theory.
In a previous paper we developed highly accurate and stable numerical methods, which in conjunction with analytical bounds, allowed us to determine the complete b−sequences and resulting Krylov dynamics up to late time for a class of strongly interacting chaotic many-body systems [14].We found that the Krylov space in SYK 4 [25][26][27] for ETH operators [28][29][30] is maximal with a dimension that scales exponentially with system size, and that the associated Lanczos coefficients b n display an initial quasi-linear growth up to n ∼ S followed by a very slow nearly linear non-perturbative descent to zero at n ∼ e S .Note that SYK 4 is rare example of a strongly interacting system, in which ETH can be shown to apply both numerically and analytically, [31][32][33].The aforementioned behavior of the Lanczos sequence for ETH operators then induces a K-complexity profile that starts growing exponentially up to the scrambling time t s ∼ log S and then transitions to a long period of linear growth up to eventual saturation at a value exponentially big in system size at times of the order of the Heisenberg time t H ∼ e S .The qualitative shape of the profile obtained was in agreement with that of the time dependent profile expected from bulk considerations [7,[34][35][36].While the emphasis of that work was put on the behavior of chaotic models, it is important to contrast their behavior with that of integrable systems.This line of inquiry was initiated in [14] by studying the SYK 2 model, and we found that its b−sequence terminates far below the bound D 2 − D + 1 (here D is the dimension of Hilbert space) which in turn leads to K-complexity saturating far below the value of the chaotic SYK 4 model.In this paper we investigate whether under-saturation of K-complexity is a ubiquitous phenomenon in more generic integrable models, allowing thus to distinguish their complexity signature from chaotic ones.The main thrust of the argument follows from studying the b−sequence and its associated one-particle hopping problem in integrable theories, where we argue that the presence of near-degeneracies resulting from the Poissonian nature of the spectrum of these systems implies an erratic structure akin to that encountered in random hopping models that show wave function localization.More precisely, we will map the Krylov chain dynamics to an Anderson random hopping model with off-diagonal disorder, rather than the diagonal disorder of the original Anderson problem [37].Our main finding is the observation of a disordered structure in the Lanczos sequence of a local operator in the XXZ spin chain, a canonical example of strongly-interacting in-tegrable system solvable via the Bethe Ansatz, which induces some degree of localization that prevents maximally efficient propagation of the particle in Krylov space, resulting in a K-complexity saturation value at late times that is lower than that obtained in previous studies for SYK 4 [14].We confirmed that this effect was due to disorder by studying a phenomenologically-inspired model with a disordered Lanczos sequence with adjustable fluctuation strength, yielding results that reproduce closely those obtained from XXZ.
This paper is organised as follows: Section 2 summarizes the notions of Krylov space and K-complexity, as well as the mapping between Krylov chain dynamics and a onedimensional hopping problem, provides a heuristic argument for the expected disorder in the Lanczos coefficients of integrable systems, and presents the framework of the Anderson problem with off-diagonal disorder, introducing some tools that we will use along the way.Section 3 presents and discusses the XXZ spin chain, which will be used for the main computations in this paper.Section 4 presents the numerical results obtained from XXZ and a comparison with SYK 4 .Section 5 shows results obtained from phenomenological models of disordered Lanczos sequences that display similar results to XXZ.Finally, our results raise some open questions in the realm of quantum many body physics, as well as some interesting puzzles within the framework of the AdS/CFT correspondence.These are discussed in Section 6.
2 Dynamics in the Krylov chain

Krylov space and the Lanczos algorithm
During its time evolution, an operator O gradually explores a subspace of the space of operators of a given system.This subspace is referred to as the Krylov space, and its structure as well as the rate at which typical operators explore it can be elucidated by constructing a suitable orthonormal basis, the Krylov basis, according to the Lanczos algorithm [12,38,39].For convenience, we shall use the smooth-ket notation |O) to denote the element of operator space corresponding to an operator O.In the Heisenberg picture, time evolution is generated by the Liouvillian operator3 L, which satisfies the relation , where H is the system Hamiltonian, so that: This motivates the definition of the Krylov space associated to the operator O, H O , as the linear span of all powers of the Liouvillian acting on it or, in other words, the linear span of all nested commutators of H with O: For a physical system with a Hilbert space of finite dimension D the operator space has dimension D 2 .Therefore, even if the set of all powers of the Liouvillian acting on O contains infinitely many elements, only a finite number of them can be linearly independent, yielding a finite dimension for the Krylov space defined in (2.2), which we denote K.This can be shown explicitly using a spectral decomposition of O in terms of the energy basis: where {E a , |E a } D a=1 are the eigenvalues and eigenstates (respectively) of H.In fact, one can show that the elements appearing inside of the sum in (2.3)  (2.4) Therefore, as discussed in [14], the Krylov dimension K is equal to the number of nonvanishing projections of O over eigenspaces of the Liouvillian, since each application of L does not mix the projections of O over each eigenspace.This implies the upper bound which arises from the fact that the zero eigenvalue ω aa = 0 (or zero phase, as the eigenvalues of L appear as phases in the spectral decomposition of the two-point function of O) is always at least D times degenerate.Given an inner product in operator space, we can build a basis for Krylov space fulfilling simultaneously the requirements of orthonormality and of retaining some notion of ordering suitable for time evolution.In this work we shall make use of the infinite-temperature inner product, or Frobenius inner product: where • is the induced norm.An orthonormal basis for the space (2.2) is constructed by performing the Lanczos algorithm, for an operator O time-evolving under a Hamiltonian H, given the inner product (2.6).This basis, in which the Liouvillian is tridiagonal, is known as the Krylov basis.The general structure of the Lanczos algorithm is as follows: 2. For n = 1: 3. For n > 1: The algorithm terminates once all independent directions in H O are exhausted, and its output consists of an orthonormal basis for H O , {|O n )} K−1 n=0 known as the Krylov basis, and a set of coefficients, {b n } K−1 n=1 , known as Lanczos coefficients.As announced, the Krylov elements are orthonormal and retain a notion of ordering suitable for time evolution, since |O n ) is a linear combination of up to n nested commutators of H with O, and therefore they are effectively probed gradually one after the other thanks to the structure of (2.1).Note that this basis can have an infinite number of elements in a Hilbert space with an unbounded spectrum, given an operator with projections over an infinite number of energy differences.However, we will be interested in long-time behaviour of operator time evolution, and consider Hilbert spaces of finite dimension D ∼ e S , where S denotes generically system size or entropy.
In the Krylov basis, the restriction of the Liouvillian over Krylov space admits a matrix representation L mn ≡ (O m |L|O n ) which takes a tridiagonal form: This feature can be exploited to connect the dynamics over Krylov space to a one-dimensional hopping problem, as we shall discuss in Section 2.2.The tridiagonal Liouvillian matrix (2.7) can be diagonalized with eigenvectors |ω i ) and associated eigenfrequencies ω i , (2.8) Note that these eigenvectors and eigenfrequencies are not all the eigenvectors and eigenfrequencies of the full Liouvillian introduced in (2.4), but the particular direction of each eigenspace (i.e. the space spanned by eigenvectors with the same eigenvalue) selected by the operator by a non-zero projection.In particular there are only K of the |ω i ) instead of D 2 which is the number of |ω ab )'s.For a precise description of the relationship between |ω i ) and |ω ab ) see [14].

Krylov space dynamics as a hopping problem
By construction, Krylov space is the minimal subspace of operator space that contains the time-evolving operator O(t) at all times.Thus, in order to keep track of this time evolution, we can decompose O(t) in terms of the Krylov elements: where φ n (t) can be viewed as wave functions over the Krylov basis or, as we shall refer to it, the Krylov chain.This terminology is justified from the fact that, as stated in (2.7), L has a tridiagonal form in the Krylov basis, which one can write more suggestively as: We note that (2.10) resembles the Hamiltonian of a one-dimensional tight-binding model on a finite chain with K sites, with zero potential energies and hopping amplitudes given by the Lanczos coefficients.The statement is that the time evolution of the operator can be mapped to such a one-dimensional quantum-mechanical problem just by constructing its Krylov basis and studying its dynamics with respect to it [12], by solving the Schrödingerlike equation for the wave functions: with the initial condition φ n (0) = δ n0 (we assume that the operator is normalized, for simplicity) and the boundary conditions b 0 = b K = 0 which ensure finiteness of the chain.We can thus understand the label n as a position on the Krylov chain, and the Krylov elements |O n ) as localized states on that chain.In general, we always begin with a localized state |O 0 ) which over time can become delocalized depending on the dynamics of the Liouvillian dictated by the hopping amplitudes given by the Lanczos coefficients b n .In [24] this perspective was used to show a relation between non-integrability and de-localization of the operator on the Krylov chain.
Two enlightening probes of the dynamics on the Krylov chain are K-complexity, C K (t) [12] and K-entropy, S K (t) [13].The former gives the average position of the propagating packet φ n (t) over the Krylov chain, (2.12) while the latter is nothing but another name for the Shannon entropy of the wave packet, and measures how spread or localized it is: (2.13) In [12] it was proposed, exploiting the relation between the asymptotics of the Lanczos coefficients and other quantities such as the two-point function, that maximally chaotic systems should display, in the thermodynamic limit, an asymptotically linear profile for the Lanczos coefficients, b n ∼ αn, which in turn would imply an exponential growth for K-complexity C K (t) ∼ e 2αt .Subsequent work [13] proposed that, for operators in finite systems satisfying the eigenstate thermalization hypothesis, ETH [30,40], the linear growth of the Lanczos sequence should transition to a constant plateau phase around n ∼ S, where S denotes entropy or system size; this implies a transition from exponential to linear growth of K-complexity at times of order of the scrambling time t s ∼ log S, and this linear growth should persist up to times of the order of the Heisenberg time t H ∼ e S , when all the directions of the Krylov space have been explored and K-complexity would therefore saturate.
Performing numerical simulations, in [14] we were able to verify these facts for the Sachdev-Ye-Kitaev model [25][26][27] with 4-site interactions, encountering the additional feature of a non-perturbative descent that corrects the plateau of the Lanczos sequence making it terminate at b K = 0, as it should.This work also observed that typical operators4 in maximally chaotic systems should saturate the Krylov dimension upper bound (2.5) due to the absence of degeneracies in the spectrum of the Hamiltonian, whose level-spacing statistics satisfy distributions of the Wigner-Dyson type and therefore feature level repulsion [41,42], and to the fact that these operators are dense (in the sense of not having any zero elements) in the energy basis.In these systems, the wave packet φ n (t) spreads in time efficiently over the Krylov chain and at times of the order of t H it tends to a uniform distribution , a value wich is exponentially big in system size because K saturates the bound (2.5), and thus K ∼ D 2 ∼ e 2S .Conversely, free or quadratic systems present in general a large number of degeneracies in the spectrum of the Hamiltonian, as well as potentially rational relations, implying degeneracies in the spectrum of the Liouvillian; this, together with abundant selection rules due to symmetries imposing the vanishing of operator matrix elements in the energy basis, reduces the Krylov dimension and therefore the upper bound for K-complexity.This was exemplified in [14] taking the case of SYK with 2-site interactions, where the Krylov space of a single Majorana scales only linearly with system size, and hence so does the K-complexity upper bound.
This paper intends to explore the remaining middle ground: there exists a large class of strongly-interacting integrable systems which feature no exact degeneracies but still exhibit Poissonian level-spacing statistics [43] and are usually solvable by the means of the Bethe Ansatz [44,45].In these systems, typical operators5 need not be sparse in the energy basis [46], thus admitting a Krylov space whose dimension can potentially be still exponential in system size.The difference between these systems and maximally chaotic ones should therefore appear at the dynamical level and, in this spirit, this work intends to explore whether the statistical properties of the Lanczos coefficients of the former systems may induce some form of localization preventing the effective diffusion of the wave packet through the Krylov chain, thus lowering the long-time value of the expectation value of its position, which is nothing but K-complexity.

Disorder in Lanczos-sequences of integrable systems
As a hopping model over the Krylov chain, the properties of the Lanczos coefficients will determine the dynamics.As we shall show below, the Lanczos coefficients depend on the Liouvillian frequencies and on the spectral decomposition of the initial operator.We will argue heuristically that the abundance of energy differences smaller than the mean level spacing ∆ in integrable systems induces an erratic behavior in their Lanczos coefficients; this is a result of the Poissonian nature of their spectrum.Conversely, chaotic systems have fewer energy differences below ∆ due to level repulsion, and we will argue that one would expect a less erratic behavior in the b-sequence of typical operators in these systems.We have observed this effect comparing the Lanczos sequences obtained from typical operators in SYK 4 [14] with those obtained from local operators in the XXZ spin chain, see subsection 2.4 and Figures 1 and 2.
The elements b n of the Lanczos-sequence can be written as ratios of Hankel determinants of the moments of the operator two-point function (see for example [47]).The moments are given by: where is the spectral decomposition of the initial operator, and it can be proved that the Lanczos coefficients satisfy the following recursion: where are the so-called Hankel determinants, with 1 for a normalized initial operator.Note that all odd moments are zero since the operator is Hermitian and the Liouvillian spectrum is symmetric around zero.In Appendix A we show that such a Hankel determinant can be expressed directly in terms of the spectral decomposition O i and the Liouvillian frequencies ω i , where the sum is over all K n possible ways of choosing n frequencies out of the K frequencies in the spectrum of the Liouvillian, and the final product is the square of a Vandermonde determinant constructed out of these n frequencies.Given Poisson statistics in the level spacing of the Hamiltonian, one expects to find in the Liouvillian frequencies which are smaller than the mean level spacing, see Section 3.With enough such frequencies, for some n every choice of n frequencies out of K will include some very small ω a − ω b in the Vandermonde determinants, reducing the value of D n−1 compared with D n−2 (this will also depend on the projections O i over the eigenspaces corresponding to those particular sets of frequencies).Specifically near the center of the Liouvillian spectrum, for integrable systems, we expect smaller frequencies and small frequency differences than for chaotic systems.Accordingly, b 2 n given by the ratio (2.15), will undergo a fluctuation, becoming either larger or smaller than the previous one, depending on whether D n−1 is in the numerator or the denominator.Conversely in chaotic systems the smaller amount of very small frequencies, would result in a less erratic b n -sequence.

Some results from Anderson localization
As argued around equation (2.10), dynamics on the Krylov chain constitute a one-dimensional hopping problem where the hopping amplitudes are given by the Lanczos coefficients.Thus, the erratic nature of the b-sequence induced by the abundance of small energy differences in integrable models can be effectively described by off-diagonal disorder in the Anderson sense.We will therefore now collect some useful results pertaining to Anderson localization that will inform our subsequent discussion.Such problems were studied as an electron localization problem (together with diagonal disorder) in [48]; as a pure off-diagonal system in [49], where it was established that the center of the band state was localized with a localization length increasing as √ K, where K is the chain size; also in [50], where localization was shown by studying the transmission coefficient; and including long-range correlations in [51].For a review, see for example [52].
For a maximal Krylov dimension (which is attained whenever no exact degeneracies are present in the Hamiltonian spectrum and the operator is fully dense in the energy basis), the Liouvillian always has a zero frequency eigenstate.This eigenstate, which is the center-of-the-band state, can be reconstructed from the Lanczos-sequence.In particular, expanding the Liouvillian eigenstates in the Krylov basis as the eigenvalue problem can be solved exactly for a given b-sequence for ω = 0.The even elements of the eigenvector are given by: and all odd elements are zero.
In the absence of long-range correlations between the hopping amplitudes, it was first proposed in [49] that in the case of disordered b's this state is localized.The localization length was shown to be proportional to √ K and inversely proportional to σ, l loc ∝ √ K/σ, where σ 2 is the variance of the logarithm of ratios of consecutive hopping amplitudes: Since b n are fluctuating, the values of x i will fluctuate around zero, and will be randomly distributed with mean 0. The variance σ 2 will depend on the details of the b n distribution.
In Figure 1 we depict the Lanczos sequences of a one-site operator for various instances of the XXZ spin chain (for more details, see Sections 3 and 4).The statistics of their fluctuations were analysed by computing (2.21) in each case 6 .The results for σ are depicted in Figure 2, where they are also compared to SYK 4 systems of similar size using data from [14].Although disorder strength is seen to decrease with the system size for both XXZ and SYK 4 , our results indicate that indeed XXZ exhibits larger disorder in its Lanczos sequences.In particular, for comparable system sizes the disorder parameter σ is bigger in XXZ than in SYK 4

Probes of late-time behavior
In this section we introduce some of the quantities used to study and observe localization in [53], and relate them to K-complexity and its long-time-average.The time-dependent coefficients φ n (t) given in (2.9) can be written in the basis of Liouvillian eigenvectors, Figure 2: Using Lanczos-sequence data for various XXZ and complex SYK 4 systems of different Krylov dimensions, the standard deviation of the logarithm of ratios of consecutive Lanczos coefficients is represented by σ on the y-axis, while the x-axis represents the Krylov dimension.The smoother the Lanczos-sequence is, the closer σ is to zero.The lines connect the average values of each set of points for a given system size and are shown to help the reader.For each system size, each point represents a single realization of SYK and a single J zz value of XXZ.In particular we present the following data.SYK: 3 realizations for L = 8 fermions, 3 realizations for L = 9 and 5 realizations for L = 10 (their Lanczos sequences can be found in [14]); XXZ: 3 different J zz coefficients for N = 9 spins, 5 for N = 10, 5 for N = 11 and 3 for N = 12.The specific values chosen for J zz for each system size can be found in Figure 1.
where we expanded the Liouvillian eigenstates in the Krylov basis |ω j ) = K−1 m=0 |O m )(O m |ω j ).From (2.9), using the orthonormality of the Krylov elements (O m |O n ) = δ mn , we recognize (2.23) We now give |φ n (t)| 2 the interpretation of a transition probability from |O 0 ) to |O n ) at time t: This expression suggests that we interpret K-complexity as the average location on the Krylov chain to which the operator will transition at each time slice t: To study the long-time behaviour of K-complexity, which is the main goal of this paper, we perform a long-time-average over |φ n (t)| 2 : By construction, the spectrum of the restriction of the Liouvillian to Krylov space has no degeneracies, and hence the phase-differences average out so that only the diagonal terms with i = j in (2.24) contribute to (2.26): (2.27) The long-time-averaged K-complexity C K can now be interpreted in two ways: 1.The long-time-averaged expectation value of the position to which |O 0 ) can "hop": (2.28) 2. The spectral average of the K-complexities of the Liouvillian eigenstates, weighted by the overlap of such eigenstates with the initial condition: where K is the K-complexity of the eigenstate of phase ω i , that is, the average position on the Krylov chain of this stationary state.Eigenstates with smaller C (i) K will be more localized towards the left of the chain and in general will have more significant overlap with the initial condition.
The transition probability, and hence also C K , are influenced by localization effects in the eigenstates.The following two limiting cases are instructive: • Full delocalization: for which the eigenstates |ω i ) are fully delocalized over |O n ), and Applying either (2.28) or (2.29) we obtain, consistently: (2.30) • Full localization: for which the eigenstates |ω i ) are fully localized over |O n ), and In this case Q 0n = δ 0n for any n and C (i) K = i for all i, and we obtain: Notice that in this extreme case C K (t) = 0 for all times because |(O n |ω i )| 2 = δ ni implies that |ω 0 ) = |O 0 ) (up to a phase), so that the initial condition is itself a stationary state and therefore it doesn't propagate.
For the case of a constant Lanczos-sequence, the eigenstates and eigenvalues of the Liouvillian as well as C K can be computed analytically.This case is shown in detail in Appendix C.
We shall now turn to discuss the Hamiltonian of the XXZ spin chain, an interacting model whose integrable nature makes it an ideal candidate for studying the erratic structure of the Lanczos coefficients in this kind of systems and the consequent imprints that localization leaves on Krylov space dynamics.
3 Complexity in integrable models: the XXZ spin chain Motivated by the arguments in Subsection 2.2, we intend to study operator dynamics in systems with an exponentially big Krylov dimension for typical operators but which still fall within the category of integrable systems.These conditions are fulfilled by stronglyinteracting systems that are integrable via the Bethe Ansatz.In fact, in the framework of the algebraic Bethe Ansatz it is possible to show that there exist extensively many commuting conserved charges that can be retrieved from the series expansion in the spectral parameter of the transfer matrix [45,54,55].However, these symmetries do not necessarily imply an extensive number of degeneracies in the spectrum of the Hamiltonian: rather, they enforce an intricate block-diagonal structure in which eigenstates carry several quantum numbers corresponding to different conserved charges (even though the underlying global symmetry might not be immediately apparent and the charges might not even be local).As a result of this, energy eigenvalues are effectively uncorrelated and therefore the nearest-neighbor level spacing statistics become Poissonian, implying that eigenvalues in the spectrum don't feature level repulsion: They can be arbitrarily close to each other, even though not exactly degenerate.Additionally, local operators in these systems can be dense in the energy basis, even though the statistical properties of their matrix elements don't follow the ETH, as exemplified in [46].These features have the effect that the Krylov space dimension for typical operators in this kind of systems is close to the upper bound (2.5) or, at least, scales exponentially in system size as it would do for an ETH operator in a chaotic system.The distinguishing features of interacting integrable systems should be found in the impact on the actual dynamical evolution along the Krylov chain of the underlying spectral statistics and of the structure of the particular operator.The XXZ spin chain [45,[56][57][58] is a canonical example for this kind of models and we shall adopt it as the main workhorse of this project.Its Hamiltonian is given by: where J is an overall dimensionful energy scale 7 and J zz is the dimensionless coupling strength.Note that the special values J zz = 0, 1 correspond to the XY and XXX models, respectively, and that the summation limits signal implicitly open boundary conditions (OBC), as opposed to periodic (PBC).Unless stated otherwise, we will use OBC in our numerical analysis.The Hamiltonian (3.1) is integrable via Bethe Ansatz.In order to unveil the Poissonian nature of its spectrum, it is necessary to first mod out the discrete symmetries that do cause exact degeneracies and restrict to a Hilbert space sector whose dimension will still be exponential in system size.H XXZ commutes with the two following global symmetries (see for example [59]): • Total spin projection along the z-direction (or total magnetization): where S z n = 1 2 σ z n denotes the spin projection along the z-axis at site n.
• Parity, defined as the transformation that exchanges sites as n ←→ N + 1 − n: where P nm is the permutation operator between sites n and m.As defined, parity is a space inversion with respect to the chain center.
There are other global symmetries that can be deduced from the algebraic Bethe Ansatz, but numerical checks indicate that only parity induces exact degeneracies when OBC are chosen.For instance, we can consider charge conjugation, which takes the form of a spin flip on each site or, up to a global phase, a rotation of angle π: It is possible to show that [H XXZ , R] = 0 = [P, R], however, [R, S z ] = 0. Since [S z , P ] = 0, we choose the set of commuting observables {H XXZ , S z , P }, with which it is possible to block-diagonalize H XXZ and restrict ourselves to a Hilbert space sector of fixed parity and total magnetization.The total Hilbert space of the system, H, has dimension 2 N and, thanks to S z -symmetry, it can be split as the direct sum of sectors with a fixed number M of spins up (i.e.fixed magnetization sectors): In each of these sectors, whose dimension is D M = N M , the total magnetization operator is proportional to the identity, S z = M − N 2 and the restriction of the Hamiltonian on each of them is closed.In future discussions, we will consider without loss of generality only sectors with M ≤ N 2 , since the M -sector is mapped to the N − M -sector by R symmetry 8 .Furthermore, the analysis based on Bethe Ansatz shows [45] that, in each sector, the spectral width scales with M for fixed interaction strength J zz , since the energy of each eigenvalue can be expressed as the sum of M magnon dispersion relations, which are bounded.Therefore, in order to have a spectrum whose width doesn't scale with system size, we choose to normalize the Hamiltonian in a fixed-magnetization sector by the number of up spins; that is, we will work with H (M ) given by: where XXZ denotes the restriction of the XXZ Hamiltonian H XXZ to a fixed magnetization sector with M spins up.The spectrum of H (M ) can present a number of exact degeneracies due to parity symmetry.In fact, each magnetization sector can be further split into positive and negative parity subsectors: whose dimensions D ± M are computed in Appendix B and summarized in Table 1.Importantly, they scale asymptotically exponentially in system size.As numerical checks indicate, once we restrict ourselves to a sector H P M , the spectrum of the Hamiltonian doesn't feature exact degeneracies due to any other discrete symmetries, and we can extract the universal Poissonian behavior of the spectral statistics, as ilustrated in Figure 3, where the level spacing distribution of the positive-parity sector in XXZ with both OBC and PBC are compared to that of an instance of complex SYK 4 with the same Hilbert space dimension.This serves to illustrate a point that could be misunderstood: the integrable character of the system reflects itself in the lack of correlation between the energy eigenvalues, yielding a Poissonian distribution characterized by a high probability of finding energy differences smaller than the mean level spacing, as opposed to the characteristic level repulsion displayed by chaotic systems.Exact degeneracies are not necessarily a defining feature of integrable systems.
As a word of caution, we shall note that, for chains with an even number of sites N , there exists a zero-magnetization sector (M = N 2 =⇒ S z = 0) in which accidentally S z and R do commute, since R doesn't change the magnetization sector because there are as many spins up as spins down.Numerically (see for instance Figure 3) we don't observe extra degeneracies produced by this symmetry enhancement.In fact, it is conjectured in the literature [60] that eigenstates within this magnetization sector are themselves eigenstates of R. In conclusion, for chains of even length, in the Hilbert space sector with M = N 2 and fixed P we can classify eigenstates according to their R-eigenvalue (±1), and the parity sectors can be split into R-subsectors: for N even (3.8) Figure 3: Nearest-neighbor level spacing distributions for XXZ with a specific J zz coupling with N = 14 spins in the sector M = 7, P = +1, both with open and periodic boundary conditions.For comparison, we include the corresponding distributions for complex SYK 4 with L (SY K) = 13 sites and occupation number N (SY K) = 6, whose Hilbert space dimension is equal to that of the positive parity sector of the above-mentioned XXZ instance.All distributions are normalized so that the area under each curve is equal to one, and all energy differences have been normalized by the mean level spacing ∆ in each system.
whose dimensions are also computed in Appendix B and can be written compactly as: for N even. (3.9) Despite not granting exact degeneracies in the spectrum of the Hamiltonian, this symmetry enhancement can yield some selection rules for operators charged under R, which impose constraints on their matrix elements in the energy basis and hence can reduce the dimension of their associated Krylov space, as also discussed in Appendix B.

Numerical results
In this section we present evidence that in the case of H XXZ discussed in the previous section and for a single-site operator which we introduce below, K-complexity indeed saturates at values below ∼ K/2, which is the expected saturation value for chaotic systems such as SYK 4 studied previously in [14], where the Lanczos sequences featured less disorder as illustrated in Figure 2. We will work with H (M ) defined in (3.6) in the sector of fixed magnetization and parity H + M and will study the time evolution of the linear operator M odd Table 1: The dimension of each Hilbert space sector of fixed parity and total magnetization is , where D M = N M and the expression for A, which depends on whether N and M are even or odd, is shown in this Table .where m can be any site on the spin chain 1 ≤ m ≤ N .This operator was chosen because it respects the sector symmetries, namely total magnetization conservation and parity, in addition to being a local one-site observable.We choose m to be near the center of the spin chain, so that boundary effects will not play a role.In general, all the matrix elements of this operator in the energy basis (2.3) are non-zero and therefore, given that we work in sectors with no degeneracies, the Krylov space dimension for such an operator is expected to saturate the upper bound of (2.5).We tested this assertion by studying the matrix elements of this operator in the energy basis to very high precision.It would be interesting to also study other local operators; preliminary checks performed at small system sizes for linear combinations of one-site operators, two-site operators and linear combinations of two-site operators display a similar phenomenology to what will be described in this section for the operator (4.1).Some of the Hamiltonians and operators were generated using the open-source package [61], while others were produced with programs of our own implementation, yielding coincident results.
Using high precision arithmetic and the re-orthogonalization algorithms we developed for the study of SYK and which are described in [14], we computed the Lanczos coefficients associated to operator (4.1) for various instances of XXZ, see Figure 1.Then, following the discussion in Section 2.5, we plot in Figure 4 the time-averaged transition probabilities Q 0n for several XXZ systems with different numbers of spins N , different magnetization sectors M and various choices of the J zz coupling.The weighted average of Q 0n , which is equivalent to the long-time average of K-complexity via equation (2.28), is plotted as a vertical line for each system.In each case, we plotted also the transition probability and C K for a similar sized complex SYK 4 system studied in [14].To make the results clearer we normalized the x-axis with respect to the Krylov dimension K for each system separately, as well as normalizing the transition probability in the y-axis with respect to K −1 .We find that for XXZ the transition probability is biased towards the left side of the chain, whereas for SYK it is almost horizontal and of the order of K −1 .The left-biased transition probability for XXZ moves the saturation value of C K to a value smaller than K/2, while the almost flat profile of Q 0n for SYK gives C K a value close to K/2.Equation (2.29) expresses the long-time average of K-complexity, C K , in terms the overlap of the Liouvillian eigenstates with the initial condition |(ω i |O 0 )| 2 and the spectral complexities C (i) K .The corresponding quantities computed numerically from XXZ and SYK are presented in Figures 5 and 6.We find that two factors cooperate to decrease the value of K-complexity in XXZ compared with SYK: (1) In XXZ the spectral profile of the initial operator is peaked towards the band center9 , while for SYK the profile is flatter and of order K −1 as expected from ETH [40].(2) In XXZ, K-complexities of the individual eigenstates near the band center reach smaller values (Figure 6) due to enhanced localization, compared with those of SYK where most eigenstates have K-complexities of around K/2.These eigenstates are particularly important, since the initial condition is localized to the left of the chain and therefore will have a larger overlap with them.In general, the eigenstates in XXZ which have a stronger overlap with the initial condition, K 's than those of of SYK, thus reducing the value of C K computed via (2.29).
Time-dependent results of K-complexity, K-entropy and snapshots of the absolute value of the wavefunction at various time scales are displayed in Figures 7, 8 and 9. Kcomplexity is seen to saturate below K/2 as predicted by its late-time average.The time scale at which it saturates is of order K. Since K ∼ e 2S we find that the saturation time for XXZ is similar to that observed in chaotic models such as SYK.K-entropy saturates as well at time scales of order K, at a value smaller than log(K) ∼ N -the value signaling total delocalization of the operator wave function as seen in SYK in [14] -which is another indication that the wave function is more localized than in SYK.Again, the saturation time is of order e 2S as seen in SYK.The snapshots of the wave function at increasing time scales show that, although the wave packet is spreading towards the right of the chain, there is always a stronger support on the initial site of the Krylov chain, and in general some bias towards the left region.
Note: In the case of an even number of spins in the zero magnetization sector M = N/2, both P and R preserve total magnetization and in Appendix B.2.1 we argue that O of the form (4.1) is charged under R, which creates a selection rule and reduces the Krylov space dimension.Numerical results and discussion for N = 10, M = 5 are given at the end of Appendix B.2.1.

Phenomenology of erratic Lanczos sequences
We have computed Lanczos sequences for various instances of XXZ and verified, in agreement with our initial expectation, that they feature disorder superimposed on the ascentdescent profile.Such Lanczos sequences play the role of disordered hopping amplitudes in the Krylov chain, and we observed that for those systems K-complexity saturated at late times at values smaller than half the Krylov space dimension, K 2 , as a consequence of the long-time averaged transition probability being biased towards the left side of the chain (i.e. the side closer to the localized initial condition).In the framework of Anderson localization due to off-diagonal disorder, it is usually difficult to establish analytically features of    time evolution, and customarily one focuses on properties of stationary states [49,50,52].For this reason, in this section we will present numerical computations with heuristically motivated Lanczos sequences whose main ingredients will be an ascent-descent profile with fluctuations on top, and we will show that they reproduce very closely the XXZ results, backing up the claim that the K-complexity phenomenology observed in this model is due to the disorder of its Lanczos coefficients.We will construct such Lanczos sequences for Krylov chains of size K in successive levels of sophistication: First and for reference, we compute transition probabilities and K-complexity resulting from a flat b-sequence with fluctuations on top, which is nothing but the canonical setup for off-diagonal disorder; then, we will use two different Ansätze for a more realistic underlying profile displaying ascent and descent regimes.

Pure off-diagonal disorder
In this canonical model for off-diagonal disorder, we generate a random sequence of Lanczos coefficients (hopping amplitudes) {b n } K−1 n=1 taking all the coefficients to be independent and identically distributed (i.i.d.) Gaussian random variables with mean J and standard deviation W J, where J is to be thought of as a dimensionful parameter setting the energy scale 10 , and W is a dimensionless parameter controlling the disorder strength.Figure 10 depicts various random realizations of Lanczos sequences with different disorder strengths, as well as the corresponding long-time averaged transition probabilities and K-complexities as a function of time computed out of them.We observe that transition probabilities are more and more biased to the left as disorder increases, and in fact decay roughly exponentially, yielding a K-complexity profile that saturates at late times at values that decrease with disorder strength.We note that, in the disorder-free case, K-complexity shows per- sistent oscillations before saturating at K 2 , consistent with the wave packet bouncing back and forth between the edges of the Krylov chain before diffusing efficiently into a uniform distribution; in Appendix C we show analytically that the time-averaged transition probability is indeed a constant (up to edge effects), yielding ∼ K 2 as the K-complexity long-time average.The addition of disorder contributes efficiently to destroying the coherence of the propagating packet, thus washing out complexity oscillations apart from reducing its saturation value (c.f. Figure 10 -bottom right).

Disordered sequence with ascent and linear descent
This time we add i.i.d.Gaussian fluctuations on top of a Lanczos sequence b n that increases linearly up to b * at n * ∼ log K and then decays linearly to (almost) zero at n = K − 1.That is: where, as announced, W n are i.i.d.Gaussian random variables with zero mean and standard deviation given by the dimensionless parameter W . Sequences generated with different disorder strengths as well as the transition probabilities and K-complexities computed out of them are depicted in Figure 11.The numerical values of b 0 and b * were chosen so that the resulting profiles reassembled those of XXZ 11 .
In the absence of disorder, we observe that K-complexity saturates above K 2 .This is due to the asymmetric shape of the Lanczos sequence: Since the hopping amplitudes are smaller on the right side of the chain, the packet tends to spend more time in that region on average, as reflected in the transition probability Q 0n being slightly biased towards the right.However, as soon as disorder is turned on, the late-time saturation value of complexity drops below K 2 to some finite fraction of K.This localization phenomenon is also reflected in the fact that the transition probability decays with n at a seemingly exponential rate (even though with a small exponent), for the disorder strengths explored.
One might naively think that under-saturation of K-complexity in the XXZ results is due to the descent in the Lanczos sequence.The model presented here contradicts that statement: We have observed that, in complete absence of disorder, the ascent-descent profile results, conversely, in "over-saturation" (in the sense that the late-time average is bigger than half the Krylov dimension K 2 ).It is thanks to the addition of disorder that the late-time complexity value can drop down, as we have explicitly demonstrated (c.f. Figure 11).
There is one slightly unsatisfactory feature of this phenomenological model: The profile of the descent on top of which disorder is added is linearly decaying, but the variance of the fluctuations doesn't depend on n.Hence, the fluctuations become stronger with respect to their mean value as n increases, amplifying the decay of the transition probability.Furthermore, some coefficients towards the end of the sequence can incidentally become very close to zero, effectively "breaking" the chain.This all produces a big drop in the transition probability, as can be observed in Figure 11 (bottom left).Localization is, however, still operating throughout the chain, as Q 0n decreases with n even on the left side of the chain when disorder is sufficiently strong.A further refinement of the model, to be presented next, tries to get rid of the above-mentioned edge effect by adding some convexity towards the end of the descent.

Disordered sequence with ascent and quasi-linear convex descent
If the descent of the Lanczos sequence had a slightly convex profile, the region where the strength of the fluctuations becomes comparable to the mean value on top of which they are added would be pushed towards the right.In fact, the analytical expression for the Lanczos sequence in RMT at large size has this feature [19]: the quasi-linear descent gets eventually modified by a square-root behavior.With this motivation, we use the following Ansatz for the toy Lanczos sequence: For n * < n K the descent is linear, but towards the edge n ∼ K the square-root behavior becomes dominant.
Figure 12 depicts the generated b-sequences as well as Q 0n and C K (t).Their striking similarity to the XXZ results presented earlier in this paper backs up the claim that disorder is responsible for the obtained transition probabilities and the lower saturation value of Kcomplexity.A weak disorder W = 0.1 is able to reduce the complexity saturation value from above K 2 to ∼ K 4 .It is, again, interesting to note that the total absence of disorder would imply long standing coherent oscillations of complexity and a saturation value slightly above K 2 due to the asymmetric profile of the b-sequence that has already been discussed.This feature should remain there at larger system sizes, even though less pronounced.Hence, the fact that in a single random realization of SYK K-complexity does saturate around K 2 (in fact, even a bit below) and doesn't feature strong oscillations (c.f.[14]) indicates that We note the similarity to the results depicted in Figure 4. Bottom right: K-complexity as a function of time, C K (t), for each disorder strength.Large oscillations disappear when disordered is turned on, and the complexity saturation value decreases as disorder strength increases.even in this system the small disorder of its Lanczos sequence is operating.From this perspective, the difference between SYK and XXZ is more quantitative than qualitative: The b-sequence of the latter is more disordered than that of the former, but they both display localization in the Krylov chain to some extent.This is consistent with usual discussions in Anderson localization: in one spatial dimension there is no critical disorder strength, and this phenomenon is always present for any value of the variance of the fluctuations, but the localization length can vary drastically.

Conclusions
In this work we studied the behavior of K-complexity for local operators in the XXZ spin chain, a strongly-interacting integrable system solvable via the Bethe Ansatz.We found that the associated Krylov space dimension is exponential in system size, however the disorder in the sequence of Lanczos coefficients induces a localization effect that prevents maximally efficient exploration of the Krylov chain by the operator wave function, resulting in a saturation value of K-complexity at late times below K 2 .It should be noted that the starting Hamiltonian, i.e. an XXZ spin chain with constant coupling strength, was not defined through a set of disordered couplings: disorder appears in the Lanczos sequence and it is enhanced by the Poissonian statistics of its spectrum; localization operates therefore in Krylov space rather than in direct space.Results also indicate that the difference between XXZ and SYK is more quantitative than qualitative: in both cases there is some disorder in the Lanczos sequence, but it is stronger in the integrable model; since the dynamics on the Krylov chain are described by a one-dimensional hopping problem, where there is no critical value for the disorder strength, localization should be operating in both cases, to a different extent in each of them.In turn, the enhanced disorder in the Lanczos sequences for XXZ causes a 'tilt' in the transition probability from the initial operator |O 0 ) to |O n ) as can be seen both in our numerical results for XXZ (Figure 4) as well as in our phenomenological models, thus reducing the value of the long-time average of Kcomplexity.This is also manifest in the presence of Liouvillian eigenstates with smaller K-complexities in XXZ, especially around the band center (Figures 5 and 6), which is again due to stronger localization as a result of the larger disorder in the Lanczos sequences.
More work needs to be done in order to determine how generic this phenomenon is.In particular, it is necessary to formalize more strictly the relation between Poissonian statistics (and hence, integrability) and the erratic structure of the Lanczos coefficients, in a way that allows to predict the strength of the fluctuations of the b-sequence from data of the initial Hamiltonian (and perhaps the operator).Eventually, it would also be enlightening to push this towards an analytical estimate for the saturation value of Kcomplexity, which so far can only be accessed numerically.It should be noted, however, that this last point translates to a generically complicated problem in the conventional framework of Anderson localization.
Last but not least let us comment on the relationship between our findings and what role complexity may play in holography, the very question that has been the source of motivation to study the phenomenology of Krylov complexity in detail in this work.As we mentioned in the introduction it is of course of great interest to push ahead in the direction of identifying precise notions of complexity on either side of the duality, in an attempt to solidify the corresponding dictionary entries.It was our goal to add to the detailed understanding of K-complexity and its associated physics, as has been summarized above.Let us now venture into more speculative territory, by addressing a few puzzles and potential resolutions raised by our results.
Going beyond the necessity to juxtapose integrable and chaotic behavior in order to understand what properties of complexity can be associated unambiguously to chaotic operator spreading, and which ones cannot, general holographic wisdom would suggest that we are only interested in the behavior of K-complexity in chaotic systems.This is in accordance with the general expectation that semi-classical black holes are thermodynamic systems and the dynamics of the theory is (maximally) chaotic, as evidenced, for example by the maximal value of the OTOC-Lyapunov exponent [62], which is determined by the red-shift factor of the classical black hole.On the other hand, there is ample evidence that the field-theory side of the duality (e.g.N = 4 SYM in d = 3 + 1) is quantum integrable -at least in the planar sector, and for certain classes of observables [63].Accepting, for argument's sake, that this integrability extends to the quantities of interest in this paper, this raises a conundrum: how can both of these statements be correct at the same time?One potential route out of this impasse lies in timescales.Our methods and results have targeted very large timescales, exponential in the entropy, where we have found quantitative and qualitative differences between integrable and chaotic behavior.The semiclassically chaotic behavior of the bulk black hole is manifest at much earlier timescales, around the scrambling time, so that our results (leaving aside the obvious fact that we are talking about the XXZ spin chain, and not N = 4 SYM) are not in immediate conflict with the early growth of complexity corresponding to the black hole geometry.It would therefore be of great interest to study further the issue of earlier (∼ scrambling-) time complexity dynamics in order to see whether any significant changes can be discerned between chaotic and integrable cases.In this regard it is interesting to remark that [17] presents evidence that the early growth phase of K-complexity is not in fact a reliable indicator of quantum chaos.On the bulk side it would then be highly interesting to elucidate what other geometries, besides the black holes contribute at later time scales [64][65][66], and whether those geometries can account for the difference between chaotic and integrable cases observed in this work.Recent work [67] has investigated a bulk proposal of the late-time non-perturbative behavior of complexity in chaotic two-dimensional gravitational systems.It would be very interesting to see if a similar analysis can be made in integrable cases and whether our suppressed complexity saturation could be seen from the bulk.
Note: during the final stages of preparation of this paper, the related article [68] appeared on arXiv.It studies localization on the Krylov chain of disordered systems featuring many-mody localization, in the thermodynamic limit.They find that, in the MBL phase, the operator wave function is localized on the Krylov chain and that the growth of K-complexity is suppressed.While they studied the thermodynamic limit of systems with disordered coupling strengths, we have considered late-time dynamics of systems with a fixed, spatially uniform coupling strength, where the Lanczos coefficients appear to be disordered due solely to their integrable nature.We thank the authors of [68] for drawing our attention to their work.

A Moments and Hankel determinants
In this appendix we show how to arrive at (2.17) starting with (2.16).Firstly, let us look at Given (2.14), the matrix in the determinant can be written as where we have assumed the operator is normalized In this form the determinant is immediate where the final product is the square of the Vandermonde determinant of the Liouvillian frequencies Hankel determinants D n with n < K − 1 are a bit more complicated since the Vandermonde matrices in this case are not square.The matrix in the determinant of D n , is the product of matrices with sizes (n + 1 × K) and (K × n + 1), namely The determinant of a matrix which is a product of two rectangular matrices, can be reduced via the Cauchy-Binet formula (see for example [69]) to the form (2.17).

B Hilbert space dimension of XXZ sectors
This Appendix presents the details of the computation of the dimensions of Hilbert space sectors of XXZ with fixed magnetization, parity and R-charge (whenever the latter is also a conserved quantum number within a given sector).

B.1 Parity sectors at fixed magnetization
We start by considering the sector H M of fixed magnetization with M spins up, whose dimension is given by: and we shall proceed to split it in sectors of fixed parity, H M = H + M ⊕ H − M , whose dimensions we denote by D + M and D − M , respectively.In order to study these parity sectors, it is convenient to use the natural tensor product basis of the spins in the XXZ chain.For the sake of language economics, we shall refer to it as the "computational basis", identifying ones with spins up and zeroes with spins down.If we examine individually each element |ψ of this basis, we find the two mutually exclusive possibilities: a) P |ψ = |ψ , i.e. |ψ is already invariant under parity, hence it already belongs to H + M .b) ψ| P |ψ = 0, i.e. |ψ is not invariant, so the action of parity gives a different element in the computational basis, and therefore their overlap vanishes due to orthogonality.For the states in this class, it immediately follows that |ψ + P |ψ has positive parity while |ψ − P |ψ has negative parity.
We denote A, B as the number of states in the computational basis that satisfy (a), (b), respectively.B needs to be always even because if |ψ belongs to (b) then so does P |ψ ; the fact that P 2 = 1 groups the states in the class (b) pairwise and would lead to a contradiction if B is odd.We note that the H + M sector receives A states from class (a) and B 2 states from class (b), where the factor of 12 comes from the fact that half of the independent linear combinations of pairs constructed in (b) have positive parity and half have negative parity.To summarize: Since (a) and (b) are mutually exclusive, we have that: and therefore: So all that's left is to determine A, i.e. the number of states in the computational basis that are invariant under parity.Since parity performs a mirror transformation with respect to the chain center, the general philosophy is to count the number of states in (a) by finding the number of ways to arrange half of the ones in the left half of the chain, as the requirement of invariance under parity will determine the position of the remaining ones in the other half.However, depending on whether N and M are even or odd, some subtleties need to be taken into account.

B.2.1 Selection rules and Krylov dimension for operators charged under R
Because of the R-symmetry enhancement, operators charged under R are constrained by selection rules that enforce nullity of a part of their matrix elements in the energy basis, which has the effect of reducing the dimension of their associated Krylov space.The reason why these selection rules exist is that eigenstates of the Hamiltonian are themselves R-eigenstates, rather than being degenerate, as proposed in [60] and as we have confirmed performing numerical checks.We shall now prove explicitly those selection rules and give the expression for the Krylov dimension that they imply in each case.
• Negatively charged operator.We consider an observable O acting on the sector H P N 2 that verifies: Since eigenstates in the sector of zero magnetization and fixed parity P of this chain of even length are also eigenstates of R with eigenvalue ±1, it follows that the action of O on these states has the effect of changing the sign of the R-eigenvalue.To show this, we shall label eigenstates in the sector explicitly with their energy and R-eigenvalues, |E, r with r = ±1.Using (B. for the sake of notational simplicity).To give the upper bound for the Krylov dimension in this case, let's assume that all matrix elements that are not forced to vanish by the selection rule are non-zero.Then, the number of non-zero matrix elements E , r |O|E, r (and hence the Krylov space dimension, as there are no degeneracies in the spectrum of the Hamiltonian in this sector) is given by: where, again for notational simplicity, R ± ≡ D P,R=± .
• R-invariant operator.In this case, we consider an observable O acting on the sector H P .It is apparent from this perspective that (B.25) is nothing but the combination of the upper bounds of the form (2.5) applied separately to each R-sector, where the operator is dense in the energy basis.
Consider the operator we study, O = σ z m + σ z N −m+1 in an XXZ system with an even number of spins N in the zero magnetization sector M = N/2 and P = +1.As discussed in Section 3, in this sector both P and R are symmetries of the Hamiltonian.It can be shown, using the property σ j σ k = δ jk + iε jkl σ l of the Pauli matrices at a fixed site, that the operator O chosen here satisfies (B.18), thus being negatively charged under R. We therefore expect it to span (in the absence of other degeneracies in the spectrum) a Krylov space of dimension (B.21).We confirm this expectation for the case of XXZ with N = 10, M = 5 in the P = +1 sector with the operator O = σ z 5 + σ z 6 , where we find the Krylov space dimension to be K = 7810, as indeed predicted by (B.21), once we work with high enough accuracy.In Figure 13 we show the Lanczos sequences as well as matrix plots of the operator in the energy basis, for two different J zz coefficients.

C Analytical results for a constant b-sequence
To contrast with the disordered systems we study in the main text, we discuss here a flat, smooth b-sequence with no disorder.Such systems were considered in the context of Anderson localization, for example in [53].
Consider the eigenvalue problem (2.19) with b n = 1 (in some units) for all n, i.e.where we want to find ψ n and ω which solve this equation.Assuming boundary conditions ψ −1 = ψ K = 0, it can be solved via the transfer matrix recurrence equation: The solution in terms of the transfer matrix and boundary conditions is Finding the eigenvalues and eigenvectors of the transfer matrix, it can be written as where λ 1,2 = 1 2 ω ± √ ω 2 − 4 .This allows us to write the solution for any n in terms of ψ 0 and λ 1,2 , as follows: Now we can use the other boundary condition, namely ψ K = 0: It is also clear now from (C.5), that the eigenvector element ψ ni for eigenvalue ω i is: where we have absorbed the denominator in the normalization constant because it doesn't depend on n (recall that the label i designates the eigenstate).C i can now be fixed from normalization, and the normalized eigenstates are from which we can compute K-complexity for the pure no-disorder system: (C.12) The transition probability (2.27) is 12 and the time-averaged K-complexity is given by, This analytical result shows that for a flat and smooth Lanczos-sequence K-complexity will saturate at K/2.
12 To arrive at this result we used N n=0 e .See [53].
are eigenstates of the Liouvillian: Denoting ω ab ≡ E a − E b and |ω ab ) ≡ |E a E b |, we see that L|ω ab ) = ω ab |ω ab ) .

Figure 1 :
Figure 1: The Lanczos sequences for the various XXZ cases studied in this paper.For the details of the XXZ model see Section 3.For details on the one-site operator used and further numerical results, see Section 4.

Figure 5 :
Figure 5: The spectral decomposition of the initial operator (top) and K-complexities of the individual Liouvillian eigenstates (bottom).Data for 4 random realization of cSYK 4 with 10 fermions are superimposed in the leftmost column.Data for 3 choices of J zz coupling for XXZ with N = 12, M = 4, P = +1 are superimposed in the second column from the left.Data for 3 different choices of J zz couplings for XXZ with N = 11, M = 5, P = +1 are shown in the third column from the left; and data for 5 different choices of J zz couplings for XXZ with N = 10, M = 4, P = +1 are shown superimposed in the right column.The frequencies are normalized according to ∆ which is the mean level spacing computed for for each system separately.The dashed red line represents the value of K −1 for each system.

Figure 6 : 2 XXZ with N = 12 M = 4 = z 6 + z 7 J
Figure 6: Same results as in Figure 5 focused at the band-center.Note that the values of C (i) K for XXZ reach smaller values than those for SYK.

Figure 7 : 7 J
Figure 7: Time-dependent results for K-complexity.Shown here are results for XXZ with N = 12 in the M = 4, P = +1 sector at 3 different J zz couplings.Each subplot shows a different time scale.K-complexity is seen to saturate below K/2 at time scales of order K (up to some dimensionful prefactor).

Figure 8 :
Figure 8: Time-dependent results for K-entropy.Results for XXZ with N = 12 in the M = 4, P = +1 sector at 3 different J zz couplings are shown.K-entropy is seen to saturate below log(K) ∼ N at time scales of order K.The general lower value of K-entropy is also a sign of localization.

Figure 9 :
Figure 9: Snapshots of the absolute value of the wave-function at selected times for XXZ with N = 12 in the M = 4, P = +1 sector at 3 different J zz .At t = 0 the wave-function is a delta function on the first site on the Krylov chain.At late times the wave function shows clear signs of localization on the first site (i.e. the initial condition), and is in general biased towards the left side of the Krylov chain.

Figure 10 :
Figure 10: Top: Lanczos sequences generated by drawing each Lanczos coefficient from a normal distribution with unit mean and standard deviation W , following Subsection 5.1.Each color corresponds to a single random realization for a fixed disorder strength, as indicated in the legend.In all cases we studied a Krylov chain of length K = 2000.Bottom left: Transition probabilities computed out of each Lanczos sequence.Bottom right: K-complexity as a function of time.

Figure 11 :
Figure 11: Top left: b-sequence for the phenomenological model (5.1).One random realization for different values of the disorder strength is depicted.Top right: Same sequence, plotted with logarithmic scale along the horizontal axis.Bottom left: Timeaveraged transition probabilities Q 0n computed out of each sequence.As disorder increases they start to follow a profile that decreases with n, signaling localization.Bottom right: K-complexity as a function of time for each disorder strength.We observe that the late-time saturation value decreases monotonously with W .

Figure 12 :
Figure 12: Top left: Lanczos sequences with a slightly convex shape towards the end of the descent, with fluctuations controlled by W . Top right: Idem, with logarithmic scale along the horizontal axis.Bottom left: Transition probabilities Q 0n computed out of each b-sequence.We note the similarity to the results depicted in Figure 4. Bottom right: K-complexity as a function of time, C K (t), for each disorder strength.Large oscillations disappear when disordered is turned on, and the complexity saturation value decreases as disorder strength increases.
18) it is possible to show that: RO |E, r = −r O |E, r .(B.19)That is: While |E, r is an eigenstate of R with eigenvalue r, O |E, r is an eigenstate of R with eigenvalue −r.This in turn implies the selection rule: E , r |O|E, r = 0 if r = r .(B.20)But (B.20) are the matrix elements of the observable in the energy basis, so they are direcly related to the Krylov dimension, and the fact that many of them vanish identically because of the selection rule has the effect of lowering the Krylov dimension with respect to its upper bound D 2 − D + 1 (where we use D instead of D P N 2

Figure 13 :
Figure 13: For XXZ with N = 10 in the sector M = 5, P = +1 for which the Hilbert space dimension is D = 126 according to Table 1, the operator O = σ z 5 + σ z 6 is not dense in the energy basis.The left panel shows the Lanczos sequences which terminate at K = 7810, and the middle and right panels show matrix plots of the operator in the energy basis for the two J zz couplings for which the Lanczos sequences were computed.