Operator complexity: a journey to the edge of Krylov space

Heisenberg time evolution under a chaotic many-body Hamiltonian $H$ transforms an initially simple operator into an increasingly complex one, as it spreads over Hilbert space. Krylov complexity, or `K-complexity', quantifies this growth with respect to a special basis, generated by $H$ by successive nested commutators with the operator. In this work we study the evolution of K-complexity in finite-entropy systems for time scales greater than the scrambling time $t_s>\log (S)$. We prove rigorous bounds on K-complexity as well as the associated Lanczos sequence and, using refined parallelized algorithms, we undertake a detailed numerical study of these quantities in the SYK$_4$ model, which is maximally chaotic, and compare the results with the SYK$_2$ model, which is integrable. While the former saturates the bound, the latter stays exponentially below it. We discuss to what extent this is a generic feature distinguishing between chaotic vs. integrable systems.


Introduction
The concepts of state and operator complexity lie at the intersection of quantum information, condensed matter physics, quantum field theory (QFT) and black hole physics. One is thus faced with a diversity of concepts and methods on how to precisely define and actually measure complexity. The identification of very large time scales in black hole physics and corresponding phenomena in QFT (via the AdS/CFT correspondence) [1][2][3][4][5][6][7][8] has motivated the search for a form of complexity whose time dependence has the following distinctive characteristics: In fast-scrambling systems with finite entropy S, complexity should grow exponentially in time until t s ∼ log(S) (known as the scrambling time) when it reaches a value of order O(S); it then switches to a linear-in-time behavior until times of order exp(O(S)) when it settles around a plateau of order exp(O(S)); finally, after times of order exp(exp(O(S))) it is expected to start performing the Poincaré dance. The growth of the Einstein-Rosen bridge was connected in [4] with complexity growth of the quantum state in the dual CFT. A notion of complexity known as Krylov complexity, or K-complexity, was introduced in [9] as a useful diagnostic of quantum chaos in the thermodynamic limit. In [10] it was argued for the first time that K-complexity of generic operators in finite-size systems exhibits the profile of complexity described above. This paper reports on a complete numerical analysis of K-complexity for all time scales, including scrambling and well beyond, in concrete models demonstrating the above-described time-dependence. This is the first time that such large time scales are studied with numerical methods, and in concrete models. We uncover non-perturbative effects which were not initially anticipated. The study overcame known numerical instabilities and employed powerful algorithms executed on large computing clusters.
Roughly speaking the time evolution of complexity quantifies how quickly a reference operator grows under Heisenberg evolution in a pre-defined basis. K-complexity, denoted C K , measures this growth with respect to the Krylov basis which, as we shall see below, is well-adapted to capture Heisenberg evolution efficiently. K-complexity depends on the Hamiltonian of the system and a chosen reference operator. It has an advantage over other notions of complexity: its definition doesn't require an arbitrary tolerance parameter. For circuit complexity 1 , to cite but one example, such a tolerance parameter must be introduced, and its presence is crucial in establishing its boundedness [13][14][15], whereas K-complexity is naturally bounded from above 2 [10]. It is known to be bounded by the dimension of the Hilbert space of operators; below we prove a stronger bound.
C K provides a fine-grained notion of complexity, as it manifestly distinguishes all linearly independent operators up to a given, fixed length. These features make K-complexity a natural choice for quantifying the evolution of operators at late times (see e.g. Sec. 2.1 in [17] for a careful definition of chaotic time scales) which is this paper's focus.
In studying K-complexity at these time scales one should either explicitly construct, or obtain bounds on, the full sequence of Lanczos coefficients, b n , which characterize the Krylov basis of an operator. These coefficients arise in the process of orthonormalizing 3 the Krylov basis [18,19]. In [9] it was hypothesized that in a chaotic quantum system in the thermodynamic limit, the Lanczos coefficients grow at a linear rate. In [10] it was suggested that, in finite chaotic systems, the Lanczos-sequence could exhibit plateau-like features over a certain range. The study of finite chaotic systems reported in this paper constructs the full Lanczos sequence to reveal a non-perturbative decay of this plateau, making an initially small but eventually persistent correction to the "cliff-ended plateau" conjectured in [10]. We refer to the resulting profile as "the Descent".
The structure of this paper is as follows: In Section 2 we describe the Krylov space and an analytical upper bound on its dimension. In Section 3 we discuss the Krylov 1 Circuit complexity of the SYK4 model, the main work-horse of this paper, has previously been studied in [11,12].
2 Operator size complexity [16] shares this advantage of not requiring a tolerance parameter, but differs from CK in other crucial respects. For example, it is directly bounded by the total system size, while the upper bound for K-complexity scales exponentially with system size. 3 The process of orthonormalization requires a choice of inner product. In this work we use the standard Frobenius inner-product (see Section 4). For other choices of inner product inherited from finite-temperature correlation functions see discussions in [9,18].
space dimension for integrable systems, focusing on SYK 2 . In Section 4 the Lanczos algorithm by which an orthonormal basis for Krylov space as well as the Lanczos-sequence are constructed, is described. In Section 5 we present the definitions of K-complexity and Kentropy as well as expectations for their saturation values. Section 6 shows our numerical results for the Lanczos-sequences, K-complexity and K-entropy for complex SYK 4 systems with 8, 9 and 10 fermions. In Section 7 we summarize the main points of our results, remark on some of the differences between K-complexity and circuit complexity and on possible holographic aspects. Appendix A presents numerical results for the Lanczos-sequence in Random Matrix Theory. Appendix B discusses the Krylov basis for SYK 2 and Appendix C shows explicitly the numerical algorithms used to produce the data in the paper.

Time evolution of operators and physical properties of Krylov space
We begin by providing the definition of Krylov space and stressing its relation to the structure of the spectrum of the system under consideration (and hence to its integrable or chaotic character in a convenient matrix representation, which we now define. Let us choose the basis |ω ab ) ≡ |E a E b |, on H, which is naturally induced by the eigenbasis of the Hamiltonian on H. Then each nested commutator takes the form: where we have defined the phases 4 ω ab := E a − E b , which are eigenvalues of the Liouvillian acting on |ω ab ). Note that the ω's related to pairs of the form (a, a) are zero, ω aa = 0 for all a = 1, ..., D. Displaying the coordinates of (2.4) as rows one can construct a matrix representation of (2.3), which turns out to be a Vandermonde matrix. The rank of this matrix is at most D 2 , as it has D 2 columns, so we shall keep the same number of rows to make it a square Vandermonde matrix: To find its actual rank, we can compute its determinant, which is given by: where ∆ ({ω ab }) is the Vandermonde determinant of the phases in the matrix. The expression (2.6) will be zero if any of the phases are degenerate, and also if any of the matrix elements in the energy basis vanish, so the corresponding columns should be removed in order to be left with a matrix of maximal rank. Hence, the Krylov dimension K can be estimated using the following algorithmic prescription: K is equal to the number of distinct phases corresponding to the indices of non-zero matrix elements of the operator in the energy basis. The zero phase ω aa = 0 is always at least D times degenerate and therefore, since it can only be counted once, the Krylov dimension is bounded by: for any non-vanishing operator. In general, if the operator O has a non-vanishing projection over several eigenstates of the Liouvillian that share the same degenerate eigenvalue, then this phase only contributes one dimension to the Krylov space. So a legitimate question is to wonder what particular linear combination of those eigenstates is actually contained in the Krylov space. In order to answer this, one can consider the form of the operator L n |O) given in (2.4). Suppose that, for some set I of pairs of indices, the eigenvalue is degenerate: i.e. ω ab = ω for all (a, b) ∈ I. Assume also that ω ab = ω for any (a, b) / ∈ I. Inserting (2.8) in (2.4) one finds: It is manifest now that the direction of the ω-eigenspace that contributes to the Krylov space is precisely the projection of the operator O over such an eigenspace: Let's call |K ω ) the eigenspace representative for the phase ω. The structure of the Krylov space is now fully understood: Each eigenspace of the Liouvillian over which the operator O has a non-vanishing projection contributes one Krylov dimension, through a linear combination of the basis of the eigenspace of the form (2.10). We can thus redefine the Krylov space as: where σ(L) denotes the spectrum of L. Finally, the Krylov dimension K is simply equal to the number of non-zero eigenspace representatives |K ω ). For instance, if one considers a system whose Liouvillian has a spectrum with no degeneracies (other than the unavoidable null phases) and an operator that is dense in the energy basis, the Krylov dimension will be maximal, K = D 2 − D + 1. The only source of degeneracy in the phases will be the universal one due to diagonal phases ω aa = 0. We can now note explicitly that the part of the operator algebra that is left out of the Krylov space belongs to the space subtended by the projectors on the energy eigenstates, since those correspond to zero phases, and the only combination of them that contributes to H O is the representative of the ω = 0 eigenspace, in this case: (2.12) To summarize, K is determined by exploiting the advantages of the basis in operator space induced by the eigenbasis of the Hamiltonian; crucially, the determinant of a suitable matrix containing the representation of these nested commutators in the given basis reduces to a Vandermonde determinant of energy differences. The following algorithmic prescription is derived: K is equal to the number of eigenspaces of the Liouvillian over which O has non-zero projection. The eigenvalues of the Liouvillian are precisely all possible energy differences, ω ab = E a − E b . The zero phase is always at least D times degenerate, thus the Krylov dimension is bounded by (2.7), for any non-zero operator. The more degeneracies there are in the spectrum of the Liouvillian, the lower will be the Krylov dimension. The expectation is that the spectrum of a chaotic system will not feature degeneracies apart from those induced by the presence of extra symmetries; we therefore conjecture that typical 5 operators in chaotic systems saturate the upper bound in (2.7); for integrable systems we expect K to be significantly smaller than the bound. We have confirmed both expectations by studying the chaotic SYK 4 model (see numerics below) as well as RMT (see Appendix A), and the integrable SYK 2 model, see Section 3 below.

Integrable systems and SYK 2
We have conjectured that the Krylov space for the time evolution of a typical operator for integrable systems will be significantly smaller than the upper bound in (2.7). Before verifying our conjecture in the case of SYK 2 , we verify it in the simpler example of the quantum harmonic oscillator H = ω(a † a+1/2) with the position operatorx = 1 2ω (a+a † ), for which the Krylov space dimension is K = 2, even though the Hilbert space is infinite 6 . In this case, the smallness of the Krylov space is due to the fact that the position operator has a non-vanishing projection over only two eigenspaces of the Liouvillian (those corresponding to the energy differences ω and −ω).
We now test our conjecture in SYK 2 , an integrable system with more structure than the harmonic oscillator. For simplicity, let us use the Majorana (real) version of the model [21,23]: where the coupling strength m ij ∈ R is antisymmetric, m ij = −m ji , and each independent matrix element is drawn from a Gaussian distribution with zero mean and variance E(m 2 ij ) = m 2 /L. The number of sites is even, L = 2M , and the Majorana fermions satisfy the relations The key point in the SYK 2 case is the fact that operators do not grow after commutation with the Hamiltonian, as also observed by [16,33]. Hence the subspace containing their time evolution will be, at most, that subtended by operators of fixed equal size. For example if the operator of study is a Majorana on some site with index A, O ≡ χ A , commutation with the Hamiltonian will give [χ i χ j , χ A ] = δ Aj χ i − δ Ai χ j , which is again a single-site operator (more details on the construction of the Krylov space in this case are given in Appendix B). For one-site operators, this implies an upper bound on K: It should be possible to reach this conclusion by studying the degeneracy structure of the spectrum of SYK 2 , which features Poisson level spacing statistics, and it is natural to expect other integrable systems to have a small Krylov space, due to the expected degeneracies in the spectrum of the Liouvillian and to the sparseness of simple operators in the energy basis. In general, interacting integrable systems feature Poisson level spacing statistics, implying, at least, the existence of quasi-degenerate levels in the spectrum of the Hamiltonian, and therefore also in that of the Liouvillian. We expect these systems to behave as if the degeneracies were exact to a good degree of approximation, hence admitting a description in terms of an effectively smaller Krylov space. This argument applies, and accounts for a significant effective reduction of Krylov space, even if the operator is not sparse in the energy basis, as long as the invariant subspace of the Liouvillian over which it has a non-trivial projection contains degenerate or quasi-degenerate levels.

Lanczos algorithm
Once the Krylov space adapted to time evolution of an operator with a constant Hamiltonian is identified, one would like to construct an orthonormal basis for it, given a certain scalar product (·|·) on operator space. This is achieved with the Lanczos algorithm, which is a particularization of the Gram-Schmidt procedure: 5. if b n = 0 stop; otherwise set |O n ) = 1 bn |A n ) and go to step 3. In this work we make use of the Frobenius product (A|B) = 1 D Tr A † B . The algorithm will construct an orthonormal set {O n } K−1 n=0 , the Krylov basis, and the Lanczos coefficients {b n } K−1 n=1 . Each Lanczos step produces an element |A n ) orthogonal to all previous |O m ) with m < n, so it is either zero or a new direction in Krylov space. For n < K, |A n ) = 0 because the set that is being orthogonalized with this procedure has rank K (in particular, |A n ) contains terms with up to n nested commutators of H with O). However, |A K ) is orthogonal to {O n } K−1 n=0 , which is already a complete orthonormal basis of H O , so it must therefore vanish, just as b K = (A K |A K ) = 0. We conclude that the Lanczos algorithm must terminate by hitting a zero once all directions in Krylov space have been exhausted. This is accounted for in Step 5 above.
The representation of the Liouvillian over the Krylov space in such a basis simplifies to a tridiagonal matrix (O m | L |O n ) = T mn , whose entries are given by the Lanczos coefficients: The eigenvalues of this matrix are all non-degenerate, and given precisely by the phases corresponding to the eigenspaces of the Liouvillian that span the Krylov space (see Section 2). The original Lanczos algorithm described above is known to suffer from numerical instabilities which can be overcome using the re-orthogonalization algorithms FO and PRO [34,35] described in Appendix C.

K-complexity and K-entropy
Using the sequence of Lanczos coefficients (which we will also call the b-sequence), one can reduce the analysis of the time-evolution of an operator O into the solution of a differential recurrence equation. The time-evolving operator can be expanded in the Krylov basis as: where ϕ n (t) are time-dependent coefficients which describe how the operator is distributed over the Krylov basis (they may be thought of as the "wavefunctions" in n). Given the Heisenberg equation dO dt = i[H, O], they satisfy the differential recurrence equatioṅ Here, ϕ n (0) = δ n0 , since for a normalized operator O(0) = O 0 . We set ϕ −1 (t) ≡ 0 in order to make the recurrence (5.2) consistent with the definition (5.1). Also, from unitarity K−1 n=0 |ϕ n (t)| 2 = 1. The Lanczos coefficients {b n } K−1 n=1 can be understood as hopping amplitudes for the initial operator O 0 to explore the "Krylov chain" and the functions ϕ n (t) can be visualized as wave-packets travelling on it [9].
We now display two quantities which highlight broad features of the distribution ϕ n (t), viz.
• K-complexity, which computes the average position of the distribution on the ordered Krylov basis: • K-entropy, which computes how randomized the distribution is: Assuming that at very late times the operator is evenly distributed over Krylov space, there exists a saturation time t sat for which We can then get a rough estimate for the values of K-complexity and K-entropy at very late times by plugging (5.5) into the formulas for C K (t) and S K (t): Since for chaotic systems K ∼ D 2 and in general D ∼ e S , where S is the entropy of the system (in the sense of "log of the number of states"), we find that the saturation value of K-Entropy is essentially of order S, while the saturation value of K-complexity is of order e 2S . If C K (t) grows linearly after scrambling, the saturation time will be roughly t sat ≡ t K ∼ e 2S , in agreement with the expectation in [10]. These properties are confirmed in our numerical results.

Numerical results
In the following we present numerical results for complex SYK 4 , whose Hamiltonian is schematically given by: where i, j, k, l = 1, 2, ..., L (L is the system size). J ij;kl is a complex matrix whose independent elements follow normal distributions with: where the overline denotes average over random realizations. It is worth recalling that the Hilbert space of real (Majorana) SYK with L sites scales as 2 L/2 , while in the complex case it scales as 2 L . However, the advantage of using this version of SYK [36] lies in the fact that the total number operator n := L i=1 c † i c i commutes with the Hamiltonian, allowing us to work in sectors of the Hilbert space with fixed occupation number, denoted by N (eigenvalue of n). For our numerical computations we will take N = L 2 . The (hermitian) observable chosen for numerical computations is the hopping operator between sites L and L − 1:  Full time range computed. Saturation occurs at time scales of order t K ∼ K ∼ e 2S with value near K/2 = 31626.5 (this is the time scale at which the wave-packet ϕ n (t), which propagates at roughly constant -but slowly decreasing -velocity, reaches the edge of the Krylov chain). Right: Zoom in at early times. Note the initial non-linear growth transforming into linear growth, as verified by the linear fit (in red). Relevant time scales are indicated in the plots, although their exact location may depend on a dimensionful prefactor (the Lyapunov exponent, as discussed in [10]).
Any other hopping operator h ij should give results equivalent to (6.3), due to the nonlocal character of the Hamiltonian (6.1). In general, one could choose other non-extensive operators, such as the number operator at a particular site n i := c † i c i . Both n i and h ij have been shown numerically to satisfy ETH [36].
We studied samples with L = 8, 9 and 10 sites. The numerical results for SYK 4 with L = 10 are shown in Figures 1, 2 and 3. For details on the scaling properties of the bsequence, K-complexity and K-entropy with system size, see Figures 4, 5, 6 and 7, as well as the summarizing Table 1.
The global picture emerging from these results can be summarized in the following bullet points: • The Descent and its associated Lanczos sequence. The length of the Lanczos sequence saturates its upper bound. It features a period of initial growth up to n ∼ S, followed by a regime of slow decrease to zero with roughly constant negative non-perturbative slope of order ∼ −e −2S , the Descent.
• K-complexity features a transition from exponential growth at very early times to linear increase. At exponentially late times it saturates at half of the Krylov dimension, since by then the operator is uniformly distributed over the Krylov basis.
• K-entropy grows linearly up to scrambling time, and then transitions to a logarithmic growth phase that continues until saturation around S K ∼ S at exponentially late times.
The relationship between the Lanczos sequence b n and quantities like C K (t) and S K (t) is highly non-linear, which is why disorder averages need to be performed with caution. The Lanczos sequence is not really a physical observable and averaging is just used as a tool to gain knowledge about the envelope of its profile. However, in order to obtain averaged K-complexity or K-entropy as a function of time, one should compute these quantities separately for each random realization of the Hamiltonian and average over the outcomes in a final step. Averaging over the Lanczos sequence before computing C K and S K results in a smoothing of the b-sequence that stops the wave packet from randomizing efficiently before reaching the edge of the Krylov chain, being therefore reflected due to the boundary conditions. These rebounds would be reflected by large, un-physical oscillations in the profiles of C K and S K , an effect that has been confirmed by numerical simulations.

Discussion
On the analytical side, we have found an algorithmic expression for the dimension of Krylov space K. This quantity is very sensitive to the degeneracy structure of the spectrum and is always strictly smaller than the dimension of operator space D 2 , since it satisfies K ≤ D 2 − D + 1. Typical operators in chaotic systems are expected to saturate this bound, due to the absence of degeneracies, whereas in integrable systems K will typically be significantly smaller than its upper bound (at least effectively up to late times), as verified analytically for SYK 2 . Furthermore, numerical observations in SYK 4 and RMT suggest that the b n -profile featuring a slow decrease to zero given by a non-perturbative slope of order e −2S -the Descent -may be a generic feature of quantum chaotic systems. It is worth noting that, since K-complexity is nothing but an average position in Krylov space, the above bound on K translates directly into an upper bound for C K (t). We note that although K-complexity features the profile expected from circuit complexity, there are notable differences between them: Krylov complexity does not depend on an arbitrary tolerance parameter, and it does depend on both the Hamiltonian of the system and the operator whose K-complexity is measured over time. This may have impact in situations that involve more than one operator insertion. We hope to address this and other related  issues in the future.
We end by returning to an issue we emphasized at the outset, namely that K-complexity is naturally bounded. This is a direct consequence of the finite dimensionality of the Hilbert space, and thus has the same origin as the plateau in the spectral form factor and correlation functions. We conjecture that other late-time universal characteristics of quantum chaos, namely the dip-ramp-plateau structure [7], also leave their imprint on K-complexity. The putative bulk realization of K-complexity should therefore be sensitive to bulk Euclidean wormholes and baby universes [37], in ensemble-averaged systems, such as SYK, as well as in individual quantum chaotic systems, as has been emphasized in [17].
Note added: in the final stages of preparation of this paper, the work [38] appeared on the arxiv studying, among other things, some aspects of K-complexity in SYK. Up to the timescales they consider, our results are consistent with theirs.

Acknowledgments
We would like to thank José Barbón as well as D. Aharonov, M. Ben-Or, N. Katz, P. Nayak, H. Neuberger, J. Rougemont and M. Vielma for enlightening discussion. ER would like to thank NHETC at Rutgers Physics Department and CCPP at NYU for hospitality. The numerical computations in this paper were performed on the Baobab HPC cluster at the University of Geneva and on the Landau cluster at the Hebrew University. This work has been supported by the SNF through Project Grants 200020 182513, the NCCR 51NF40-141869 The Mathematics of Physics (SwissMAP). The work of ER and RS is partially supported by the Israeli Science Foundation Center of Excellence. We thank the organisers of the conference "Frontiers in Holography" (May 2020, Moscow), where the results of this work were presented.

A Lanczos sequence in RMT
Some preliminary numerical checks have indicated that RMT reproduces qualitatively the features observed in complex SYK with q = 4, that is: saturation of the upper bound for Krylov space, K = D 2 − D + 1, and slow decrease of the b-sequence to zero after initial growth, with a non-perturbative slope of order ∼ − 1 K ∼ −e −2S , S being the entropy (system size). Figure 8 depicts the Lanczos sequence of a system whose Hamiltonian H is drawn from a GUE with potential: where D is the Hilbert space dimension and J is the coupling strength that sets energy scales. A more detailed numerical and analytical study of the Lanczos sequence in RMT constitutes work in progress, but these preliminary checks make it natural to conjecture that the discussed features are universal in chaotic systems.

B Krylov space for SYK 2
Given the SYK 2 model described in 3.1, we can proceed to apply the Lanczos algorithm analytically to the operator O = χ A . We will use the Frobenius scalar product defined in the main text, recalling that the dimension of the Hilbert space of states in this case is D = 2 M . It is not difficult to prove that: ( • O 0 : Making use of (B.1) one finds that our starting operator is not normalized, since (O|O) = 1 2 , so: • O 1 : One first needs to compute A 1 : We make use of the commutator of the fermionic bilinear with the single Majorana: so that finally The first Lanczos coefficient is: and recalling (B.1) one finds: Thus, the next Krylov element is: The form of (B.7) allows the computation of the ensemble average of the first Lanczos coefficient, since the calculation will turn out to be simple if one makes use of spherical coordinates in sample space: where we defined b(x) ≡ √ x 2 , and note that there are L−1 random variables because (B.7) only involves L−1 independent coupling strengths. The probability distribution is a product of Gaussians: where, in agreement with E(m 2 ij ) = m 2 /L, the standard deviation is given by σ = m √ L . As promised, the use of spherical coordinates greatly simplifies the computation and the result is: (B.11) • More elements: It can be an interesting task to compute in closed form the rest of the Lanczos sequence, but we can already anticipate what its length will be: As shown in (B.4), the commutator of a bilinear with a single-site operator will still yield a linear combination of Majoranas on a single site, so still be a one-site operator, and so will the rest of the Krylov elements. Hence the maximum possible number of linearly independent Krylov elements will be given by the number of sites of the system, that is to say, in this integrable system the Krylov dimension K is bounded by: i.e. the Krylov space scales at most linearly with entropy (system size), instead of with the operator space dimension (exponential in entropy).

C Numerical algorithms
It is known that the original Lanczos algorithm, as presented in the main text, features an important numerical instability [34,35,39]: the construction of each Krylov element makes use of the two previous ones, so errors due to finite-precision arithmetic accumulate dramatically and orthogonality of the Krylov basis is soon lost in numerical computations. Residual overlaps between Krylov elements grow exponentially (or even faster) with the iteration number n, which makes the Lanczos coefficients unreliable after a few iterations.
In particular, the original Lanczos algorithm doesn't feature the termination of the sequence by hitting a zero at n = K when run with finite precision, and instead it generally outputs a Lanczos sequence that oscillates wildly around some constant value, whose disorder average yields a completely flat sequence (after initial growth) in complex SYK; this is purely a product of numerical precision errors and needs to be corrected in order to shed light on the structure of the b-sequence along the full length of the Krylov chain. Some numerical algorithms need to be implemented in order to cure this instability, allowing observation of the slow decay to zero of the Lanczos coefficients after the initial growth in complex SYK. Such algorithms are described below.

C.1 Full Orthogonalization (FO)
This algorithm [35] performs a brute-force re-orthogonalization of the newly constructed Krylov element with respect to the previous ones at every iteration of the Lanczos algorithm, ensuring orthonormality of the Krylov basis up to machine precision ε M . The algorithm is not very efficient time-wise, and is also costly in terms of memory, since the whole Krylov basis needs to be stored and is used in every iteration, it is therefore used mainly when one wants to compute only part of the Lanczos-sequence (see for example the recent work [40]). However, for small samples, FO can be used safely and one can check that the results yielded agree with the theoretical predictions described previously in this article (length of the Lanczos sequence, termination by hitting a zero, eigenvalues of the tri-diagonal matrix matching those corresponding to the eigenspace representatives that span the Krylov space). The FO algorithm amounts to performing explicit Gram-Schmidt at every iteration in the Lanczos algorithm to ensure orthogonality (up to machine precision). For numerical purposes, it is usually optimal to perform Gram-Schmidt twice every time: 2. For n ≥ 1: Compute |A n ) = L |O n−1 ).
5. Set b n = (A n |A n ).
6. If b n = 0 stop; otherwise set |O n ) = 1 bn |A n ) and go to step 2.

C.2 Partial Re-Orthogonalization (PRO)
This algorithm [34] allows the residual overlaps between Krylov elements to grow up to a certain threshold, and re-orthogonalization is only performed when the threshold is crossed. For a machine precision ε M , the threshold is typically taken to be √ ε M .
In our notation, the recursion relation for the Krylov basis is, including finite-precision errors: b n |O n ) = L |O n−1 ) − b n−1 |O n−2 ) + |ξ n−1 ) (C.1) where |ξ n−1 ) accounts for some spurious vector generated by accumulated numerical errors. All the objects denoted above represent the quantities numerically computed, rather than the actual (analytically exact) Lanczos coefficients and Krylov elements. Acting with (O k | from the left: where we recall that (O k | L |O n−1 ) = T k,n−1 , T being the tri-diagonal symmetric matrix built out of the Lanczos sequence, see (4.1). We now define the matrix W , with items W kn = (O k |O n ), that is: Even though we don't write all the entries in the matrix, it's hermitian by definition of the scalar product (·|·), so W = W † ⇐⇒ W kn = W * nk . In the PRO algorithm, however, we construct iteratively the entries written explicitly in (C.3): For a given n we'll want to estimate W kn , for k ≤ n. We do so by noticing that (C.2) is nothing but: b n W kn = T k,n−1 − b n−1 W k,n−2 + (O k |ξ n−1 ) . (C.4) Renaming indices k ↔ n − 1 (i.e. k → n − 1 and n → k + 1): b k+1 W n−1,k+1 = T n−1,k − b k W n−1,k−1 + (O n−1 |ξ k ) . (C.5) Computing (C.4) − (C.5), recalling that T is symmetric and solving for W kn : (C.6) We want to use (C.6) to determine, for a fixed n, all W kn with k ≤ n given that we know all {W ij , i = 0, ..., j, ∀j = 0, ..., n − 1} (according to what we depicted in (C.3): we compute iteratively each upper-diagonal column making use of the previous ones).
We note that W nn is not determined by (C.6) in terms of previous upper-diagonal columns, but we can set it to W nn = 1 because in the n-th Lanczos step, |O n ) is explicitly normalized to unity. Likewise, W n−1,n is not determined from (C.6) in terms of previous columns, and the Lanczos recursion (C.1) doesn't explicitly orthogonalize |O n ) against |O n−1 ), so we will need to orthogonalize them explicitly, and then set W n−1,n = O (ε M ) (i.e. zero up to machine precision).
Also, an estimate for (O k |ξ n−1 ) − (O n−1 |ξ k ) is needed. One can take something of order of the machine precision times the norm of the Liouvillian: where L should be the norm of the Liouvillian induced by the scalar product (·|·) in the Hilbert space of operators. For practical applications, this contribution can just be ignored, since in any case each iteration of the algorithm will already generate spurious errors of order ε M . All in all, the LanPRO algorithm reads: |O).
-Compute the a-priori Lanczos coefficient: b n = (A n |A n ).
Some comments are in order: • In every iteration n ≥ 2, only the two previous upper-diagonal columns {W k,n−1 , k = 0, ..., n − 1} and {W k,n−2 , k = 0, ..., n − 2} of the matrix W are required. This is why, whenever re-orthogonalization is required, one only carries it out for |A n−1 ) and |A n ).
• Whenever re-orthogonalization is needed, it is optimal to perform it twice, in the same way Gram-Schmidt is applied twice at every step of the FO algorithm.
In the case of complex SYK 4 , for the system sizes studied in this article and a typical floating point precision of ε M ∼ 10 −15 , PRO reduces the number of re-orthogonalizations required by approximately a factor of 10, as compared to FO.