Ultimate speed limits to the growth of operator complexity

In an isolated system, the time evolution of a given observable in the Heisenberg picture can be efficiently represented in Krylov space. In this representation, an initial operator becomes increasingly complex as time goes by, a feature that can be quantified by the Krylov complexity. We introduce a fundamental and universal limit to the growth of the Krylov complexity by formulating a Robertson uncertainty relation, involving the Krylov complexity operator and the Liouvillian, as generator of time evolution. We further show the conditions for this bound to be saturated and illustrate its validity in paradigmatic models of quantum chaos. In quantum isolated system, operator growth can be quantified by Krylov complexity. Here, the authors establish a rigorous bound on the Krylov complexity growth rate based on the uncertainty principle and show that the presence of quantum chaos is not strictly necessary to saturation of the bound.

Q uantum speed limits (QSL) impose fundamental constraints on the pace at which a physical process can unfold. Since their conception 1,2 , they have been formulated as bounds on the minimal time at which a distance between quantum states can be traversed. The freedom in the choice of the distance can be used to sharpen the discrimination between quantum states, and with it, the notion of the speed of evolution 3,4 . Additional efforts have been devoted to exploring the role of the underlying dynamics, generalizing early results from isolated systems to open [5][6][7][8] and classical processes 9,10 . The resulting speed limits have become a useful tool in various branches of physics, ranging from information processing 11 to many-body physics 12 , quantum control 13 , and quantum metrology 14 . However, traditional QSL are too conservative in estimating the relevant time scales in many processes, such as thermalization 15 . This has motivated the development of speed limits suited for specific measures and observables 16 , as in the pioneering work by Mandelstam and Tamm 1 . In this sense, certain speed limits follow from generalized uncertainty relations such as those derived by Heisenberg and Robertson 17 .
In parallel with the study of QSL, quantifying the complexity of a physical process is a central task for the advancement of fundamental physics and quantum technologies. Lloyd pointed out that the computational complexity of physical processes is limited by QSL 18 . Analogously, the circuit complexity of a quantum state 19 , defined as the number of elementary operations required to generate it from a reference state, can be characterized in terms of conventional QSL [20][21][22][23] . A complementary approach for manybody quantum systems focuses on the buildup of complexity in the time evolution of an initial local observable, known as operator growth [24][25][26][27][28] . The intuition is that simple operators unitarily evolve into increasingly complex ones. Quantum information initially encoded in a few degrees of freedom is thus scrambled over the system in the course of evolution, making it impossible to recover it through local measurements and giving rise to thermalization. The unambiguous description of this scrambling process remains an open problem. One possibility is to probe it via an out-of-time-ordered correlator 29,30 that may be used to identify an analog of the Lyapunov exponent, providing a connection with classical chaos, e.g., the butterfly effect. Such quantum Lyapunov exponent obeys a universal upper bound 30 , which helps refine the notion of maximal chaos, is saturated by black holes and is further tied to the eigenstate thermalization hypothesis 31,32 . A related approach, which we shall pursue in this work, is to study the dynamical evolution of operators in Krylov space, exploited in numerical techniques such as the recursion method 33 . In this context, operator growth is quantified by the so-called Krylov complexity, a measure of the delocalization of the time-dependent operator in the Krylov basis [34][35][36][37][38] . The authors of 34 made a conjecture on the universal operator growth, namely, that Krylov complexity can grow at most exponentially, and it does so in generic non-integrable systems. Remarkably, its growth rate upper bounds the Lyapunov exponent, establishing a connection with the bound on out-of-time-ordered correlators 30,39 . Further studies have shown that exponential operator growth is possible in free and integrable systems 40 , while the role of the interaction graph in a quantum network has been explored in 41 .
Here, we characterize the growth of Krylov's complexity by deriving a fundamental limit on its rate of change and by studying analytically the conditions under which this bound is saturated. Our results show that saturation, which is also found to correspond to a particular notion of minimum uncertainty, occurs whenever the dynamical evolution of the system has the underlying structure of a three-dimensional complexity algebra, which was introduced by 42 . In this setting, the unitary evolution of an operator can be represented as the displacement of generalized coherent states 42 , which display classical-like behavior 43 . As demonstrated in several paradigmatic examples, the saturation of the growth rate may be possible in some chaotic systems, but quantum chaos is not required for it.

Results and discussion
Quantum dynamics in Krylov space. Consider an isolated quantum system in which the time evolution of an observable O is generated by a time-independent Hamiltonian H according to the Heisenberg equation of motion ∂ t OðtÞ ¼ i½H; OðtÞ, setting ℏ = 1. The solution to this equation with the initial condition Oð0Þ ¼ O is given by OðtÞ ¼ e itH Oe ÀitH . In terms of the Liouvillian superoperator given by L ¼ ½H; Á, the Taylor expansion of the time-evolving observable OðtÞ ¼ ∑ 1 n¼0 ðitÞ n n! L n O shows that its dynamics is contained in the complex linear span of the operators fL n Og 1 n¼0 . This span is completely determined by the Hamiltonian and the initial observable and is known as the Krylov space.
From now on, we consider the restriction of each operator and superoperator to the Krylov space. To highlight the vector space structure, we make use of the bracket notation A j Þ when expressing operator A in an equation. We choose to equip the Krylov space with an inner product satisfying the properties An example of a family of inner products satisfying these two properties is given by ðAjBÞ ¼ e βH=2 A y e ÀβH=2 B β . The bracket Á h i β denotes the thermal expectation value with respect to the equilibrium Gibbs state e −βH /Z and thus (A|B) reduces to the Hilbert-Schmidt inner product when β = 0, up to a normalization factor. It follows from the second property of the inner product that the operators O and LO are orthogonal. Let b 0 ¼kOk and b 1 ¼kLOk, where ∥⋅∥ is the norm induced by the inner product. By starting from the normalized vectors O 0 ¼ O=b 0 and O 1 ¼ LO=b 1 , we can construct an orthonormal basis fO n g DÀ1 n¼0 for the Krylov space by applying the Lanczos algorithm. This algorithm works as follows: given the first n + 1 basis vectors, one constructs the orthogonal vector jA nþ1 Þ ¼ LjO n Þ À b n jO nÀ1 Þ, where b n = ∥A n ∥ and then normalize it to obtain jO nþ1 Þ. We call the constructed basis the Krylov basis. It is possible that the Krylov dimension D is infinite, in which case the Lanczos algorithm never halts. We remark that the Lanczos algorithm is only guaranteed to construct an orthonormal basis if the Liouvillian is self-adjoint, i.e., the first property of the inner product is satisfied. Generally, the Lanczos algorithm involves a third term on the right-hand side of the equation for jA nþ1 Þ. This term is, however, always zero whenever the second property of the inner product is satisfied. Thus, with our chosen inner product, the action of the Liouvillian on the Krylov basis takes a specific form LjO n Þ ¼ b nþ1 jO nþ1 Þ þ b n jO nÀ1 Þ. As pointed out in ref. 42 , this motivates one to consider abstract raising and lowering operators that we denote by L þ and L À , respectively. Their action on the Krylov basis is given by L þ jO n Þ ¼ b nþ1 jO nþ1 Þ and L À jO n Þ ¼ b n jO nÀ1 Þ. The Liouvillian can then be expressed as their sum.
It is further convenient to introduce the real-valued functions φ n (t), which appear in the expansion of OðtÞ as jOðtÞÞ ¼ 1 kOk ∑ DÀ1 n¼0 i n φ n ðtÞjO n Þ. We will refer to these functions as the amplitudes of the observable. These amplitudes evolve according to the recursion relation ∂ t φ n (t) = b n−1 φ n−1 (t) − b n φ n +1 (t) with the initial conditions φ 0 (0) = 1 and φ n (0) = 0 for n > 0. Thinking of the Krylov basis vectors as forming the sites of a onedimensional lattice, b n can be interpreted as a hopping amplitude, see, e.g., 34,35 . In this sense, one can think of O as a onedimensional discrete wave function that is initially localized and then spreads out over the lattice as time evolves. An increase in the population of the sites further away from the origin reflects a greater increase of complexity of the observable. In order to quantify this, it is natural to consider the Krylov complexity of OðtÞ, defined to be The main task of our work is to bound the growth of Krylov's complexity. Due to unitary dynamics, the norm of the evolution is preserved and the Krylov complexity is unchanged if one normalizes the operators studied. We will, therefore, without loss of generality, consider O to be normalized. By introducing the complexity operator K ¼ ∑ DÀ1 n¼0 njO n ÞðO n j, which plays the role of the position operator in the Krylov lattice, it is possible to express Krylov complexity as the "expectation value" of K with respect to OðtÞ. More precisely, if K h i t ðOðtÞjKOðtÞÞ then Dispersion bound on Krylov complexity. If the Krylov space forms an inner product space in which A and B are self-adjoint superoperators, then there ought to exist a Robertson uncertainty relation given by ΔAΔB is the dispersion of A with respect to some state A j Þ. When the Krylov dimension is infinite, it is necessary that A j Þ is contained in the intersection between the domains of AB and BA, otherwise the inequality might not hold 44 In other words, the growth of Krylov complexity is upper bounded by a constant times the dispersion of the complexity operator. By defining a characteristic time scale τ K ¼ ΔK= ∂ t KðtÞ , one obtains τ K b 1 ≥1/2, which takes the form of a Mandelstam-Tamm bound, and emphasizes the role of b 1 ¼k LO k as a norm of the generator of evolution in Krylov space. To avoid confusion with the uncertainty relation for observables, we will refer to this bound as the dispersion bound. We note that no bound tighter than (2) can be found by considering the more general Schrödinger uncertainty relation, as the extra term given by the anticommutator identically vanishes, as shown in Methods. It is not self-evident that saturation of the dispersion bound can be achieved under the unitary dynamics of the observable. There are very specific relations between L, O and K that need to hold: the Liouvillian is required to be tridiagonal in the eigenbasis of the complexity operator and the initial state of the observable is required to be parallel to the eigenvector with the lowest eigenvalue. The conditions for the saturation of the dispersion bound are thus highly constrained and differ from those known for saturation of a Robertson uncertainty relation in general. The required conditions admit a geometrical interpretation, elaborated in Methods. The bound is saturated if and only if the evolution curve moves along the gradient of the Krylov complexity. This requires that the dynamics is directed along the direction that maximizes the local growth of complexity; see Methods. The only exception involves extremal points in which any direction away from the extremal point leads to saturation. This is indeed the case for t = 0. Indeed, there exists Liouvillians of the form L ¼ L þ þ L À for which the tangent of the generated path will be parallel with the gradient for all times.
Saturation of the dispersion bound. Time evolutions saturating the dispersion bound are characterized by a unique algebraic structure. Define the superoperator B ¼ L þ À L À . Following 42 , we consider their simplicity hypothesis: namely, the assumption that L, B and the commutatorK ¼ ½L; B close an algebra with respect to the Lie bracket. It was shown in 42 that this forcesK to be related to the complexity operator viaK ¼ αK þ γ, where α; γ 2 R. We show in Supplementary Note 2 that γ is a positive number and α is a real number satisfying the condition α ≥ 0 for infinite Krylov dimension and α ¼ À 2γ DÀ1 for finite Krylov dimension. Moreover, the only possible closure of the algebra is given by the commutation relations Given this algebra, the evolving observable can be interpreted as a curve of generalized coherent states evolving according to the displacement operator DðξÞ ¼ e ξL þ ÀξL À , where ξ = it. Moreover, the initial state is the highest weight state of the representation, which is annihilated by L À by construction. Coherent states can be viewed as the states closest to the classical ones in the sense that they typically minimize an uncertainty relation. It is for example known that coherent states of the Harmonic oscillator saturate the Robertson uncertainty relation for the pair of observables of position and momentum. Building on this intuition, we could expect that the dispersion bound is saturated for the simplicity hypothesis. It turns out that this intuition is indeed correct. In fact, as we show in Supplementary Note 2, the dispersion bound is saturated if and only if the simplicity hypothesis holds. The saturation of the dispersion bound dictates the evolution of the Krylov complexity, where three different scenarios are possible, as shown in Fig. 1a. The growth of complexity at the speed limit is described by the differential equation with the conditions that K(0) = 0 and K(− t) = K(t). For finite Krylov dimension, saturation of the dispersion bound sets the complexity growing according to . In this, case, the corresponding complexity algebra (3) reduces to the SU(2) algebra. By contrast, for infinite Krylov dimension there are two distinct scenarios for the complexity growth: for α > 0 one finds KðtÞ ¼ 2γ α sinh 2 ffiffi α p t 2 , while for α = 0 the solution reads KðtÞ ¼ γ 2 t 2 . The complexity algebra in these two cases reduces to SLð2; RÞ and the Heisenberg-Weyl algebra (HW), respectively. Reference examples maximizing the Krylov complexity growth rate at all times are discussed in Supplementary Note 1. One such example with α > 1 is the Sachdev-Ye-Kitaev (SYK) model 45 , a paradigm of quantum chaos. However, the saturation of the bound does not require quantum chaos and can indeed be achieved by a single qubit, with α = 0 (Supplementary Note 1). Together with the time-dependence of K(t) and the complexity algebra, the value of α also determines the growth of the Lanczos coefficients in the Krylov lattice. As proven in Supplementary Note 2, the dispersion bound is saturated if and only if the Lanczos coefficients grow according to exhibiting three different scalings as a function of α, see Fig. 1b. That the simplicity hypothesis implies (5) has already been pointed out in 42 . For α > 1 and large n, this dependence captures the linear growth b n ¼ ffiffiffi α p n conjectured by Parker et al. to hold in generic non-integrable systems, maximizing the Krylov complexity growth 34 .
Krylov complexity in generic systems. We next discuss the Krylov complexity growth in generic systems not fulfilling the simplicity hypothesis. We can use Eq. (5) to estimate when and at what time scale a generic system deviates from the bound. By expanding Krylov complexity up to fourth order, we find that Since we can always find a value on α and γ such that b 1 and b 2 satisfy (5), we conclude that the bound (2) is saturated up to the third order in time. By expanding the Krylov complexity up to sixth order, we find that the Lanczos coefficient b 3 will appear in the last term, and since we are not guaranteed to be able to find a value on α and γ such that b 1 , b 2 , and b 3 satisfy(5), we conclude that the system can only start deviating from the bound (2) as a result from fifth-order terms in the expansion. We can estimate this time scale by finding the value of t for which the third order coefficient of ∂ t K(t) is equal to its fifth-order coefficient. We will call this time the deviation time, denoted by τ d , and it is explicitly given by To get an understanding of the complexity growth in a generic setting, we next illustrate the Krylov dynamics of a system described by a random matrix Hamiltonian. Specifically, we consider the Krylov complexity of an ensemble EðHÞ of random matrix Hamiltonians, a paradigm of quantum chaos 46 . We sample the Hamiltonian matrices H from the Gaussian Orthogonal Ensemble GOE(d), where d is the dimension of the Hilbert space. We then calculate the Lanczos coefficients {b n } with partial reorthogonalization 36,47 . Specifically, we consider samples of real matrices H = (X + X ⊺ )/2, where all elements x 2 R of X are pseudo-randomly generated with probability measure given by the normal distribution, expðÀx 2 =ð2σ 2 ÞÞ=ðσ ffiffiffiffiffi 2π p Þ. In order to study the general behavior of Lanczos coefficients, we choose an initial observable, which is represented as the normalized vector jOÞ ¼ ð1=d; 1=d; ; 1=dÞ T , expressed in a fixed eigenbasis of the Liouvillian. However, the following results do not depend strongly on the choice of O, provided it is dense in the eigenbasis of the Hamiltonian. Figure 2a shows the squares of the Lanczos coefficients for a single realization and the average hfb n gi EðHÞ over 100 different Hamiltonians of dimension d = 32, sampled from GOE(d) with standard deviation σ = 1. Operator growth is displayed by the time-dependent amplitudes, which are found by solving the recursion relation and exhibit diffusion-like dynamics on the Krylov basis, shown for a single realization in Fig. 2b. The corresponding time evolution of Krylov complexity and its growth rate are shown in panels c and d, respectively. Hamiltonians sampled from GOE(d) behave as a generic system, given that the Lanczos coefficients do not, in general, grow according to (5), as shown in Fig. 2a. As a result, the growth rate starts deviating from the dispersion bound around the time scale τ d in Eq. (6), indicated by the vertical line in Fig. 2 c, d. In short, while GOE Hamiltonians provide a useful paradigm in the description of quantum chaotic systems, the dynamics generated by them do not maximize the growth of Krylov's complexity for t > τ d .
Our results establish the ultimate speed limit to operator growth in isolated quantum systems. Specifically, the dispersion bound governs the growth rate of Krylov complexity, playing the role of a Mandelstam-Tamm uncertainty relation in operator space. This bound is saturated by quantum systems in which the Liouvillian governing the time evolution fulfills a simple algebra. The latter arises naturally in certain quantum chaotic systems, such as the SYK model. However, other paradigmatic instances of quantum chaos, such as random matrix Hamiltonians, do not maximize the growth of Krylov complexity. Indeed, a saturation of the bound does not require quantum chaos and can be achieved, e.g., by a single qubit.

Methods
Vanishing of the anticommutator contribution in the Robertson uncertainty relation for K and L. We establish a universal feature of Krylov complexity, valid for any physical system: namely, that its anticommutator with the Liouvillian L has vanishing expectation value over the evolved operator jOðtÞÞ. The relevance of this result relies on the fact that this quantity enters the Schrödinger uncertainty principle for the two operators K and L where Let us now demonstrate that the anticomutator term in Eq. (9) is identically zero. By expanding jOðtÞÞ over the Krylov basis, we obtain which, by performing the sums over k and n, yields Since the amplitudes φ n and the coefficients b n are real quantities, comparing Eqs. (9) and (12) we immediately conclude that Let us note that the key condition to obtain this result is the fact that the Liouvillian connects only states that are nearest neighbors on the Krylov lattice so that we are left with a purely imaginary phase (−i) m (i) m±1 = ±i. It is this peculiar property that allows the Liouvillian to be interpreted as a sum of generalized ladder operators L ± 42 . However, let us point out that here we are not making any assumption regarding the commutation rules between these operators: we are considering the structure of Krylov space in full generality.
Moreover, from Eq. (8), we immediately obtain the relation between the anticommutator ½K; L and the complexity rate ∂ t K: Therefore, the Schrödinger uncertainty relation (7) can be recast as the dispersion bound (2) on the growth of Krylov complexity: Geometrical interpretation of the saturation of the bound. For the geometrical interpretation of the saturation of the bound, we assume the Krylov space to be of finite dimension. However, the results could potentially be extended to infinitedimensional Krylov spaces as well. The Krylov space is isomorphic to a 2D-dimensional real vector space, and we can therefore consider the Euclidean metric g, given by the real part of the inner product. The evolution curve of O will then be restricted to the unit sphere of the Krylov space. This unit sphere forms a Riemannian manifold and we can consider the Krylov complexity as a function on this manifold defined by KðAÞ ¼ ðAjKAÞ; for any element A j Þ in the Krylov space with a unit norm. In this sense, when we write K(t) we simply mean KðOðtÞÞ which is consistent with how we defined complexity for the evolution. The differential of Krylov complexity will be denoted by dK and its action on any tangent vector _ A at A is given by dKð _ AÞ ¼ ð _ AjAÞþ ðAj _ AÞ. This differential together with the metric can be used to define the gradient of Krylov complexity. It follows from the theory of differential geometry that the gradient of Krylov complexity at A, denoted by ∇KðAÞ, is the unique vector satisfying the expression gð∇KðAÞ; _ AÞ ¼ dKð _ AÞ for all tangent vectors _ A at A 48 . It can be checked that the gradient must then be given by ∇KðAÞ ¼ 2ðK À K h iÞA, which indeed is tangent to the unit sphere at A. The change of Krylov complexity along the curve OðtÞ, generated by the Liouvillian, is given by ∂ t KðtÞ ¼ gð∇KðtÞ; ∂ t OðtÞÞ, where ∇ K(t) is the gradient at OðtÞ. Applying the Cauchy-Schwarz inequality on the right-hand side gives us the inequality The right-hand side of this inequality is exactly 2b 1 ΔK and we note that it is saturated if and only if the tangent vector of OðtÞ is parallel to the gradient of Krylov complexity. We also note that the gradient is the zero vector at time zero and so the dispersion bound is always initially saturated. The unitary orbit of O is the set of all points U y OU, where U is a unitary operator. We emphasize that this is a proper subset of the unit sphere in Krylov space which, in contrast, is the set of all points UO, where U is a unitary superoperator. The gradient we have considered is with respect to the unit sphere, and it is therefore not obvious that this gradient will ever be tangential to the unitary orbit of O. However, the gradient is indeed tangential to the unitary orbit at time zero and at all times, provided the simplicity algebra is fulfilled.
On the closure of the complexity algebra. Here we show the proof that the only possible closure of the complexity algebra introduced by 42 is given by Eq. (3). The (anti-Hermitian) operator B ¼ L þ À L À "conjugated" to the Liouvillian can be expanded in Krylov space as We note that one can establish a formal analogy with the harmonic oscillator: L plays the role of the position of the harmonic oscillator, while iB corresponds to its momentum. However, in general, the commutator between L and B is not proportional to the identity, indeed: where it is understood that b 0 has to be replaced with 0. Let us now investigate the conditions under which L, B andK form a closed algebra with respect to the operation [,]: the so-called complexity algebra 42 . This happens if and only if the commutators ½L;K and ½B;K can be written as linear combinations of the operators L, B andK themselves. These commutators can be expanded over the Krylov basis as follows: where we have defined Now, it is clear that the commutator (19) between L andK cannot contain any element of the complexity algebra other than B ¼ ∑ DÀ1 n¼0 b nþ1 ½jO nþ1 ÞðO n j À jO n ÞðO nþ1 j, while the commutator (20) can only contain L ¼ ∑ DÀ1 n¼0 b nþ1 ½jO nþ1 ÞðO n j þ jO n ÞðO nþ1 j. Moreover, the only possibility for the algebra to be closed is that the discrete function f(n) is a constant. By looking at Eq. (21), we conclude that f(n) is constant if and only if for some constants α and γ (the factors 2 are included for convenience). Again, b 0 has to be replaced with 0, so that Eq. (22) holds for n ≥ 1, while 2b 2 1 ¼ α þ 2γ. Then, the function f(n) takes the constant value f = −α/2, so that the only possible closure of the complexity algebra is given by: ½L; B ¼K; ½K; L ¼ αB; ½K; B ¼ αL: Moreover, from Eq. (22), we immediately conclude that Therefore, if α ≠ 0, the Krylov complexity is related toK by a shift. Conversely, if α = 0, there is no simple relation between the Krylov complexity and the operator K. In this case,K is proportional to the identity and the complexity algebra reduces to the Heisenberg-Weyl algebra 43 , being ½L þ ; L À ¼ γ1.
Possible scenarios under the closure of the complexity algebra. As already discussed, if L, B and their commutatorK closes an algebra, then the only possible commutation relations are given by (23). This complexity algebra is then reduced to the Heisenberg-Weyl algebra whenever α = 0. We next show that for the cases α < 0 and α > 0, the complexity algebra reduces to the SU(2) algebra and the SLð2; RÞ algebra, respectively. Let us introduce the operators J + and J − , which are defined by νJ þ ¼ L þ and νJ À ¼ L À , where ν is a strictly positive scaling parameter. We can then write L ¼ νðJ þ þ J À Þ and B ¼ νðJ þ À J À Þ. Let us also introduce the operator J 0 defined by J 0 ¼ À 1 2ν 2K . By substituting these operators into (23), one can rewrite the commutation relations as By choosing the scaling parameter such that 2ν 2 = α, we find that the algebra (23) is equivalent to ½J þ ; J À ¼ J 0 ; ½J 0 ; J ± ¼ ± J ± α < 0 SUð2Þ; ð26Þ ½J þ ; J À ¼ J 0 ; ½J 0 ; J ± ¼ ÇJ ± α > 0 SLð2; RÞ: ð27Þ What we have shown is that, whenever the simplicity hypothesis holds, then the algebra generated by L, B and their commutator can always be reduced to either SU(2), SLð2; RÞ or the Heisenberg-Weyl algebra, and for which of these it reduces to depends on the value of α.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author upon reasonable request.

Code availability
The codes generated and used during the current study are available from the corresponding author on reasonable request.