Krylov complexity and orthogonal polynomials

Krylov complexity measures operator growth with respect to a basis, which is adapted to the Heisenberg time evolution. The construction of that basis relies on the Lanczos algorithm, also known as the recursion method. The mathematics of Krylov complexity can be described in terms of orthogonal polynomials. We provide a pedagogical introduction to the subject and work out analytically a number of examples involving the classical orthogonal polynomials, polynomials of the Hahn class, and the Tricomi-Carlitz polynomials.


Introduction
A key concept of (quantum) information theory is complexity, which incarnates, in all of its variants, the intuitive understanding that an object is the more complex the more ingredients are required to assemble, create, or describe it. It has been recognized in recent years that complexity embraces areas as distant as computer science, chaos, emergent phenomena in many-body systems, and black holes. Many of these efforts within quantum field theory and quantum gravity, as well as an outline of some remaining challenges, have been summarized in the snowmass white paper [1].
In a classic paper, Kolmogorov [2] proposed to define as the "relative complexity of an object y with a given x, [. . . ] the minimal length l(p) of the program p for obtaining y from x." Kolmogorov also pointed out that such an algorithmic approach transcends the notion of entropy for quantifying information content. The minimum number of operations necessary to implement a given task is now known as computational complexity in computer science and information theory [3]. Similarly, quantum complexity is defined as the minimum number of elementary operations (quantum gates) needed to build a state |ψ from a given reference state |ψ 0 . A quantum circuit in operation increases the complexity of an initially simple state scrambling its information over a large number of degrees of freedom. The rate of complexity growth is expected to satisfy a physical bound proportional to the portion of energy that is available for information processing [4,5]. Black holes are considered the fastest scramblers in nature [6][7][8] and are deeply connected to quantum chaos [9,10]. Computational complexity has also been added to the holographic dictionary [5,8,11].
The same physical mechanism that scrambles information is responsible for the emergence of irreversible macroscopic behaviour such as thermalization or hydrodynamics [12][13][14][15]. Precisely quantifying the growth of operators [16,17] is important to understand such emergent phenomena.
A frequently used definition of operator size builds upon out-of-time-ordered correlators (OTOCs), which measure the amount a (time-dependent) operator fails to commute with other, simple, operators and provide a diagnostic of quantum chaos [18][19][20]. OTOC-based operator size has been discussed in [21,22] for the Sachdev-Ye-Kitaev model [23][24][25], which has a classical gravity dual and is a benchmark model for a fast scrambler. 1 Krylov complexity, or K-complexity for short, was introduced in [27] as a measure of operator size with respect to a basis that is well adapted to Hamiltonian time evolution. In this context, the authors of [27] proposed a Universal Operator Growth Hypothesis for systems with chaos, according to which the spectral density decays exponentially as e −|ω|/ω 0 for large frequencies |ω|, with a decay constant ω 0 related to the Lyapunov exponent. Within the SYK model and when working in the thermodynamic limit, Krylov complexity is similar to other notions of operator size. However, whereas operator size fails to be a good measure of complexity beyond the thermodynamic limit, i.e., at time scales much larger than the scrambling time, Krylov complexity continues to show the 1 For a pedagogical introduction to SYK, see [26].
typical features of complexity [28,29]. In particular, it grows linearly in the post-scrambling period until it saturates because of the ultimately finite dimensionality of the operator space, of order e O(S) .
The main advantage of Krylov complexity over OTOCs and (quantum) computational complexity is that it is free of ambiguities and provides a truly intrinsic measure of operator and complexity growth. For example, quantum complexity crucially depends on the choice of elementary gates (or operations) and a tolerance parameter , which describes the accuracy with which the target state |ψ is reached. Similar ambiguities exist in Nielsen's geometric approach to quantifying complexity [30][31][32], which requires a choice of an information metric penalizing hard vs. easy operations.
OTOCs depend on a reference operator. In contrast, the only freedom in Krylov complexity is the choice of the inner product in operator space, but this freedom is usually limited by the physical context. Once the inner product has been defined, the construction of Krylov space proceeds by means of the Lanczos algorithm, also known as the recursion method [33], and Krylov complexity may be understood as the effective dimension of Krylov space.
The seminal work [27] has spurred quite a number of applications, which we would like to mention here. Krylov complexity has been studied for the SYK model, besides in [27,29], also in [34]. The relation to random matrix theory and black holes has been investigated in [35]. The authors of [35] suggested a microcanonical inner product, which can be used to define a thermal Krylov complexity. We point out that the complexity renormalization group framework introduced in [35] resembles the introduction of a square root terminator function in the recursion method [33].
It has also been found that quantum chaos can be understood by delocalization in Krylov space [36].
In contrast, Krylov complexity is suppressed in integrating integrable models of finite size [37], which has been called Krylov localization. A number of simple systems, which are symmetry generated and allow an exact treatment using generalized coherent states, were considered in [38,39].
Further tests, numerical calculations, applications and generalizations include chaotic Ising chains and systems with many-body localization [40,41], a variety of exemplary systems, including 1d and 2d Ising models as well as 1d Heisenberg models [42], operator growth in conformal field theory [43,44], lattice systems with local interactions under Euclidean time evolution [45], the emergence of bulk Poincaré symmetry in systems exhibiting large-N factorization [46], topological phases of matter [47], integrable models with saddle-point dominated scrambling [48], cosmological Krylov complexity [49], a condition for the irreversibility of operator growth [50], complexity in We add the exact calculation of complexity, whenever possible. Examples involving the orthogonal polynomials of the Hahn class, i.e., the Charlier, Krawtchouk and Meixner polynomials, are worked out in section 5. The complexities for these cases cannot be understood from a scaling regime and, somewhat surprisingly, are very similar to each other. Moreover, they appear to resemble the numerical results of systems with many-body localization [41]. The non-standard Tricomi-Carlitz polynomials are discussed in section 6. These are interesting, because in spite of the fact that the Lanczos coefficients obey a scaling law the continuum approximation is suspected to fail. Unfortunately, we are not able to obtain a useful analytical result, but we leave a numerical investigation for the future. Finally, in section 7, we review and discuss several operator inner products.
2 Operator growth in Krylov space

Recursion method
In the Heisenberg picture, an initial observable O evolves with time into is the Liouvillian superoperator. Clearly, O(t) describes a trajectory in the vector space of all hermitean operators and becomes, generically speaking, increasingly more complex. In order to quantify this growth of complexity it is necessary in the first place to describe in more appropriate terms the space O(t) moves in. An intrinsic approach for doing this is to construct an operator basis using the Lanczos algorithm, also called the recursion method in physics [33]. 2 The recursion method can be summarized as follows. A necessary preliminary step is to promote the space of operators to a Hilbert spaceĤ by endowing it with an inner product (A|B) satisfying In addition, the Liouvillian must be hermitean with respect to the inner product, (A|LB) = (LA|B) . (2.4) For simplicity, we assume O to be normalized, The choice of the inner product is quite important and will influence the outcome of the recursion method, but we do not need to be more specific for the purpose of developing the general method.
We will postpone the review of a variety of inner products to section 7.
Expanding the Heisenberg evolution (2.1) into a formal power series in time shows that the time dependent operator |O(t)) is constructed out of the sequence of states {|L n O), n = 0, 1, 2, . . .}.
The Lanczos algorithm essentially performs a Gram-Schmidt orthogonalization on this sequence, producing an (ordered) orthogonal set of states {|O n ), n = 0, 1, 2, . . .} called the Krylov basis. The space spanned by the Krylov basis is called the Krylov (sub)space, which is a subspace ofĤ.
In table 1, we spell out two variants of the Lanczos algorithm, which we call the orthonormal and the monic versions, with reference to the kinds of polynomials they produce. Both versions can be found in the literature and produce, obviously, equivalent results. Let us add a few comments.
First, the positive numbers b n are called the Lanczos coefficients. The relation between ∆ n and b n is simply Second, if the iteration stops, then the Krylov subspace is finite-dimensional. Third, it is known that the Lanczos algorithm is numerically unstable. Since our focus is on analytical results, we will ignore this issue and refer interested readers to [29] and references therein.
In either version of the recursion method, one has by construction with (O n |O n ) = 1 in the orthonormal version. More importantly, the procedure gives

Lanczos algorithm (recursion method)
orthonormal version monic version where P n (L) denotes a certain polynomial of degree n in L, with real coefficients. It is evident from the table that these polynomials satisfy a three-term recurrence relation, 3 Because the eigenvalues of L are real (recall that L is required to be hermitean with respect to the inner product), by Favard's theorem [57,58] there exists a measure µ(L) on R such that the P n are orthogonal with respect to µ(L), dµ(L) P m (L)P n (L) = h n δ mn (h n > 0).
The support of µ(L) is identified with the spectrum of L, and integrals such as the one appearing in (2.11) are intended as integrals over the eigenvalues of L. Consider the (normalized) measure The fact that there is no term with Pn on the right hand sides of (2.9) and (2.10) derives from the fact that the initial operator O is assumed to be hermitean. It follows by induction that also (iL) n O is hermitean, from which one can show that (L n O|L m O) = 0 for m + n odd. 4 Favard's theorem does not guarantee that the measure µ(L) in (

Wave functions
Given the Krylov basis, define a semi-infinite chain of wave functions by 5 which are real, because both i n O n and O(t) are hermitean operators. 6 These wave functions have the equivalent forms By virtue of (2.9), they satisfy the Schrödinger equation In other words, the wave functions φ n constitute a semi-infinite tight-binding model, for which the Lanczos coefficients represent the hopping amplitudes [36]. Moreover, they satisfy n φ n (t) 2 = 1 . (2.17) The easiest way to see this is to note that (2.17) holds by definition for t = 0 and that the derivative with respect to time vanishes by virtue of (2.16).
In contrast to the construction of the Krylov basis, (2.16) cannot be solved by iteration. A formal solution can be given as follows [33]. Consider the Laplace transform of φ n (t), In terms of these, (2.16) translates to 5 hn has been included to make the φn independent of the chosen normalization. 6 Indeed, whereas [33] uses (2.14), the complex conjugate definition is often found in the literature (e.g., in [27,35,38,44]) and gives, obviously, the same result.

After introducing
where (2.6) has been used. R 1 (z) is also called the memory function of the system [33]. In turn, (2.21) implies the continuous fraction representation of c 0 (z) (2.22) The functions c n (z) are closely related to the functions of the second kind [59,60] Q n (z) = dµ(L) More precisely, which also provides the analytic continuation of (2.18) to z < 0. In particular, is called the resolvent. It contains the spectral function The functions of the second kind satisfy, except for n = 0, the same three-term recurrence relation as the polynomials P n , which is easily proven from (2.23) and (2.10). In the case of monic P n , one We conclude this subsection by remarking that the entire information about the operator O(t) is contained equivalently in any of the following: • the set of Lanczos coefficients (b n or ∆ n ), • the fluctuation function φ 0 (t) or its Laplace transform c 0 (z), which is equivalent to the function of the second kind Q 0 (z), • the spectral density, which is related to Q 0 (z) by (2.26), • the moments of the spectral density, which coincide with the coefficients of the expansion of φ 0 (t) in a series in t 2 . There are recursive algorithms translating the set of moments into the set of Lanczos coefficients and vice versa. These are useful for numerical calculations with known spectral functions, but we shall not need them here. Interested readers are referred to [33].

Consequences of the even spectrum
In this subsection, we will list some simple, but important, properties, which derive from the fact that the recurrence relations (2.9) and (2.10) do not contain a term with P n on their right hand sides. This implies that the measure µ(L) is even, µ(L) = µ(−L). To be specific, let us work with the monic polynomials.
A general property of orthogonal polynomials with respect to an even measure is that P 2n (L) and P 2n+1 (L) can be expressed in terms of orthogonal polynomials in L 2 with appropriate measures.
In this case, one formally gets a contribution dμ(L 2 ) = 2aδ(L 2 ) dL 2 + · · · . In order to avoid double counting, one can adopt the convention we also have We remark that the measuresμ andμ are normalized, i.e.,h 0 =h 0 = 1.
Both,P n andP n , satisfy three-term recurrence relations, which can be obtained from (2.10).
Multiplying (2.10) by L and using it once more one finds Moreover, they are interrelated byP The functions of the second kind (2.23) are where the definition ofQ n andQ n is analogous to (2.23). They are interrelated by Finally, consider the wave functions (2.15). With (2.28) and (2.29), they become where (2.32) has been used.

Measures of operator growth
Krylov-or K-complexity has been defined in [27] as the time-dependent mean position on the semi-infinite chain of wave functions, This definition can be generalized to an entire class of measures of complexity growth. We define Krylov complexity of degree k as and the Krylov entropy [28,29], The recurrence relation of the wave functions (2.16) allows to obtain which, taking into account the boundary conditions It is useful to define the generating function from which the Krylov complexities of various degrees can be obtained by λ is formal parameter, which one can take negative in order to ensure that the sum in (2.46) is convergent.
, it is often more useful to consider its Laplace transform, With the help of (2.15), this can be expressed as which, after using (2.23), becomes The term independent of λ is just the Laplace transform of unity, i.e., the norm (2.17). The term linear in λ gives the Laplace transform of the complexity, (2.53) Applying to (2.53) the same trick as above, one gets This is, of course, nothing but the Laplace transform of (2.45).

Generic features of operator growth
In this section, we shall discuss a number of generic features, which follow directly from general properties of orthogonal polynomials or can be derived easily from the results of the previous section.

Measures with finite, bounded and unbounded support
We must distinguish three types of systems according to whether the support of the measure µ(L) is finite, bounded, or unbounded. In the case of infinite spectra, which are associated with a semi-infinite chain, the late-time complexity depends crucially on the asymptotic behaviour of the recursion coefficients ∆ n for large n. Several mathematical facts are known, some of which we will state here.
Let I = [ξ 1 , η 1 ] be the true domain of orthogonality of the polynomialsP n , which is defined as the smallest closed interval containing all zeroes y n,i of allP n [57,58]. It is also the smallest closed interval containing the support ofμ(L 2 ), i.e., supp(μ(L 2 )) ⊆ I ⊆ [0, 4E 2 ]. We recall from (2.33) that the polynomialsP n satisfy a three-term recurrence relation with coefficients Then, the following mathematical statements hold [57,58].
• Because I is bounded, the sets of both coefficients,b n andc n , are bounded. Therefore, also the set of ∆ n is bounded.
• If, for n → ∞,b n →b andc n →c, then supp(μ) is bounded with at most countably many points outside of the interval J = [b − 2 √c ,b + 2 √c ], andb ± 2 √c are limit points of supp(μ).
Consider the simplest case in which ∆ n → ∆ in the large n limit. It follows from (3.1) that b n → 2∆ andc n → ∆ 2 , so that J = [0, 4∆] with ∆ ≤ E 2 . Moreover, 0 and 4∆ are limit points of supp(μ). Hence, because 0 is a limit point, there cannot be a "gap", if by "gap" we intend the absence of states with arbitrarily small positive energy.
A single "gap" occurs, if the Lanczos coefficients for even and odd n approach different limiting values, i.e., ∆ 2n → ∆ e and ∆ 2n+1 → ∆ o . In this case, More generally, if the support of L 2 contains several disconnected intervals ("bands") within [0, 4E 2 ], then with increasing n the coefficients ∆ n will continue to oscillate in a predictable fashion [55].
Unbounded spectrum An unbounded set of coefficientsc n is sufficient to have a spectrum with unbounded support [57,58]. Then, because of (3.1), also ∆ n andb n grow without bounds. We may characterize the average growth using a power law [33], where 0 < λ ≤ 2. λ = 2 is associated with chaos, whereas systems with 0 < λ < 2 are integrable [27].
In physical applications, ∆ n may oscillate considerably around the mean behaviour. The lower limit λ = 0 corresponds to a bounded spectrum, while the upper limit λ = 2 arises from physical considerations on the spectral function [27,33], which decays as for large frequencies. With (3.3) one can associate an exponential growth of φ 0 for imaginary times [33], Mathematically, everything that can go wrong in the general formalism goes wrong for λ ≥ 2.
The fluctuation function φ 0 (t) is not an entire function (of complex time), the continuous fraction representation of c 0 (z) (2.22) does not converge, and the measure cannot be uniquely determined from the sequence of ∆ n . Therefore, with the exception of the limiting case λ = 2, we shall exclude these cases from our considerations.
We shall explicitly see in the next subsections and in the examples of section 4 that the special cases λ = 0, λ = 1 and λ = 2 are associated with linear, quadratic and exponential late-time growth of complexity, respectively.

Simple solutions for the complexity
There are only a few cases, in which one can provide an explicit and exact formula for the complexity In this subsection, we shall present three exceptionally simple cases, which coincide with the three cases related to generalized coherent states discussed in [38]. Our new derivations of (3.5) Therefore, the location of a pole of K O (z) encodes exponential or oscillating behaviour, and for the latter the poles must come in complex conjugate pairs. Moreover, the degree of that pole is related to a power law correction, or a pure power law, if the pole is at z = 0.
The simplest example is the purely linear sequence ∆ n = n∆. Let us use (2.54) to calculate The sum in (3.6) is the same as in (2.50) with λ = 0, so that the result can be read off from (2.52).
Thus, one finds This is the Heisenberg algebra case of [38]. We shall meet it again in subsection (4.4) in connection with the Hermite polynomials.
The next case is the limiting case of a quadratic plus linear growth, ∆ n = α 2 n 2 + βn. Now, The sum involving the term (α 2 + β) can be done as in the example above, while the sum with α 2 n becomes again K O (z) by virtue of (2.51). Hence, one obtains which is solved by Performing the inverse Laplace transform yields Clearly, with β = α 2 (2h − 1), the result of the SU (1, 1) case of [38] is reproduced. We remark that the orthogonal polynomials for this case are Meixner-Pollaczek polynomials with an even measure [61,62], where L = 2αx. These are orthogonal in x ∈ (−∞, ∞) with respect to the (normalized) measure For the case h = 1 2 see [63,64]. Finally, we may consider ∆ n = −α 2 n 2 + βn, which is consistent only if the spectrum is discrete, i.e., ∆ N = 0 for some N ≥ 1, which implies β = α 2 N . The calculation of K O (z) goes as in the previous case after replacing α → iα. Hence, Here, we recognize the SU (2) case of [38] with 2j = N − 1. The polynomials associated with this simple case are (shifted) Krawtchouk polynomials [61],

Continuum limit
If the Lanczos coefficients b n (or ∆ n ) can be approximated by sufficiently simple functions of n, then it is possible to estimate the late-time behaviour of the complexity using a continuum limit.
In this subsection, we shall illustrate this approach for some simple cases following [28].
To start, let us introduce with 1, where we assume that both, b(x) and φ(x, t) are sufficiently smooth functions of x.
It is clear that this assumption does not hold for the wave function φ(x, t) at t = 0, but it is sufficient for our purposes that it holds for some t > t 0 . For the rest of this subsection, we will shift t − t 0 → t. With (3.13), the discrete Schrödinger equation (2.16) can be approximated by the first-order partial differential equation where we omitted the O( 2 ) terms. Equation (3.14) can be solved by means of the method of characteristics. Using s = t to parameterize the characteristic curves, (3.14) translates into the Formally, the general solution can be expressed as where τ (x, t) is an integration constant arising in the first equation of (3.15), e.g., To take some simple specific examples, consider the scaling behaviour In the special case λ = 2, which represents a fast scrambler, the two equations in (3.15) are solved by x = x 0 e 2αt and φ = ψ e −αt , respectively, where ψ is now an arbitrary function of the integration constant x 0 . The constant α is related to the Lyapunov exponent. Thus, we can write down the i.e., Therefore, the complexity is, approximately, This agrees with the exponential growth at late times in (3.8).
For λ < 2, the solution of the first equation in (3.15) is 8 where τ is a non-negative constant. Thus, the general solution (3.16) reads Let us assume that, for t = 0, the support of φ is restricted to , which implies ψ(τ ) = 0 for τ ≥ T > 0. Therefore, the norm constraint (3.20) becomes, after a change of variable, Because the norm must be time-independent, one has the additional constraint ψ(τ ) = 0 for τ < 0.
Taking this into account, the complexity (3.22) is easily found to be For large enough t, this behaves as where (3.26) has been used.
Let us briefly comment on these results. One can observe from (3.24) that the length of the interval on which φ(x, t) has support grows with time, if 0 < λ < 2, whereas it shrinks for λ < 0.
More precisely, if φ(x, 0) has support in x ∈ [0, X) then, at late times, the size of the domain of support at time t is to leading order in T /t. This does not shrink for λ ≥ 0, which implies that we can trust the continuum approximation. As a nice check, for λ = 1, (3.28) agrees precisely with the exact result (3.7). In the limiting case λ = 0, φ(x, t) describes a wave moving with constant velocity 2 α, which implies K O (t) ≈ 2αt. As we shall see in the subsections 4.1, 4.2, and 4.3, this reproduces only the qualitative linear growth of the complexity with time, whereas the velocity is not generic. The slope 2α appears to be an upper limit. The situation is different for λ < 0. Although the solution of the differential equation remains formally correct, the continuum approximation must fail at a certain point, and we should not trust the result (3.28).
The arguments used in these simple examples can be extended to the general solution (3.16)- (3.29) and the complexity The last approximation holds for t T , and τ 0 ∈ (0, T ). Remember, however, that the time variable has been shifted in the continuum approximation with respect to the exact initial value problem. This suggests that τ 0 simply describes this shift. Extending (3.30) to model the complexity also at early times, we set τ 0 = 0, such that x 0 (0) = 0 from (3.17).
Let us illustrate this idea with two examples, which will be treated in more detail in later which gives rise to the Gegenbauer polynomials, see subsection 4.3. In this case, (3.17) yields an .
Clearly, we must restrict to β ≥ 1, as otherwise the integrand would be imaginary for small x.
For large n, the b n satisfy the scaling relation (3.18) with α = 1 2 . In agreement with the previous discussion, (3.32) shows that n(t) ≈ t for t large, but one would expect that this overestimates the velocity of complexity growth. Indeed, this will be confirmed in subsection 4.3. The integral in (3.32) can be solved in terms of elliptic integrals, but numerical integration is sufficient for our purposes. Some cases of β are shown in figure 1.
Another example is the sequence of Lanczos coefficients which gives rise to the Tricomi-Carlitz polynomails, see section 6. In this case, the characteristic curves n(t) have the implicit form where α must be restricted to α ≥ 1. As before, (3.34) can be expressed in terms of elliptic integrals, but numerical integration suffices for our purposes. Some cases are shown in figure 2. We remark that the characteristics in this case behave like n(t) ∼ t 2 3 for large t, but we should not trust the continuum approximation, as discussed above.
As one last example, let us consider the case in which the Lanczos coefficients b n approach two different constants, depending on whether n is even or odd. This case is clearly not included in the simple scaling behaviour discussed above, but it can be analyzed analytically, as we shall see in subsection 4.5. Consider the Schrödinger equation From (2.16), one easily finds the second-order differential equation which holds for both, even and odd, n (n > 1). The continuum approximation of (3.36) is the massive Klein-Gordon equation Without further calculation, we can conclude that the "pulse" φ n (t) propagates like a massive particle, with constant velocity smaller than the "speed of light". Hence, we obtain an approximate upper bound for the late-time growth,

Examples with classical orthogonal polynomials
In this section, we will illustrate the recursion method and the calculation of Krylov complexity using systems described by the classical orthogonal polynomials. Following a pedagogical spirit, we present the examples in order of increasing complexity. All of the examples we present here have appeared elsewhere, see [33,35,55] and references therein, to which we only add the calculation of the Krylov complexity. Analytic examples are also important for the construction of terminator functions, which provide suitable approximations of the continued fraction representation of the function c 0 (z) [33]. Such terminator functions have a physical interpretation as describing the complex inaccessible environment, into which the operator dynamics dissolves in the late-time limit [35]. The position N c , at which the terminator function is put sets a renormalization scale.
One may then study the renormalization group flow followed by the complexity as a function of N c .

Chebyshev polynomials of the second kind
The simplest case is given by the Chebyshev polynomials of the second kind, which are associated with constant ∆ n = ∆. They are a special case of the Gegenbauer polynomials, which shall be discussed in subsection 4.3. Here, we will present two derivations of the Krylov complexity.
Our starting point is a set of constant Lanczos coefficients. For simplicity, we use which one can achieve by re-scaling L to a dimensionless variable x. Here and in the following examples, we will use x in place of L in order to remind ourselves of this fact.
The recurrence relation (2.10) now reads which is solved by the (monic) Chebyshev Polynomials of the second kind [61] P n (x) = 2 −n U n (x) . They are orthogonal with respect to the measure and have the norms h n = 4 −n .
The wave functions (2.14) are The first way to calculate the complexity is by solving the differential equation (2.45). The sum on the right hand side of that equation simplifies to the single summand with n = 0, so that the problem to be solved is The solution is This has the late-time asymptotic behaviour Another way to calculate the complexity is via its Laplace transform (2.54). For this, we need the functions of the second kind, Q n (z), or equivalently c n (z). These are easily obtained from the recursion (2.21) by noting that one can set R n (z) = R(z) by virtue of (4.1). The recursion leads to the quadratic equation which is solved by 9 Then, (2.21) and (2.20) give It is easy to check the measure (4.4) from c 0 (z). The functions of the second kind are obtained from (2.24), Consider now the Laplace transform of the complexity (2.54). As before, only the n = 0 term contributes to the sum, so that After substituting the measure (4.4) and using the fact that it is even, this becomes (4.14) 9 We only write the physical solution.
We note that (4.14) should be the Laplace transform of (4.7), but this is not obvious to us. For small z, (4.14) reduces to 10 This agrees with the late-time behaviour (4.8).
We end this subsection by noting that the function R(z) (4.10) is nothing but the square root terminator that is often used for sequences of ∆ n that approach a constant ∆ for large n. This terminator function was advocated in [35] as the universal terminator function for chaotic systems with a parametrically long Lanczos plateau. With this terminator function implemented at level k, the operator "environment", i.e., the basis states |O n ) with n > k, become irrelevant for the Krylov complexity. This follows directly from (2.54), in which the sum terminates at n = k, because the square root terminator corresponds to constant ∆ n . In the example discussed in this subsection, this simplification was maximal.

Chebyshev polynomials of the first kind
The second simplest case are the Chebyshev polynomials of the first kind. They correspond to an parameter exceptional value, which is excluded from the Gegenbauer polynomials that will be discussed in the next subsection. The Chebyshev polynomials of the first kind are orthogonal on x ∈ (−1, 1) with the (normalized) weight The monic polynomials are and we have 11 The O(z 2 ) terms give rise to divergent integrals, which is a hint that the expansion should continue with z 2 ln z, but we shall not pursue this issue. 11 For the Chebyshev polynomials of the second kind (Gegenbauer polynomials with β = 1) one has ∆n = 1 4 ∀n.
The wave functions (2.14) are The unitarity relation (2.17) is simply represented by the known sum formula [61,65] To calculate the Krylov complexity we solve again the differential equation (2.45), which in this case reads The solution is which has the asymptotic behaviour K O (τ ) ≈ 2 π τ . We remark that (4.22) establishes the following sum It is also possible to calculate K (2) using Neumann's addition theorem [61], For completeness, we consider also the Laplace transform of the complexity (2.54). As above, two terms contribute to the sum in (2.54), because ∆ 1 − ∆ 0 = 1 2 and ∆ 2 − ∆ 1 = − 1 4 . Therefore, we need the functions Q 0 (z) and Q 1 (z). These are obtained after solving the recursion (2.21), which is quite similar to the case of the previous subsection. Indeed, one simply has R n (z) = R(z) for n ≥ 2 and R 1 (z) = 2R(z) with R(z) given again by (4.10). Hence, so that .
With this information, (2.54) becomes This should be the Laplace transform of (4.22), which is again not obvious to us. To leading order in z, the integral is This value translates into the late-time complexity K O (τ ) ≈ 2 π τ in agreement with the result obtained above.

Gegenbauer (or ultraspherical) polynomials
The Gegenbauer polynomials are orthogonal on (−1, 1) with respect to the even weight function which we have normalized for convenience. The cases β = 0 and β = 1 correspond to the Chebyshev polynomials of the first and second kinds, respectively, and have been discussed in the preceding subsections. Other special cases are the Legendre polynomials (β = 1 2 ). The monic Gegenbauer polynomials are where (a) n = Γ(a + n)/Γ(a) denotes the Pochhammer symbol. For these, we have [61] .
The wave functions (2.14) can be easily obtained from the Gegenbauer formula [61, 10.23.9] where J ν denotes a Bessel function. One finds With the help of sum formulas [61,65], it can be checked that (2.17) holds.
To calculate the Krylov complexity, we need the following sum [66],  It is interesting to note that the limiting case β = − 1 2 gives which is (3.11) for N = 2. This agrees with the fact that the system collapses to the trivial finite system with N = 2, because ∆ 2 = 0 for β = − 1 2 . At late times, the asymptotic expansion of (4.34) gives [69] The function 2F3 has reduced to 1F2 because b = β + 1 2 . For β ≥ 0, the last term can be omitted, and we can rewrite (4.35) as The late-time growth rateK ∞ is not universal. It is a monotonously increasing function of β, starting with zero at β = − 1 2 and reaching unity at β = ∞. However, the larger β the later the linear growth regime is reached (τ 0 grows with β). For − 1 2 < β < 0 the last term in (4.35) gives an oscillating sub-dominant contribution. The complexity (4.34) is plotted in figure 3 for various values of the parameter β.

Hermite polynomials
The last example of classical orthogonal polynomials with for an even measure are the Hermite polynomials. In contrast to the previous examples, the associated spectrum is not bounded. This is, actually, the simplest case and corresponds to the Heisenberg algebra case considered in [38].
We have already met this case among the examples discussed in subsection 3.2.
Hermite polynomials are orthogonal with respect to the (normalized) weight function The wave functions (2.14) are To characterize the operator growth, consider the complexity generating function (2.46), It follows that the complexities (of degree one and two) are

The simplest case with a gapped spectrum
As discussed in subsection 3.1, a gapped spectrum can be associated with a sequence of Lanczos coefficients ∆ n that approach different limits depending on whether n is even or odd. The simplest case, which can be treated analytically [33], is given by For completeness, we first report the solution from [33] and then add some further details in order to prepare the ground for the more general case involving the Jacobi polynomials, which will be the subject of the next subsection.
Our first goal is to obtain the functions c n (z). For this purpose, consider the second equation of (2.21) for even and odd n, . (4.45) This system is solved by 13 (4.47) The functions c n then follow from the first equation of (2.21) and by applying (2.20) recursively, with the result where we have introduced The spectral density can be obtained from c 0 (z) using (2.26). One finds Notice that the delta function at ω = 0 appears only if ∆ e < ∆ o .
Alternatively, one may consider separately the spectral densities for the even and odd polynomials. For the even ones, first use (2.36) to obtain 14 (4.52) From (4.52) we can read off the spectral density 15 The fact that Re and Ro must vanish for large z fixes the sign in the solution of the quadratic equation. 14 We retain the factors of i 2 under the square root in order to keep track of the phase. 15 The comment made in footnote 7 applies here. Since ω 2 ∈ [0, ∞), there is an extra factor of 2 in the term with the delta function. This factor must be dropped, if we one defines x = ω 2 and extends its domain to x ∈ (−∞, ∞).
It is convenient to introduce 54) in terms of which (4.53) reads It is easy to check that (4.55) is normalized on y ∈ (−∞, ∞).
Similarly, to get the measure for the odd polynomials, we calculatē This yields the spectral density These polynomials are monic in ω 2 . One easily obtainsh n = (∆ e ∆ o ) n , which with (2.32) yields The even polynomials follow immediately from (2.35), By consistency, we haveh n = h 2n = ∆ e h 2n−1 = (∆ e ∆ o ) n .
Let us try to find an analytic expression for the complexity. Because we do not know the wave functions φ n (t), but we have their Laplace transforms c n (z) as well as the polynomials P n , we shall look at the generating function K O (λ; z) in the form (2.50). Hence, consider and look separately at the terms with odd and even n. For odd n, using (4.59) and (4.48), the summands become where we recall the definition of T (z) from (4.49). The terms (4.62) can be added up using the definition of the generating function of the Chebyshev polynomials [61]. This results in n odd (4.63) Similarly, for even n, we rewrite the summands in (4.61) using (4.60) and (4.48), which gives, for while one simply gets c 0 (z + iL) for n = 0. Then, summing up the terms with even n yields n even (i e λ ) n √ h n P n (L)c n (z + iL) = c 0 (z + iL) 1 − e 2λ ∆e ∆o T (z + iL) 1 + 2y e 2λ T (z + iL) + e 4λ T (z + iL) 2 . (4.65) Thus, after adding (4.63) and (4.65) one has To make progress, we can use the identities and which follow from the recurrence relations. It is then easy to check that setting λ = 0 gives the expected result Σ(0; z, L) = 1/z after using (4.67), (4.68) and (4.54).
To find the (Laplace transform of) complexity, we need With the help of the identities (4.67) and (4.68) one gets Substituting (4.70) into (4.69), the first term in the brackets does not contribute, because it is odd.
The easiest way to find the integral of the third term is by writing where the contour surrounds only the support of the measure µ(L), and then deforming the contour to pick up the residue at ω = iz 2 . (There is no contribution from ∞.) Finally, the integral of the second term can be written explicitly using the measure (4.50). This yields To find the late-time asymptotics of the complexity, consider the dominant term in (4.72) for small z. Observe that so that the 1/z term precisely cancels between the first two terms on the right hand side of (4.72).
Therefore, the dominant contribution to K O (z) is found by taking the limit z → 0 in the integral, which gives Note that (4.74) Finally, (4.73) translates into the late-time complexity It can be checked that (4.75) satisfies the bound (3.38). Moreover, if we set ∆ e = ∆ o = 1 4 , we reproduce (4.8).

Jacobi polynomials
Jacobi polynomials [61,65] are, in general, related to a non-symmetric measure. They can be used to model a more general system with a gap than that of the previous subsection while maintaining a bounded spectrum. Let us assume a measure with support in ω − < |ω| < ω + plus possibly a delta function at ω = 0. To be concrete, define so that y ∈ (−1, 1).
For large n, the right hand sides of (4.82) and (4.83) reduce toω 2 and 1/(4γ 2 ), respectively, which suggests that the ∆ n with even and odd n approach certain limiting values ∆ e and ∆ o , respectively.
These two values satisfy the relations of the simpler case discussed in the previous subsection (4.51), but we cannot tell yet which one is which.
Cross-over solutions such as those illustrated in figure 4 are models for systems with two different regimes of linear operator growth. Using the results of the previous subsection, if the sequence of ∆ n remains long enough close to the unstable fixed point then, before the cross-over, the early-time complexity is At late times, i.e., after the cross-over to the stable fixed point, the complexity grows as The case ω − = 0 is special. In this case, there is no gap, but there might be a central singularity, and the ∆ n with even and odd n approach the common limit 1 4 ω 2 + from opposite sides with terms of order 1/n. For further details, we refer the reader to [33,55].
The systems discussed in this subsection can be further generalized to spectral densities with delta function singularities at the lower and upper bounds of the spectrum. The orthogonal polynomials for these systems were derived from the Jacobi polynomials in [70]. 16 Recall the convention in footnote 7. Left: α = β = 1/2, ∆ 1 = 0.36, which corresponds to the special case discussed in subsection 4.5.
In this case, the numerical error is sufficient to drive the cross-over to the stable solution. Right:
The even polynomials are given by (2.35). They are orthogonal with respect to the measure where the constant A is determined by imposing that (4.98) be normalized, Here, Φ(α, γ; x) denotes the Kummer confluent hypergeometric function.
First, let us discuss the case without a gap, ω − = 0. In this case, (4.99) simplifies to which restricts ∆ 1 to ∆ 1 ∈ (0, ακ]. Therefore, as mentioned above, given any ∆ 1 ∈ (0, ακ], all other ∆ n are determined recursively by (4.96) and (4.97). As one can easily see from (4.96) and (4.97), the limiting case ∆ 1 = κα gives rise to the following exact sequence ∆ 2n = κn , ∆ 2n+1 = κ(n + α) . Because A = 0 from (4.100) in this case, for this sequence the central central delta function in the measure is absent. We shall calculate the complexity for the sequence (4.101) further below, but before doing so let us return to sequences with 0 < ∆ 1 < ακ.
Another formal solution of (4.96) and (4.97) when ω − = 0 is  Let us return to the special sequence (4.101). In this case, the even polynomials are simplȳ The norm of (4.103) is h 2n = κ 2n n!(α) n . The polynomials P n (ω) (n even or odd) are also called generalized Hermite polynomials [57].
To calculate the complexity, let us consider the generating function K O (λ; z) in the form (2.49).
The even and odd contributions to the sum over n can be considered separately. For the even part, one has where the summation has been performed using the formula [65, 8.976.1], and I α denotes a modified Bessel function. Recall that λ is a formal parameter, and one should take λ < 0 for the sum to be convergent. Similarly, for the odd part one obtains Thus, after adding (4.105) and (4.107), the complexity generating function is To make progress, consider the asymptotic behaviour of the modified Bessel function [65] where we have omitted the higher order terms, which give contributions of order λ 2 in (4.109).
After substituting (4.110) into (4.109) and changing an integration variable by y = w 2 , one realizes that the contributions arising from the terms with e −z in (4.110) are identical to those arising from the terms with e z upon replacing w → −w. Thus, the result of these operations is After another change of variable, w = e λ √ x + √ 1 − e 2λ v, we expand the integrand for small λ and find where the ellipses denote terms of higher order in λ. We remark that the expansion of the integrand also contains a term of order √ −λ, but this term is odd in v and, therefore, vanishes upon integration. We can now perform the integral in v in (4.112), which gives with the complexity where Ψ denotes a confluent hypergeometric function, which is related to the incomplete gamma function.
Instead of using the final result in (4.114), it is more useful to take the first line and rewrite it as This allows to perform the inverse Laplace transform using the general formula (3.5), The integral can be done using [65, 3.952.8], which yields the final result Formally, one can rewrite (4.117) as In order to obtain the late-time complexity, we use the asymptotic behaviour of the confluent hypergeometric function in (4.117). Hence, where the ellipses denote asymptotically vanishing terms. One must distinguish the following three cases: It is evident from figure 6 that, for α < 1 2 , K O (t) < 1 2 κt 2 , whereas α > 1 2 leads to K O (t) > 1 2 κt 2 . Moreover, for large enough α the complexity undergoes some oscillations before approaching the asymptotic late-time behaviour.

Charlier polynomials
Charlier polynomials [57] provide an interesting example of a system with an unbounded discrete spectrum. Allowing for only one parameter, they are the simplest polynomials of the Hahn class.
In the Askey scheme [61], they are at the same level as the Laguerre polynomials.
The monic Charlier polynomials can be defined by their generating function 17  17 We use the monic polynomials of [57], which are related to the polynomials in the notation of [61] by C They are orthogonal with respect to the Poisson distribution, e −a a x x! = n!a n δ mn (5.3) and satisfy the three-term recurrence relation To insert the Charlier polynomials in the context of the recursion method, consider the threeterm recurrence relation for the even polynomials (2.33). Obviously, (5.4) has the form of (2.33) For simplicity, we shall not rescale the energy in this example. Hence, we identifȳ The easiest way to obtain the odd polynomials is to consider their measure,w x = x aw x , which gives In order to calculate the complexity, we shall need the functions of the second kind. Consider Recall that 1 F 1 () is just the Kummer confluent hypergeometric function Φ(). The functionsQ n (z) can by found according to (5.8) by differentiating (5.11) with respect to w, Similarly, for the functions of the second kind associated with the the odd polynomials (5.7), one obtainsQ (5.14) Let us calculate the complexity. Recall from (2.54) that but we also have, from (2.51), After substituting the coefficients (5.5) into (5.15), we use (5.16) and (5.17) to eliminate all the terms involving the even polynomials. This yields where we have substituted the relevant relations from subsection 2.3 and abbreviated First, consider the combination After writing out the generalized hypergeometric series, we rearrange the sums over n and k, A further sum rearrangement introducing m = n + l gives The sum in (5.21) results in .
At this point, using x − y ± = z 2 ∓ 2iz √ x from (5.19), we calculate the combination which appears in (5.18). After substituting this result, (5.18) becomes , (5.25) where the term arising from the third term in the brackets of (5.24) has been combined with the first one after shifting the summation index.
It is possible to carry out the summation in (5.25), but the result is not very enlightening. It is more instructive to rewrite Finally, it is straightforward to compute the inverse Laplace transform of (5.26), The complexity (5.27) is shown in figure 7 for a number of values of the parameter a. The behaviour of the complexity at late times is quite surprising. Naively, because the ∆ 2n grow linearly and the ∆ 2n+1 are constant, from the results of subsections 4.1 and (4.4) one would expect a growth between linear and quadratic. In fact, the complexity does not grow at late times, but oscillates more or less erratically (depending on a) around the mean value which is found from (5.27) by replacing the sine squares by 1 2 . The same value can be obtained by looking at the coefficient of the 1/z term in (5.26).  The initial oscillations are easily understood when a is large. In that case, the Poisson weight e −a a x+1 /(x + 1)! can be approximated by a Gaussian distribution peaked at x = a − 1. Hence, the dominant oscillation is with all the oscillations with nearby frequencies starting in sync at τ = 0. As time passes, they get more and more out of sync, which produces the late-time average (5.28).
A complexity that approaches a constant is characteristic of a finite chain of wave functions.
Indeed, the Charlier polynomials can be obtained as a limit from the Krawtchouk polynomials, which have a finite spectrum. The limit is [61] lim N →∞ We will discuss the Krawtchouk polynomials in the next subsection.

Krawtchouk polynomials
The Krawtchouk polynomials [61] provide an example of a finite spectrum. They are a special case of the Hahn polynomials, and the limit (5.29) gives the Charlier polynomials. We shall work with the monic Krawtchouk polynomials, which are given explicitly by 18 and are orthogonal with respect to the binomial distribution, Moreover, they satisfy the three-term recurrence relation n−1 (x) . (5.33) 18 In the notation of [61], the Krawtchouk polynomials are Kn(x; p, N ) = 2F1 −n, −x; −N ; p −1 .
We have already met the Krawtchouk polynomials in the special case p = 1 2 of the simple finite system (3.12). Here, we shall provide a different example. For this purpose, we identify (5.33) with (2.33) by setting with x, n ∈ {0, 1, . . . , N }. Therefore, for the even polynomials, we havē For the odd polynomials, one easily gets from (2.34) The functions of the second kind can be calculated using the generating function (5.31).
Using (5.10), one obtains and (5.37) givesQ Similarly, for the odd part, one finds We have now all the ingredients that we need to calculate the complexity. We shall proceed as in the previous section, cf. (5.15)-(5.18). The expression corresponding to (5.18) is where y ± are again given by (5.19).
Consider the combination After writing out the generalized hypergeometric series, this becomes We rearrange first the sums over n and k, and then the sums over n and l, setting m = n + l, The sum over l can be performed using (5.22), which yields At this point, the sum over k can be extended to infinity, because x ∈ {1, 2, . . . , N } in (5.42), so that (1 − x) k = 0 for k ≥ N . The summation can then be carried out, with the result It is a nice check that (5.44) reduces to (5.23) in the limit pN = a, N → ∞.
To continue, we use x − y ∓ = z 2 ± 2iz √ x from (5.19) and calculate the combination We now return to (5.42), substituting (5.45) and shifting x by one in the term arising from the third term in brackets in (5.45). This gives where we have carried out similar steps as in the previous subsection. We can readily perform the inverse Laplace transform of (5.47). This yields our final result for the complexity, At late times, the complexity oscillates around the mean value 49) which is found by setting the sine squares in (5.48) to 1 2 . We illustrate the complexity for various values of N and p in figure 8. It can be seen that the complexity is qualitatively quite similar to the case with the Charlier polynomials discussed in the previous subsection.

Meixner polynomials
The Meixner polynomials [61] provide another example of an unbounded discrete spectrum. In the Askey scheme, possessing two parameters, they are at the same level as the Krawtchouk polynomials. Therefore, they are a special case of the Hahn polynomials, and reduce to the Charlier polynomials in a certain limit.
As before, we shall work with the monic polynomials, which are given explicitly by 19 with β > 0 and 0 < c < 1.
They may be defined via their generating function Moreover, they satisfy the three-term recurrence relation In order to use the Meixner polynomials in the context of Krylov complexity, we identify (5.53) with the recurrence relation of the even polynomials (2.33) by setting with x, n ∈ {0, 1, 2 . . .}. Therefore, for the even polynomials, one has For the odd polynomials, one gets from (2.34) In order to find the functions of the second kind for the even part, we consider again the generating function Using (5.10), one finds from which one gets the functions of the second kind Similarly, for the odd part, one has To calcultate the complexity, we proceed again as in subsection 5.1, cf. (5.15)- (5.18), with the where y ± are again given by (5.19).
Consider the combination The next steps are analogous to the calculations in the previous subsections. We write out the generalized hypergeometric series, and rearrange first the sums over n and k, Rearranging also the sums over n and l, with m = n + l, yields The sum over l can be performed using (5.22), so that we get Carrying out the remaining summation, we obtain One can verify that (5.65) reduces to (5.23) in the limit (5.54).
To continue, we use x − y ∓ = z 2 ± 2iz √ x from (5.19) and calculate the combination Returning to (5.63), we substitute (5.66). Shifting also x by one in the term arising from the third term in brackets in (5.66), we find where we have carried out similar steps as in the previous subsections.
Finally, after performing the inverse Laplace transform, we obtain the complexity  The result (5.68) is quite similar to the results of the Charlier and Krawtchouk cases. At late times, the complexity oscillates (more or less erratically) around the mean value which is found by setting the sine squares in (5.68) to 1 2 . The complexity is plotted in figure 9 for various values of β and c.

Other polynomials of the Hahn class
The remaining polynomials of the Hahn class [61] are the Hahn polynomials, the Meixner-Pollaczek polynomials, and the continuous Hahn polynomials. We will not provide the full analysis of these cases, but limit ourselves to a few comments.
This is a generalization of the Krawtchouk and Meixner cases. The former is reproduced in the limit α = pt, β = (1 − p)t, t → ∞, the latter by The Meixner-Pollaczek polynomials are orthogonal with respect to a weight function with support on (−∞, ∞). They have two parameters, which places them at the same level as the Jacobi, Krawtchouk and Meixner polynomials in the Askey scheme [61]. In our context, one parameter is fixed by demanding that the spectrum is symmetric (φ = π 2 ). The resulting polynomials give rise to the SU (1, 1) case of [38], as discussed in subsection 3.2.
A more general case with continuous support, x ∈ (−∞, ∞), is given by the continuous Hahn polynomials. These polynomials occur in relation with exactly solvable Quantum mechanics problems [64,71,72]. In general, continuous Hahn polynomials are parameterized by two complex coefficients a and b, or, equivalently, by four real parameters. However, only three of these real parameters are free, because one parameter can be absorbed by a shift of the independent variable x.
Moreover, demanding the weight function to be symmetric reduces the number of free parameters to two. There remain two possibilities, in both of which P n (x) ∼ p n (x; a, b,ā,b) is a continuous Hahn polynomial of degree n. First, one can choose a and b to be real and positive, in which case one has the coefficients The special case b = a + 1 2 reduces to the Meixner-Pollaczek case. Second, for b =ā, a +ā > 0, one has , (a +ā > 0) .
In both cases, (5.72) and (5.73), ∆ n grows as n 2 for large n. This implies maximal exponential growth of the complexity at late times, just as in the Meixner-Pollaczek case.

Tricomi-Carlitz polynomials
As an exotic example of non-classical orthogonal polynomials, we consider the Tricomi-Carlitz polynomials [57,73,74]. They are also known as Karlin-McGregor polynomials, as they have been studied by Karlin and McGregor in the connection with random walks [75]. In our context, the Tricomi-Carlitz polynomials are interesting, because their sequence of ∆ n converges to zero for n → ∞. By the argument given in subsection 3.1, the spectrum must consist of countably many points, i.e., it must be discrete. As we shall see below, the spectrum is bounded, but countably infinite.
Tricomi [73] studied the polynomials n denotes a Laguerre polynomial. Carlitz [74] noted that, if one sets with α > 0, then the polynomials f satisfy the three-term recurrence relation Because the sequence of ∆ n = n (n + α)(n + α − 1) (n ≥ 1) (6.5) approaches zero for large n, by the argument given in subsection 3.1, the spectrum must be discrete, as we mentioned above. Indeed, the orthogonality relation is where Ψ (α) is a step function with jumps The generating function is Furthermore, the Tricomi-Carlitz polynomials are related to the Hermite polynomials through the limit [76] lim α→∞ f (α) n √ 2x α = 2 −n/2 n! H n (x) . (6.9) From (6.5) we see that ∆ n ∼ n −1 for large n, which corresponds to the exponent λ = −1 in the continuum approximation of subsection 3.3. This would indicate a late time complexity that grows as t 2/3 , but as we pointed out in that subsection, this result should be taken with a grain of salt, because the continuum approximation is not justified at large times, if λ < 0.
Let us find the functions of the second kind (2.23) for the Tricomi-Carlitz polynomials. This time, our approach is to calculate Q 0 (z) and then use the recurrence relation (2.27) to obtain the and little more to prove that satisfy the recurrence relation (2.27). Written as a series, (6.15) reads Let us attempt to calculate the complexity. We shall use the generating function (2.50), which now reads and consider first the combination and Q n from (6.6), (6.3), and (6.16), respectively, gives After rearranging the summations, this becomes (6.20) The sum over n can be carried out now, which yields To proceed, let us define the following three sums, It is easy to see that they satisfy It is a bit more difficult to show that holds. To prove (6.24), take S 1 , substitute the identity and write m + α = (α + l − y −2 ) + y −2 + (m − l). For the first term, the sum over l is non-zero only for m = 0, which gives the first term on the right hand side of (6.24). The other two terms yield the second term on the right hand side of (6.24) after some manipulations.
After solving the system (6.23)-(6.24) one obtains (F Q)(x, y) = e λ S 1 + S 2 = 1 y + x e λ 1 + (e 2λ −1)S 3 . (6.26) Now, let us return to the generating function (6.17), which is While setting λ = 0 simply gives the completeness relation, to find the complexity we are interested in the term linear in λ. In (6.26), this term is (6.28) The double sum in (6.28) is known as a (Humbert) confluent hypergeometric series of two variables [77], [65, 9.261], so which has the following integral representation [65, 3.385] or [78], We were unable to make further progress analytically. One should carry out the summation over the spectrum in (6.27), setting y = iz − x. Obviously, odd terms in x cancel in this sum, in particular, the terms that would behave as 1/z 2 . One such term is the first term on the right hand side of (6.30), but also the Humbert series hides such a term. It comes from (α − y −2 ) m+2 in the denominator in (6.28), which has a 1/z contribution that is odd in x for x = ±x k = (α + k) −1/2 .
From the absence of 1/z 2 terms we can conclude that the late time complexity grows less that linearly in time, but in order to conclude that it approaches a constant we should know whether the resulting expression is analytic in z. We suspect that it is not. The issue is complicated by the fact that it is not possible to consider the small z limit in (6.29) before summing over the spectrum, because the spectrum contains (countably) infinitely many values of |x| that are smaller than any given, even if very small, |z|. We mention that the continuum approximation, which cannot be considered as valid in this case, would indicate a non-analytic behaviour ∼ z −5/3 . There is no evidence of such a behaviour in our results.

Operator inner products
In this final section, we shall return to an important detail, which was left out at the beginning of subsection 2.1. Given an operator O in a system with Hamiltonian H, its complexity K O (t) will crucially depend on the choice of the inner product adopted before applying the Lanczos method. Other than that, there is no additional freedom or parameter to choose, except possibly for approximations such as a terminator function [33]. The choice of the inner product has a direct bearing on the results of the entire procedure, i.e., the Lanczos coefficients, the measure µ(L), the system of orthogonal polynomials, and the complexity. It is known that different choices may even lead to essentially different qualitative behaviours of the complexity including its late-time behaviour. We shall limit ourselves to a brief review of physical inner products that have been introduced in the literature.

Frobenius inner product
If the physical Hilbert space is finite, dim(H) = D, then the easiest choice for the inner product between two operators A and B is the Frobenius (or Hilbert-Schmidt) inner product The dimension ofĤ, which is the space of linear operators acting on H, is D 2 , so the Lanczos algorithm necessarily terminates. The maximal dimension of the Krylov subspace K = span{O n } is somewhat smaller than D 2 , namely dim(K) = N ≤ D 2 − D + 1 [29].
For infinite-dimensional Hilbert spaces, besides writing D = tr [1], one would need to define the trace with a suitable regularization. For example, on may consider (7.1) as the infinite temperature limit of the thermal inner product, which we will discuss in the next subsection.
In some references, the "connected" inner product is used instead of (7.1), i.e.,

Thermal two-point functions
In quantum thermodynamics, for a system with Hamiltonian H and inverse temperature β, the thermal expectation value of an observable A is defined by where Z = tr e −βH is the partition function.
In [33] the connected version of (7.4) is used, in which A † β B β is subtracted from the right hand side. However, our remark below (7.2) applies here as well.

Microcanonical inner product
In systems with an unbounded spectrum, the Krylov complexity calculated starting from the canonical inner product (7.4) often obscures the telltale signatures of quantum chaos. In order to improve on this situation, a microcanonical inner product has been introduced in [35]. The main argument of [35] is that the mean energy between two energy eigenstates, E = 1 2 (E i + E j ), is a conserved quantity under the Liouvillian, and operator growth should be considered separately in each conserved sector. In what follows, we shall review the construction of the microcanonical inner product of [35] with some generalizations. For simplicity, we will consider a system with discrete energy spectrum in the thermodynamic limit and take the ground state energy to be zero.
If we define the density of states where the sum is over all energy eigenstates, then the thermal expectation value (7.3) takes the classical form where A E denotes the expectation value of A in the microcanonical ensemble of states with energy E, which we call M(E). 21 Specifically, 1 .
(7.8) 21 In the sense of the thermodynamic limit, the definition of the microcanonical ensemble includes an energy window 1, i.e., M(E) = {|i : |Ei − E| < }, and one should accordingly smear the delta functions in (7.6).
Notice that our normalization is such that 1 β = 1 E = 1.
The idea of [35] is to rewrite the thermal correlator (7.4) as a thermal expectation value, i.e., in the form where A|B E is the microcanonical inner product we are looking for. To do this, one starts from and introduces the pair density 22 Then, (7.4) takes the form (7.9) with the microcanonical inner product dω ρ(E, ω)γ β (ω) A|B E,ω , (7.13) where i|A|j * i|B|j . (7.14) For the Wightman inner product, for which γ β (ω) = 1, (7.13) reproduces the simple expression obtained in [35]. In general, both ρ(E, ω) and γ β (ω) are even functions of ω.
In subsection 2.1 we assumed the start operator |O 0 ) of the Lanczos algorithm to be normalized.
In the canonical or the infinite-temperature settings, the initial norm is irrelevant, because one can normalize the operator once and for all. In the microcanonical approach, the norm is energydependent and one should keep track of it. Hence, we write explicitly Clearly, the microcanonical inner product gives rise to a measure with bounded support.
The wave functions (2.14) become φ E,n (t) = i n O|O E h n O(t)|O n E . (7.18) In particular, with (7.15), Thermal Krylov complexity can then be defined as a weighted average of the microcanonical complexities K O,E (t), which are defined as usual by (2.41) with the wave functions φ E,n (t). We point out that the weight is not unique. Motivated by (7.19), the authors of [35] chose to include the energy dependent operator norm into the weight, We remark that this is a choice of the unit length on the (energy dependent) chains of wave functions rather than a consequence of operator normalization.
Because the measure (7.17) has bounded support with |ω| ≤ 2E, in the thermodynamic limit the Krylov complexities K O,E (t) grow at most as K O,E (t) Et at late times. Therefore, (7.21) gives immediately which resembles Lloyd's bound on complexity growth [4].