Krylov complexity and chaos in quantum mechanics

Recently, Krylov complexity was proposed as a measure of complexity and chaoticity of quantum systems. We consider the stadium billiard as a typical example of the quantum mechanical system obtained by quantizing a classically chaotic system, and numerically evaluate Krylov complexity for operators and states. Despite no exponential growth of the Krylov complexity, we find a clear correlation between variances of Lanczos coefficients and classical Lyapunov exponents, and also a correlation with the statistical distribution of adjacent spacings of the quantum energy levels. This shows that the variances of Lanczos coefficients can be a measure of quantum chaos. The universality of the result is supported by our similar analysis of Sinai billiards. Our work provides a firm bridge between Krylov complexity and classical/quantum chaos.


Introduction
Since quantum mechanics is the basis of the description of nature, studying chaos in the quantum realm, or quantum chaos is important both theoretically and experimentally.In recent years, quantum chaos has attracted interest in a variety of studies, including thermalization processes in thermodynamic systems and black hole dynamics.However, a clear standard for the definition of quantum chaos is still overdue.Historically, there have been several proposals to define quantum chaos.The traditional characterization is by the statistical distribution of adjacent energy level spacing in quantum systems [1].Turning to the dynamics of quantum systems, it is expected that local operators in chaotic systems will become more complex with time evolution.To measure this operator growth, the outof-time-order correlator (OTOC) [2] has attracted much attention as a possible indicator of quantum chaos.
The exponential growth rate of the OTOC, called the quantum Lyapunov exponent, is regarded as a naive quantum counterpart of the Lyapunov exponent in classical systems.Remarkably, the quantum Lyapunov exponent of a finite-temperature quantum manybody system has an upper bound determined by the temperature [3].This upper bound is expected to be an indicator of the existence of a gravity dual through the AdS/CFT correspondence.The Sachdev-Ye-Kitaev (SYK) model [4,5], which saturates the upper bound, has played important roles in the studies of quantum chaos and quantum gravity.Furthermore, the thermodynamic well-definedness of the OTOC suggests an upper bound on the energy dependence of the Lyapunov exponent for very general systems [6].
However, it has been shown that even in a classically chaotic system, the OTOC does not necessarily show exponential growth [7], 1 and conversely, even in a classically nonchaotic system, the OTOC can grow exponentially due to the presence of unstable points in the potential [8][9][10].Thus, the OTOC alone is unfortunately insufficient to characterize the quantum chaoticity.Also, the relationship between the energy level statistics, which has been previously proposed as a characterization of quantum chaos, and the OTOC is unclear.
Recently, the notion of Krylov complexity for operators, Krylov operator complexity, 2 was proposed as a new indicator to evaluate operator growth more directly and quantitatively [11]. 3Also, based on the same idea as the Krylov operator complexity, the Krylov complexity for states, Krylov state complexity, has been proposed as a quantitative measure of the complexity of a quantum state in the Schrödinger picture [13]. 4In the calculation of the Krylov complexity for operators and states, a given quantum system is reduced to a one-dimensional chain model, where the hopping is given by a series of real numbers called Lanczos coefficients.Thus the time evolution of the Krylov operator/state complexity is encoded into the set of Lanczos coefficients.
A natural question is a possible relationship between the Lanczos coefficients and chaos.In this regard, in finite-dimensional quantum systems such as spin systems and the SYK model, novel correlations between certain variances of the Lanczos coefficients, asymptotic values of the Krylov operator complexity and the energy level statistics were reported [16,17]. 5Unfortunately, these models have no classical counterparts, and their relation to classical chaoticity is not clear.In fact, if we are to define quantum chaos consistently with classical chaos, the first interest to us is the quantum mechanical system obtained by quantizing the classical chaotic system.In this paper, we consider billiard systems [19][20][21] which are one of the most-studied classical/quantum system of chaos.We investigate in the billiard systems the relationship between the classical chaoticity, the quantum energy level spacing statistics, and the Lanczos coefficients for operators and quantum states and show the existence of a correlation between them by numerical analyses.This result on the Krylov operator/state complexity and their Lanczos coefficients may also hold for broader general quantum mechanical systems, providing a key step for refining the definitions of quantum chaos.
The novel venue on the bridge connecting complexity and chaos can shed light more on quantum gravity.Complexity has attracted much attention in black hole studies, and specifically, there are several conjectures in the AdS/CFT correspondence that the complexity of the holographic quantum system probes the interior of the dual black hole [22][23][24][25][26][27].One of the issues there is the existence of several different proposals on the definition of complexity, allowing ambiguity.In contrast, the Krylov operator/state complexity is defined unambiguously.It is an interesting and important question whether the Krylov operator/state complexity describes the interior of the black hole. 6his paper is organized as follows.In Sec. 2, we first review the Lanczos coefficients, the Krylov operator/state complexity, and their possible relation with chaos.Then, we illustrate the Krylov complexity by showing some analytically calculable complexity in integrable quantum mechanical systems.In Sec. 3, we consider the stadium billiard, which is a typical chaotic system, and numerically investigate the Krylov complexity for the momentum operator.We see that a certain variance of the Lanczos coefficient is correlated with the chaoticity of the system.In Sec. 4, we numerically investigate the Krylov state complexity in the stadium billiard.Again, we see that there is a correlation between the variance of the Lanczos coefficients and chaos.In Sec. 5, we perform the same analysis for the Sinai billiard as for the stadium billiard.Sec.6 is for our summary and discussion.In App.A, we discuss the details of numerical calculations.

Review on Krylov complexity
Here, we review the definitions of Krylov operator/state complexity and their possible relation to chaoticity.
where L is the Liouvillian superoperator, we can construct an orthonormal basis for H O by the following procedure, which is known as the Lanczos algorithm: (2.4) 7 See [33] for a pedagogical introduction to Krylov operator complexity and some analytical methods.
For studies of Euclidean time evolution and Krylov complexity, see [34,35]. 8One can choose another definition of the inner product.Although there is often a normalization coefficient in (2.3) in the literature, we omitted it since such an overall constant does not change the Lanczos coefficients.
Note that φ n (t) obeys the normalization condition, (2.5) Substituting (2.4) into the Heisenberg equation yields where the dot represents the derivative with respect to time.The initial condition is φ n (0) = δ n0 ∥O∥ by definition.It is convenient to normalize the operator in the first place so that the condition simplifies to φ n (0) = δ n0 .The Krylov complexity for the operator O is defined as (2.7) The operator O n contains the nested commutator L n O, which is expected to become more complex with n.Therefore, the Krylov operator complexity roughly measures the number of the nested commutators in the given Heisenberg operator O(t).Equation (2.6) can be regarded as a hopping model on a one-dimensional chain.The Lanczos coefficients {b n } represent the hopping amplitudes.The Krylov operator complexity can be interpreted as the expectation value of the "position" n of the hopping particle on the one-dimensional chain.
Originally, the authors of [11] considered many-body systems in the thermodynamic limit.They conjectured that the Lanczos coefficient b n grows at most linearly with n, and when the system is chaotic, b n asymptotically grows linearly as b n ∼ αn + γ (n → ∞) with some constants α and γ.It was also shown that when the Lanczos coefficient behaves as b n ∼ αn + γ, the Krylov operator complexity grows exponentially as K O (t) ∼ e 2αt .Thus, the exponential growth of the Krylov operator complexity may be an indicator of chaos.Note that this argument and conjecture assume the thermodynamic limit, and do not necessarily hold in finite-dimensional systems. 9n finite-dimensional systems, the Lanczos algorithm must terminate.Therefore, even if the Lanczos coefficient grows at first, it necessarily turns to decrease at some point [37].The authors of [16] proposed a possible relation between the Lanczos coefficients and the chaoticity.Given an operator O(t), define the moments {µ n } of the two-point correlation function (O(t)|O(0)) by It is known [11,16,38] that the Lanczos coefficients {b n } of O(t) are intimately related to the moments {µ n } via where D n is defined as The authors of [16] argued that when the system is integrable, that is, the distribution of the adjacent energy level spacing obeys the Poisson statistics, there may be more possibility for D n to become small, which is in turn expected to cause fluctuation in the Lanczos coefficient b n through (2.9).Conversely, if the system is chaotic, that is, the level spacing obeys the Wigner-Dyson statistics, the behavior of b n is expected to become less erratic.The authors of [16] then introduced the variance of the Lanczos coefficient as where ⟨• • • ⟩ represents the mean value. 10This quantity measures the magnitude of the erratic behavior of b n .They confirmed the expected behavior in the XXZ spin chain model and in the SYK model.Adding an integrability breaking term in the XXZ model [17], they also investigated the correlation between the variance σ 2 and the average ⟨r⟩ of the parameters which is a popular characterization of the level statistics [39,40].Specifically, ⟨r⟩ takes the following values depending on the distribution of adjacent level spacing in the spectrum under consideration [40]: where GOE, GUE and GSE are ensembles of the Hermitian random matrices model (see [40] for their definitions).It is known that the statistical distribution of adjacent energy level spacing of stadium and Sinal billiards obeys that of GOE.If the dimension of the Krylov space is finite, the Krylov operator complexity remains at a finite value and often saturates at late time.It is proposed that this late-time saturation value may be concerned with the erratic or non-erratic behavior of the Lanczos coefficients mentioned above [16].We will discuss this in Sec. 7.

Krylov state complexity
Similarly to the Krylov operator complexity, we can define Krylov complexity for quantum state [13] and substituting (2.15) into the Schrödinger equation, we have The initial condition is ψ n (0) = δ n0 by definition.The Krylov state complexity of the state |ψ⟩ is defined as (2.17) In [13], the Krylov state complexity for the thermofield double (TFD) states in some chaotic systems was considered and found to show the characteristic peak and plateau structure, which [13] presumed was universal for the chaotic systems.Recently, this structure is studied both analytically and numerically in [12].We will investigate this feature by considering the Krylov state complexity for states other than the TFD in chaotic quantum mechanics.
Unlike the case of Krylov operator complexity, there are two types of Lanczos coefficients, a n and b n .We introduce variances of Lanczos coefficients a n and b n as i ≡ ln It is known that these Lanczos coefficients can be expressed in the same manner as (2.9) [38].
We will later numerically study the variances similar to (2.11) and investigate the possible relation with chaoticity.

Analytic examples of quantum mechanical complexity
Here, as an illustration, we present some explicit examples of the analytic calculation of the time evolution of Krylov complexity for operators and states in quantum mechanics.

Krylov operator complexity and inverse harmonic oscillator
As a typical operator we take a normalized Gaussian operator Here x is the position operator, and α is a real positive constant.The reason why we take this Gaussian operator is that it is normalizable and thus well-defined: (O 0 |O 0 ) = 1 where the inner-product is the standard trace in the quantum mechanics Hilbert space, with dx in the x-representation.
In this one-dimensional quantum mechanical system, the time evolution of the Krylov complexity for the Gaussian operator (2.19) is quite nontrivial, as the commutator of the operator with a given Hamiltonian generically generates infinite kinds of operators consisting of x and p.As a calculable example we consider the following Hamiltonian where e is a real constant.Taking the commutator of the Gaussian operator O 0 with the Hamiltonian n times generates the following form of the operator up to a constant normalization.Since the Krylov basis is given as an orthonormal set of operators, we conclude that such a set of the form above is unique and given by the Hermite polynomials, In other words, the Krylov basis in this case is identical to that of the energy eigenfunctions of a harmonic oscillator.Using this Krylov basis, we obtain the Lanczos coefficient as Therefore at the large n the Lanczos coefficient grows linearly in n.
For this form of the Lanczos coefficient, the time evolution of the Krylov complexity is explicitly obtained.Applying the formula given in [11], we find This grows exponentially for t ≫ √ α/e.In [11], it was argued that the exponential behavior is related to chaos.However, the Hamiltonian (2.20) is too simple that it is integrable and does not give any chaos.Then where can we locate a possible relation between the exponential growth of the complexity and the dynamical behavior of the current system?In fact, we can spot a behavior typical to chaotic systems: it is an exponential run-away behavior of the trajectory.The classical Hamilton equation is whose solution is exponentially growing in time, Thus the classical motion is exponential in time.The origin of this exponential behavior can be more easily understood if one makes a canonical transformation11 The transformed Hamiltonian is which is the inverse harmonic oscillator, giving the exponential growth in the motion.This inverse harmonic oscillator is not chaotic while has been used for the demonstration of scrambling in quantum models [9,10].This suggests the exponential behavior of the Krylov complexity of the present integrable system.

Krylov state complexity for a free particle and a harmonic oscillator
As another illustrating analytic examples, we calculate the Lanczos coefficients of the Krylov state complexity of Gaussian quantum mechanical states for the case of a single free particle and the case of a harmonic oscillator. 12We will see that they grow linearly in n, as a n ∼ 2b n ∝ n at a large n.We also find that the complexity time evolution for the free particle case goes as t 2 .
Let us start with the free particle in one spatial dimension (the case of the harmonic oscillator is analyzed precisely in the same manner and will be mentioned later).The Hamiltonian is where p is the momentum operator, and we consider the wave function at t = 0 as where α is a real positive number to parameterize the width of the initial wave function in p space.The Lanczos algorithm for the Krylov state complexity [13] is applied to this initial state.In the following we prove that the algorithm generates Hermite polynomials.
The algorithm basically consists of the following procedure: first, multiply the Hamiltonian to the n-th state |K n ⟩, and second, subtract a linear combination of |K n ⟩ and |K n−1 ⟩ to make the generated state be orthogonal to |K n ⟩ and |K n−1 ⟩, and third, normalize the generated state.This procedure automatically makes sure that generated states are orthogonal to each other and normalized.It is easy to notice here that in fact if we start with the Gaussian state (2.30) it automatically generates states represented by Hermite polynomials H 2n times the Gaussian.The reason is that the Hamiltonian is simply p 2 so the n-th state |K n ⟩ need to be of the form up to an overall factor.The only mutually orthogonal set of functions of this form is the Hermite polynomials, as is well known for energy eigenstates of a harmonic oscillator.A simple calculation leads to the result With this explicit Krylov basis, we can easily compute the Lanczos coefficients.
At a large n, we find n.
Therefore the Lanczos coefficients grows linearly in n and satisfies a relation a n ≃ 2b n .
As another example, let us consider the harmonic oscillator with the Hamiltonian We find that the logic for the free particle also applies to this harmonic oscillator precisely, and the Krylov basis is just again given exactly by (2.32).The reason is that the effect of the harmonic oscilator potential x 2 on any function of the form (2.32) is exactly the same as the multiplication of p 2 with some addition of a constant.The Lanczos coefficients are obtained in a similar calculation but with the harmonic oscillator Hamiltonian, At a large n, we again find the linear growth of the Lanczos coefficients as Note that for α = 1/ √ 4mω all the Lanczos coefficients disappear.This is because the initial wave function becomes the ground state wave function for that value of α and thus the Krylov basis consists just only of the ground state wave function.
Finally let us calculate the time evolution of the Krylov state complexity for the case of the free particle.The time evolution of the wave function is Then we find the Krylov components of this wave function as (2.39) Putting t ≡ t/(4mα), the explicit integration gives (2.40) Using this, the Krylov state complexity is calculated as 13 (2.43) Thus the Krylov state complexity of the Gaussian state for a free particle evolves as t 2 , with the typical time scale mα where m is the mass of the particle and √ α is the spatial Gaussian width.In summary, the Lanczos coefficients are linearly growing, while the system is not chaotic, and the Krylov state complexity is not exponentially growing in time.
In this subsection we have performed analytic evaluation of the Krylov complexities.Unfortunately, systems in which the Krylov complexity can be analytically calculated are quite limited, and in particular, chaotic systems allow only numerical evaluations.In the next section, we explain our numerical method.

Preliminary for numerical calculations
We consider numerical calculations of Krylov complexity.Since the dimension of the Hilbert space is often infinite, some regularization is necessary to perform numerical calculations.In our later calculations, we consider only a finite number of levels and ignore the others.For a Hamiltonian of the form we perform our numerical calculation by the following procedure. 14We assume that the system is bounded from below.Generalization to systems with more degrees of freedom is straightforward.

Krylov operator complexity
Suppose that we want to calculate the Krylov complexity for the momentum operator p 1 . 15or that purpose, we first numerically solve the Schrödinger equation 2) 13 The nontrivial summation in (2.43) can be performed in the following manner.First one should note that the consistency of the probability conservation 1 = n ⟨Kn|ψ(t)⟩ 2 is proven by the expansion and the identification y ≡ ( t) 2 /(1 + ( t) 2 ).Then the desired summation (2.43) is rephrased in terms of y as (2.42) We define energy eigenfunctions as ϕ n (x, y) ≡ ⟨x, y|n⟩ (n = 1, 2, • • • ).We compute them by numerical calculations. 16The matrix representation of the position operator x is given by We can evaluate above by the numerical integration.Although we can also obtain the matrix representation of the momentum operator using derivatives of eigenfunctions, the derivative of the numerical solution tends to lose numerical accuracy.Thus, we employ the technique used in [7], as described below.With (3.1), the momentum operator p 1 can be expressed as Then using (3.4), we can find the matrix element of p 1 as where We now define the truncated momentum operator as where N max is an truncation number.Then, we perform the Lanczos algorithm for P .Technically, in the numerical calculation, the naive Lanczos algorithm does not work well near the end of the algorithm due to numerical errors.In our numerical calculation, we use the full orthogonalization method [42,43] instead of the naive Lanczos algorithm to ensure orthogonality.In the current case, the algorithm is as follows: Let K P be the dimension of the Krylov space, which is determined by b K P = 0.
where N max is an arbitrarily chosen truncation number.Each component is obtained by the integral where ψ(x, y) ≡ ⟨x, y|ψ⟩.In our numerical calculation, however, we simply use the following vector as the initial state: while the other initial states give qualitatively similar numerical results.Then, we perform the Lanczos algorithm for the state Ψ.As we did for the operator case, we need to modify the naive algorithm to avoid numerical errors during the orthogonalization.Defining the algorithm goes as follows [43]: , and go to step 3. Let K Ψ be the dimension of the Krylov space, which is determined by b K Ψ = 0.Then, we numerically solve (2.16) and obtain the Krylov state complexity (2.17).
We set N max = 500 in numerical calculations of Krylov state complexity shown in the following sections.The variance in (2.18) is calculated for a n , b n with 50 ≤ n ≤ 250, where this range is chosen to avoid transient feature at small n and possible numerical error at large n ∼ N max .

Krylov operator complexity in the stadium billiard
It is well known that the stadium (circular) billiard is chaotic (integrable) both classically and quantum mechanically [44][45][46].In Fig. 1, we show the geometry of the stadium billiard.The shape of the billiard is determined by a/R.Considering the one-parameter deformation of the billiard, we study the integrable/chaos transition and correlation between the variance of Lanczos coefficients and classical/quantum chaoticity.
Classical chaoticity can be measured by the Lyapunov exponent.In Fig. 2(a), we show the a/R dependence of the Lyapunov exponent, which is consistent with [44].In the calculation, we set the area of the billiard and the velocity of the particle to the unity.When a/R = 0, the system becomes an integrable circular billiard.As a/R increases, so does the Lyapunov exponent and the system becomes chaotic.
Quantum chaoticity is traditionally distinguished by the level statistics, or more con- veniently by the ratio (2.12).In Fig. 2(b), we show the a/R dependence of the ratio. 19he orange line corresponds to the value for the Poisson statistics and the green line corresponds to that for the Wigner-Dyson (GOE) statistics (2.13).We can see that the transition similar to the classical case takes place as the a/R increases, which is consistent with [45].

The early time dependence of Krylov operator complexity
Using the method in Sec. 3, we numerically compute Krylov operator complexity for the truncated momentum operator P defined in Eq. (3.7), with truncation N max = 100. 20In Fig. 3, we show the Lanczos coefficients for a/R = 0 and a/R = 1.Although the former case is integrable and the latter is chaotic, the initial behaviors of the corresponding Lanczos coefficients are almost identical.We identify the dimension of the Krylov space K P as K P = 9900.Note that the horizontal axis in Fig. 3 is in log scale.In Fig. 4, we show the Krylov operator complexities as functions of t for the stadium billiards with a/R = 0, 0.1, 0.2, • • • , 1.While the stadium billiard with a/R > 0 is chaotic, the early time growth of the Krylov operator complexity is not apparently exponential. 21igure 4 shows that C(t) saturates by t ≲ 30 and reaches asymptotic values that depend on a/R.In Fig. 5, we show the dependence of saturation value of C(t) (the average of C(t) taken over 40 ≤ t ≤ 100).We discuss its implications in Sec. 7.

Correlation between Lanczos coefficients and chaos
The Lanczos coefficients for a/R = 0 distribute broader compared to a/R = 1 in Fig. 3.In Fig. 6 we show the variance of Lanczos coefficients (2.11) as a function of a/R.The variance becomes larger in the integrable regime compared to the chaotic regime.To compare with other chaos indicators, we show scatter plots in Fig. 7, where the points are sampled from 0 ≤ a/R ≤ 0.5. 22We can see that there are correlations between λ, σ 2 , and ⟨r⟩.Quantitatively, a correlation between two sets of data A and B can be evaluated by the correlation coefficient defined by where E[ • ] means the average value.If there is no correlation between two sets of data, the correlation coefficient will be close to zero. 23In Table 1, we show the calculated correlation coefficients.Since the coefficients are far from zero, there should be some correlations between λ, σ 2 , and ⟨r⟩.We can see that the quantity σ 2 is as good as a possible indicator of quantum chaos as the ratio ⟨r⟩.

Krylov state complexity in the stadium billiard
In this section, we summarize the numerical results on the Krylov state complexity for the stadium billiard.In Sec. 4 we found that the Lyapunov exponent, variation of the Lanczos coefficients, and the ratio ⟨r⟩ are correlated with each other for the Krylov operator complexity.We will observe similar tendency for the state complexity.The time dependence of the state complexity, however, shows behavior qualitatively different from that of the operator complexity.Figure 8 shows the Lanczos coefficient a n , b n for the state complexity on the stadium billiard in the integrable case (a/R = 0) and chaotic case (a/R = 1).We find tendency similar to that of the operator complexity, that is, the variation of the Lanczos coefficients becomes small in a chaotic system.We also find that not only b n but also a n , which appears only for the state complexity, has a smaller variation when the system is chaotic.Later in Sec.5.2, we will quantitatively show the correlation of the variation of the Lanczos coefficients with various indicators of the chaos.

The time dependence of Krylov state complexity
Figure 9(a) shows the time dependence of the state complexity for the stadium billiard with 0 ≤ a/R ≤ 1.Despite the system is chaotic for a/R ̸ = 0, the state complexity grows almost linearly in time at early time, neither exponentially nor polynomially.
The growth pattern of the complexity at early time, rather than the growth rate however, shows a clear correlation with the chaos.The complexity reaches the maximum  value at t ∼ 0.2 before it settles down to the asymptotic value.Although the asymptotic value is insensitive to the shape (a/R) of the stadium billiard (Fig. 9(b)), the peak value C max of the complexity depends on a/R rather smoothly (Fig. 9(c)).Hence, for the state complexity, the variation of the Lanczos coefficients is reflected not in the asymptotic value but in the peak value at early time.

Correlation between Lanczos coefficients and chaos
Now we turn to the correlation between the variation of the Lanczos coefficients and other indicators of the chaos, namely, the Lyapunov exponent λ and the ratio ⟨r⟩.Prior to such an analysis, let us show the dependence of the variance σ 2 a,b of the Lanczos coefficients on the parameter a/R in Fig. 10.This figure shows that the variances (both σ 2 a and σ 2 b ) decrease as a/R increases in the range 0 < a/R ≲ 0.3, then it stays at the asymptotic value for a/R ≳ 0.3.Since the Lyapunov exponent is monotonically increasing with respect to a/R, it is expected that σ 2 a,b and λ are negatively correlated.Such an expectation can be confirmed by explicitly plotting the relationship between σ 2 a,b and λ. Figure 11(a) shows the relationship between λ and σ 2 a,b , in which the negative correlation between these quantities can be observed.Since λ and ⟨r⟩ are positively correlated as explained in the previous section, ⟨r⟩ and σ 2 a,b are negatively correlated as shown in Fig. 11(b).Quantitative values of the correlations are summarized in Table 2, which clearly shows the negative correlations of σ 2 a,b with λ and ⟨r⟩.
6 Universality: the case of the Sinai billiard In this section, we analyze Krylov complexity of the Sinai billiard, which is another typical chaotic system, to show that our results in the previous section is universal.In Fig. 12, we show the shape of the Sinai billiard which we consider.The system is chaotic for l > 0 while the system becomes integrable for l = 0 [47].In Fig. 13, we show the l/L dependence of the Lyapunov exponent λ and the ratio ⟨r⟩.Although the Lyapunov exponent increases gradually with l/L, the ratio increases rapidly around l/L = 0.1.

Krylov operator complexity
Using the method described in Sec. 3, we numerically compute Krylov operator complexity for the truncated momentum operator P , with truncation N max = 100. 24In Fig. 14, we show the Lanczos coefficients for l/L = 0.05 and l/L = 1. 25 We identify the dimension of the Krylov space K P as K P = 9900.Obviously, the Lanczos coefficients for l/L = 0.05 distribute much broader compared to l/L = 1.On the other hand, their initial behaviors are almost identical.
In Fig. 15, we show the Krylov operator complexities as functions of t for the stadium billiards with l/L = 0.05, 0.1, 0.2, • • • , 1.While the Sinai billiard is chaotic for l/L > 0, the early time growth of the Krylov operator complexity is not exponential.The result for l/L = 0.05 is different from the others.This is because the behavior of the Lanczos coefficients changes abruptly below l/L = 0.1.In Fig. 16, we show the variance (2.11) as a function of l/L.The variance becomes larger in the integrable regime compared to the chaotic regime.Similarly to the ratio given in Fig. 13(b), the variance changes rapidly with l/L.
In Fig. 17, we show scatter plots, where the points are sampled from 0.01 ≤ l/L ≤ 0.2.26In Table 3, we show the correlation coefficients (4.1) calculated from these plots.Since the correlation coefficients are far from zero as easily expected directly from the plots, there should be some correlations between λ, σ 2 , and ⟨r⟩.Again, we see that the quantity σ 2 works as an indicator of quantum chaos as good as the ratio ⟨r⟩ does.

Krylov state complexity
We briefly summarize the results for the Krylov state complexity of the Sinai billiard in this section.The results are in parallel with those for the stadium billiard, which show the universality of the properties of the chaos indicators studied in this work.The numerical λ vs σ 2 -0.899970 ⟨r⟩ vs σ 2 -0.872723 λ vs ⟨r⟩ 0.924828 Table 3.The correlation coefficients between λ, ⟨r⟩, σ 2 of the Krylov operator complexity for the Sinai billiard.setup is the same as that for the stadium billiard analyzed in Sec. 5 except that the billiard table is switched to that given in Fig. 12.
The Lanczos coefficients a n , b n for the Krylov state complexity (Fig. 18), the time dependence of the Krylov state complexity (Fig. 19), and the variance σ 2 a,b of the Lanczos coefficients (20) for the Sinai billiard are qualitatively similar to those for the stadium billiard given in Sec. 5.One of the differences is that, in Fig. 19(a), the complexity C(t) becomes oscillatory when the billiard is integrable (l/L = 0), while such a periodicity was not observed in the stadium billiard (Fig. 9(a)).This feature is attributed to the fact that the energy spectrum of the Sinai billiard with l/L = 0 is commensurable, which is not the case for the stadium billiard with a/R = 0.
As shown in Fig. 21 and Table 4, the variances σ 2 a,b are correlated with the indicators of classical and quantum chaos, namely λ and ⟨r⟩.Hence the variance of Lanczos coefficients is a faithful indicator of chaos also for the Sinai billiard.
On top of it, the peak value C max of the complexity at early time in the time evolution (Fig. 19(c)) depends rather smoothly on the parameter l/L compared to the average of the complexity at late time ⟨C late time ⟩ (Fig. 19(b)).In this sense, C max is more faithful compared to ⟨C late time ⟩ as an indicator of quantum chaos.This tendency is in common with the stadium billiard, hence it is suggested that this feature is universal.

Summary and discussions
In this paper, we studied the Krylov complexity in quantum mechanics.To investigate the relationship between complexity and chaos, we numerically studied Krylov complexity in billiard systems.The stadium billiard, which is chaotic, allows a one-parameter deformation of its shape, and in a limit it reduces to the circular billiard which is integrable.We observed that under this deformation there exists a significant correlation between the classical Lyapunov exponents and the variances of the Lanczos coefficients for both the Krylov operator complexity and the Krylov state complexity.We also found a significant correlation between the variances of the Lanczos coefficients and the statistical distribution of the adjacent spacings of the quantum energy levels.This suggests that the variances of the Lanczos coefficients are a good indicator of quantum chaos as well as the energy level statistics.Similar results were confirmed for the one-parameter family of the Sinai billiard, which suggests that in more general quantum mechanical systems the variances of the Lanczos coefficients can be a good indicator of quantum chaos.Note that in our numerical analysis a finite cutoff N max in the quantum energy spectrum was necessary.The descent part of the sequence of the Lanczos coefficients, which we used in the evaluation of the variance, would be lost in the limit where N max goes to infinity since the Lanczos algorithm would not always terminate in that limit.In this point of view, the variance of the Lanczos coefficients measures the response of the system to a regularization N max , and this response distinguishes the chaoticity.
In viewing the structure of the resultant correlations, we observe one issue: when the billiard is deformed, there is an abrupt change in variances of the Lanczos coefficients and the energy level statistics, while the change in classical Lyapunov exponents is mild.This suggests that there may be some discrepancy between the notion of classical chaos and the currently proposed quantum chaos.We leave to a future work the quest for the cause of this discrepancy and a search for a better indicator of quantum chaos that more accurately reflects classical chaoticity, with a consistent understanding of how classical chaoticity manifests itself through the classical limit from quantum theory.
Here let us discuss the late-time behavior of the complexity in more detail.As briefly mentioned in Sec. 2, there was an expectation that the late-time value of the Krylov operator complexity would be larger in chaotic systems compared to integrable systems [16].This behavior may be understood [16] in terms of the Anderson localization on the one-dimensional chain (2.6). 27Intuitively, when the Lanczos coefficients behave erratically, it becomes difficult for the amplitude to spread smoothly along the chain and the amplitude will be localized around the initial position.Since the Krylov complexity represents the expectation value of the position along the chain, this means that the complexity stays at a small value.By the same argument, when variances of the Lanczos coefficients are small, the Krylov complexity is expected to take a larger value.
In the stadium billiard, as shown in Fig. 5, the saturation value of the complexity becomes smaller as the system approaches the integrable regime (a/R → 0).However, this late-time behavior is not so obvious in the Sinai billiard.In Fig. 15, the complexity for l/L = 0.05 appears to be small, but it does not show any saturation in the displayed range.It also seems that the complexity for l/L ≥ 0.05 approaches the same value.In App.A, we also consider the Krylov operator complexity for l/L = 0, whose time dependence is shown in Fig. 25(c).It is possible that the complexity remains small and oscillates for l/L = 0, but this result should not be seriously trusted because of the numerical instabilities.Therefore, our tentative conclusion is that it may be difficult to distinguish the chaoticity from the late-time behavior of the Krylov operator complexity.This also implies that the effect of the erratic behavior of the Lanczos coefficients on the one-dimensional chain (2.6) is more subtle.The authors of [49] have reported similar results in Ising spin chains that the late-time behavior of the complexity strongly depends on the choice of the operator and is not correlated with chaoticity although there was a correlation between the variance of the Lanczos coefficients and chaoticity.
Even for the Krylov state complexity, its asymptotic value is insensitive to the billiard table shape, hence its relationship with the variance of the Lanczos coefficients is subtle as it was for the Krylov operator complexity.One difference is that the peak value of the Krylov state complexity at early time shows some correlation with the classical Lyapunov exponent, which was discussed in [12,13].It would be interesting to study this property in more detail and to examine how universal this behavior is for various systems showing quantum chaos.
As concluding words of this paper, let us make a remark from a broad perspective.The numerical method in this paper is valid for general quantum mechanical systems.As in the case of OTOCs, it is important to investigate the universal nature of chaos by performing similar analyses specifically on various quantum mechanical systems such as, for example, polygonal billiards [50], coupled harmonic oscillators [51][52][53][54][55] 28 , inverted harmonic oscillators [9,10], and anharmonic oscillators [56].Evaluating the variances of the Lanczos coefficients in those popular examples, one can compare them with classical/quantum Lyapunov exponents, leading to a broad and universal view of quantum chaos and complexity.
Readers are reminded of the history that the chaos and the complexity have been intensively studied in the context of the AdS/CFT correspondence and black holes.In that viewpoint, clarifying the relationship between conventional computational complexity and Krylov complexity is another important issue.There were proposals for the gravity dual interpretation of the conventional complexity, and recently, a dual interpretation of Krylov complexity has been studied [57].A further exploration of the Krylov complexity in various quantum mechanics may reveal the mystery of the interior of black holes.billiard system is integrable for a/R = 0. Indeed, the energy eigenvalues and the energy eigenfunctions can be obtained analytically: where J k is the Bessel function of the first kind and ρ kl is the l-th zero of J k .The matrix element of the momentum operator can be obtained by numerical integration.

A.2 Sinai billiard
The shape of the Sinai billiard is determined by l/L (see Fig. 12).In Fig. 24, we show the resolution dependence of the numerical results for l/L = 0.5.Although σ 2 slightly drops at resolution = 0.0002, resolution = 0.0001 is sufficient for convergence as a whole.In Fig. 25, we show the resolution dependence of the numerical results for l/L = 0, at which the system becomes integrable.Although the Frobenius norm of P is well convergent, the variance σ 2 does not converge and increases as the resolution decreases.Correspondingly, the Krylov operator complexity does not also converge.This behavior is different from that of the stadium billiard.Notice that although the numerical calculation is unstable for l/L = 0, the variance increases and the Krylov operator complexity decreases, respectively, when the resolution becomes smaller.This is consistent with the general expectation   that for integrable systems, the variance becomes large and the late-time asymptotic value of the Krylov operator complexity becomes small.To explore the possible causes of the irregular dependence on the resolution, we consider the analytic solution of the Schrödinger equation.We can regard the analytic solution as resolution = 0.The energy eigenvalues and the energy eigenfunctions for l/L = 0 are Then, the matrix element of the momentum operator can be obtained by analytic integration.Due to the orthogonality of the trigonometric functions, there may be many zero elements in the analytically obtained P .In numerical calculations, these elements can have small non-zero values due to numerical errors.These errors do not significantly affect the Frobenius norm of P .However, when the numerically obtained P is input to the Lanczos algorithm, small numerical errors on elements that should be zero can have a significant impact on the final result due to the nonlinear nature of the algorithm.

B Truncation dependence
In our numerical calculation of the Krylov complexity, we considered only a finite number of energy levels by introducing a truncation N max to the spectrum.In this section, we will λ vs σ 2 -0.497215 ⟨r⟩ vs σ 2 -0.214267

B.2 Krylov state complexity
We summarize the numerical results for the Krylov state complexity with reduced truncation order N max = 250 and compare them with those shown in Sec. 5 for N max = 500.
Figure 28 shows the time dependence and peak value of the Krylov state complexity and also the variance of the Lanczos coefficients for N max = 250.Qualitative features are the same as those in Figs. 9 and 10, while the asymptotic value of the complexity is reduced to ∼ 125 = N max /2.Other difference is that the variance σ 2 a , σ 2 b of the Lanczos coefficients and also the fluctuation of the complexity slightly increase compared to the case with N max = 500.
Correlation between the indicators of classical and quantum chaos are summarized in Fig. 29 and Table 6.The correlation coefficients in Table 6 are as large as those for N max = 500 in Table 2, which the variance of the Lanczos coefficients σ 2 a , σ 2 b can be used as an indicator of chaos as long as N max is moderately large.

C Initial operator/state dependence of the Krylov complexity
In this section, we consider the dependence of the Krylov complexity on the choice of the initial operator/state.

C.1 Krylov operator complexity
Since we studied Krylov state complexity for the flat initial state, we consider the flat operator, which is considered in appendix B of [17], with N max = 100 in this section.In Fig. 30, we show the Krylov operator complexity and the variance of the Lanczos coefficient.For the flat operator, the Krylov operator complexity, especially its saturation value, hardly depends on a/R [17].Also, the variance of the Lanczos coefficient cannot distinguish chaoticity and integrability in this case.In Fig. 31, we cannot find any correlation.The correlation coefficients take positive values as shown in Table 7 even though they should be negative.Thus, the Krylov operator complexity and the behavior of the Lanczos coefficients depend on the choice of the initial operator.In order to distinguish chaoticity, we need to choose the initial operator properly.Although we do not know the criterion for this choice, we expect that the fundamental operators, which appear in the definition of the Hamiltonian, should be proper.

C.2 Krylov state complexity
The results on the Krylov state complexity are based on the flat initial state (3.10).To examine the dependence of the results on the choice of the initial state, in this appendix   λ vs σ 2 0.118923 ⟨r⟩ vs σ 2 0.204392 For numerical calculation below, we chose (x 0 , y 0 ) = (0.45, 0.37), σ = 0.075 and p 0 = 40 so that the distribution of the wave packet is contained well within the billiard table for any a/R.We show the profile of the wave packet and the corresponding components of the initial state Ψ = (Ψ 1 , . . ., Ψ Nmax ) T for a/R = 1 in Fig. 32.We use Eq.(C.2) for any a/R, and the energy eigenstates and the components Ψ n=1,...,Nmax with respect to them are calculated for each a/R.For any a/R, the components Ψ n take nonzero value fluctuating  around zero for n ≲ 300, while they tend to zero for n ≳ 300.
The numerical results for N max = 500 are shown in Figures 33,34 and 35.They are qualitatively the same as those for flat initial state in Sec. 5 except for the following features: (i) the Lanczos coefficients show linear growth for the first few indices (a n , b n for n ≲ 4, see Fig. 33); (ii) fluctuation of the Krylov state complexity around the late time saturation value is relatively large, and the peak structure at early time evolution is less obvious compared to that for the flat initial state (see Fig. 34(a)).
In Table 8, we show the correlation coefficients for the indicators of the chaos for the wave packet initial state with N max = 500.Although these values are slightly smaller compared to those for the flat initial state, they are still far from zero and indicate that the correlations between λ, σ 2 a,b and ⟨r⟩ are strong.For comparison, we also show the correlation coefficients for N max = 250.From these values we can observe some correlations, though they are slightly weaker than those N max = 500.This result support that the variances of the Lanczos coefficients are robust probe of the chaos once the trancation order N max is moderately large.

2. 1 . 1
Krylov operator complexityConsider a quantum system with Hamiltonian H.The Krylov operator complexity of a given operator O is defined as follows.To begin with, we can rewrite the Heisenberg operator O(t) = e iHt O(0)e −iHt as . If b n = 0 stop; otherwise set O n = A n /b n and go to step 3.If the dimension K O of the Krylov space is finite, the above algorithm ends with b K O = 0.The algorithm produces the orthonormal basis {O n } K O −1 n=0 called the Krylov basis, and a set of positive numbers {b n } called the Lanczos coefficients.Now expand the Heisenberg operator O(t) in terms of the Krylov basis as

R a 1 Figure 1 .
Figure 1.Geometry of the stadium billiard.Dirichlet boundary conditions are imposed on the boundaries.

Figure 2 .
Figure 2. The a/R dependence of the Lyapunov exponent and the ratio.

Figure 3 .
Figure3.The Lanczos coefficients for the truncated momentum operator P in stadium billiards with a/R = 0 (blue dots) and a/R = 1 (orange dots).Note that the horizontal axis is in log scale.The inset is the enlarged version, where the data are used to calculate the variance.

Figure 4 .Figure 5 . 2 Figure 6 .
Figure 4.The time dependence of Krylov operator complexity for various values of a/R.

Figure 7 .
Figure 7.The scatter plots of (a) the Lyapunov exponents and the variances, (b) the ratios and the variances, (c) the Lyapunov exponents and the ratios.The points are sampled from 0 ≤ a/R ≤ 0.5.

Figure 8 .
Figure 8.The Lanczos coefficients of the state complexity for equally-distributed initial state in stadium billiards with a/R = 0 (blue dots) and a/R = 1 (orange dots).The inset is the enlarged version, where the data are used to calculate the variance.
The time dependence of Krylov state complexity for various values of a/R.The a/R dependence of the late-time value of Krylov state complexity.The a/R dependence of the peak value of Krylov state complexity.

Figure 9 . 2 Figure 10 .
Figure 9.The a/R dependence of Krylov state complexity.Panel (a): time dependence of the complexity.Panels (b), (c): the late-time average and the peak value of the complexity.The latetime average of the complexity is taken over the time range 1 < t < 20.

Figure 12 .
Figure 12.Geometry of the Sinai billiard.Dirichlet boundary conditions are imposed on the boundaries.

Figure 13 .
Figure 13.The l/L dependence of the Lyapunov exponent and the ratio.

Figure 14 .Figure 15 .
Figure 14.The Lanczos coefficients for the truncated momentum operator P in Sinai billiards with l/L = 0.05 (blue dots) and l/L = 1 (orange dots).Note that the horizontal axis is in log scale.The inset is the enlarged version, where the data are used to calculate the variance.

2 Figure 16 .Figure 17 .
Figure 16.The variance σ 2 as a function of l/L of the Krylov operator complexity for the Sinai billiard.

Figure 18 .
Figure 18.The Lanczos coefficients of the Krylov state complexity of the Sinai billiard with l/L = 0 (blue dots) and l/L = 1 (orange dots).The inset is the enlarged version, where the data are used to calculate the variance.
The time dependence of Krylov operator complexity for various values of l/L.〈C late time 〉 (b) The l/L dependence of the late-time value of Krylov state complexity.The l/L dependence of the peak value of Krylov state complexity.

Figure 19 . 2 Figure 20 .
Figure 19.The l/L dependence of Krylov state complexity for the Sinai billiard.Panel (a): time dependence of the Krylov state complexity.Panels (b), (c): the late-time average and the peak value of the Krylov state complexity.The late-time average of the complexity is taken over 1 < t < 20.

Figure 21 .λ vs σ 2 a
Figure 21.The scatter plots of (a) the Lyapunov exponents and the variances and (b) the ratios and the variances.λ vs σ 2 a -0.741803 λ vs σ 2 b -0.757869 ⟨r⟩ vs σ 2 a -0.965785 ⟨r⟩ vs σ 2 b -0.962833 The Krylov operator complexities as functions of time for resolution = 0.0001 and 0.0002.

Figure 22 . 2 (
Figure 22.The resolution dependence of the Frobenius norm of P , the variance σ 2 and the Krylov operator complexity for a/R = 0.5.

Figure 23 .
Figure 23.The resolution dependence of the Frobenius norm of P , the variance σ 2 and the Krylov operator complexity for a/R = 0.
The Krylov operator complexities as functions of time for resolution = 0.0001 and 0.0002.

Figure 24 . 2 (
Figure24.The resolution dependence of the Frobenius norm of P , the variance σ 2 and the Krylov operator complexity for l/L = 0.5.

Figure 25 .
Figure 25.The resolution dependence of the Frobenius norm of P , the variance σ 2 and the Krylov operator complexity for l/L = 0.

2 (
b) Variances σ 2 a,b as functions of a/R.The peak value of Krylov state complexity.

Figure 28 .
Figure 28.The a/R dependence of Krylov state complexity.Panel (a): time dependence of the complexity.Panel (b): variance σ 2 a , σ 2 b of the Lanczos coefficients.(c): the peak value of the complexity.

2 (
The time dependence of Krylov operator complexity for various values of a/R.b) The variance σ 2 as a function of a/R.

Figure 30 . 2 (b) ⟨r⟩ vs σ 2
Figure 30.The a/R dependence of Krylov operator complexity and the variance of the Lanczos coefficients.

Figure 31 .
Figure 31.The scatter plots of (a) the Lyapunov exponents and the variances and (b) the ratios and the variances.

Figure 34 .Figure 35 .
Figure 34.The a/R dependence of Krylov state complexity for the wave packet initial state (C.2).Panel (a): time dependence of the complexity.Panel (b): variance σ 2 a , σ 2 b of the Lanczos coefficients.

(a) Nmax = 500 λ vs σ 2 a
. Note that the Schrödinger state |ψ(t)⟩ = e −iHt |ψ(0)⟩ is a linear combination of and go to step 3.If K |ψ⟩ is finite, then this Lanczos algorithm ends with b K |ψ⟩ = 0.The resulting orthonormal basis {|K n ⟩} K |ψ⟩ −1 n=0 is called the Krylov basis.Notice that there are two sets of Lanczos coefficients {a n } and {b n } in this case.Expanding the Schrödinger state |ψ(t)⟩ in terms of the Krylov basis as 18Then, we numerically solve(2.6)and obtain the Krylov operator complexity (2.7).Also, the variance (2.11) can be obtained from the Lanczos coefficients.In our calculation, the variance (2.11) is evaluated for the region N 2 max /200 ≤ i ≤ N 2 max /40.Let us describe how to calculate the Krylov complexity for a state |ψ⟩.Similarly to the calculation of Krylov operator complexity, Krylov state complexity can be calculated by expanding |ψ⟩ in terms of the energy eigenstates |n⟩ as a basis.Solving the Schrödinger equation numerically, we can prepare the state vector as

Table 2 .
Correlations between λ, ⟨r⟩, σ 2 a,b for the state complexity of stadium billiard.

Table 4 .
Correlations between λ, ⟨r⟩, and σ 2 a,b for the Krylov state complexity of Sinai billiard.

Table 6 .
Correlations between λ, ⟨r⟩, σ 2 a,b for the state complexity of stadium billiard.

Table 7 .
The correlation coefficients between λ, ⟨r⟩, σ 2 for the Krylov operator complexity.weshow the results for the initial state corresponding to a wave packet on the Stadium billiard table.We use a wave packet defined byΨ(x, y) = A exp − (x − x 0 ) 2 + (y − y 0 ) 2