* [3] P. Appell, Exemple de mouvement d’ un point assujettià une liaison exprimé par une relation non linéaire entre les composantes de la vitesse, Rend.Circ. Mat. Palermo 32 (1911), 48–50
  * [4] G.D. Birkhoff, Dynamical systems, New York, 1927\.
  * [5] A.M. Bloch, Nonholonomic Mechanics and Control, Springer, Berlin 2003.
  * [6] S.A. Chaplygin, On the theory of motion of nonholonomic systems. Theorems on the reducing multiplier, Mat. Sb. 28 (1911), 303–314 (in Russian).
  * [7] N.G. Chetaev, Izv. Fiz.Mat. Obshch. Kazan 6 (1932), 68–71 (in Russian).
  * [8] N.G. Chetaev, On Gauss principle Izv. Fiz.Mat. Obshch. Kazan 6 (1941), 323–326 (in Russian).
  * [9] M. Favretti, Equivalence of dynamics for nonholonomic systems with transverse constraints, J. Dynam. Differential Equations 10 (1998), 511–536.
  * [10] N.M. Ferrers, Extension of Lagrange’s equations. Quart. J. of pure and applied mathematics 12 (1872), 1–5.
  * [11] F.R. Gantmacher, , Lektsi po analitisheskoi mechanic, Ed. Nauka, Moscow,  1966 (in Russian).
  * [12] X. Gracia, J. Marin–Solano, M. Muñoz–Lecanda, Some geometric aspects of variational calculus in constrained systems, Reports on Mathematical Physics 51 (2003), 127–148 .
  * [13] I.M. Gelfand and S.V. Fomin, Calculus of variations, Ed. Prentice-Hall, INC., New Jersey 1963\.
  * [14] P.A. Griffiths, Exterior differential systems and the calculus of variations, Birkhäuser Boston-Basel-Stuttgart 1983\.
  * [15] G. Hamel, Teoretische Mechanik, Berlin, 1949\.
  * [16] H. Hertz, Die Prinzipien der Mechanik in neuem Zusammenhaange dargestellt, Ges. Werke, Leipzig, Barth. 1910.
  * [17] O. Hölder, Ueber die prinzipien von Hamilton und Maupertius, Nachtichten Kön. Ges. Wissenschaften zu Gottingen Math.–Phys. Kl. (1896), 122–157.
  * [18] P.V. Kharlamov, A critique of some mathematical models of mechanical systems with differential constraints, J. Appl.Math. Mech. 54 (1992), 683–691 (in Russian)
  * [19] V.I. Kirgetov, Transpositional relations in mechanics, J. Appl.Math. Mech. 22 (1958), 682–693 (in Russian).
  * [20] D.J. Korteweg, Über eine ziemlich verbreitete unrichtige Behand–lungsweise eines. Problemes der rollenden Bewegung, Nieuw Archiv voor Wiskunde 4 (1899), 130–155.
  * [21] V.V. Kozlov, Theory of integration of equations of nonholonomic mechanics, Uspekhi mekh. 8 (1985), 85–101.
  * [22] V.V. Kozlov, Realization of nonintegrable constraints in classical mechanics, Dokl. Akad. Nauk SSSR 272 (1983), 550–554 (in Russian).
  * [23] V.V. Kozlov, Gauss principle and realization of the constraints, Regular and Chaotic Dynamics 13, (2008), 431–434.
  * [24] V.V. Kozlov, Dynamics of systems with non-integrable restrictions I, Vestn. Mosk. Univ., Ser.I Mat. Mekh. 3 (1982), 92–100 (in Russian).
  * [25] V.V. Kozlov, Dynamics of systems with non-integrable restrictions II, Vestn. Mosk. Univ., Ser.I Mat. Mekh. 4 (1982), 70–76 (in Russian).
  * [26] V.V. Kozlov, Dynamics of systems with non-integrable restrictions III, Vestn. Mosk. Univ., Ser.3 Mat. Mekh. 3 (1983), 102–111 (in Russian).
  * [27] I. Kupka and W.M. Oliva The nonholonomic mechanics, J.Differential equations 169 (2001), 169–189.
  * [28] M. de Leon and D. M. de Diego, On the geometry of generalized Chaplygin systems, Mathemathical Proceeding of the Cambridge Philosophical Society 132 (2002) 1389–1412.
  * [29] A.D. Lewis and R.M. Murray, Variational principle for constrained mechanical systems: Theory and experiments,Internat. J. Non–linear Mech. 30 (1995), 793–815.
  * [30] J. Llibre, R. Ramírez and N. Sadovskaia, Integrability of the constrained rigid body, preprint, (2012).
  * [31] J. Llibre, R. Ramírez and N. Sadovskaia, Inverse problems in ordinary differential equations, preprint, (2012).
  * [32] A.I. Lurie, Analytical dynamics, Ed. Fisiko–matematisheskoi literatury, 1961\.
  * [33] C.M. Marle, Various approaches to conservative and nonconservative nonholonomic systems, Reports on Math. Physics 42 (1998), 211–229.
  * [34] J.M. Marushin, A.M. Bloch, J.E. Marsden and D.V. Zenkov, A fiber bundle approach to the transpositional relations in nonholonomic mechanics, J. of Nonlinear Sci. 22 (2012), 431–461.
  * [35] Ju.I. Neimark and N.A. Fufaev,  Dynamics of Nonholonomic Systems, American Mathematical Society, Rhode Island, 1972\.
  * [36] V.S.Novoselov, Example of a nonlinear nonholonomic constraints that is not of the type of N.G. Chetaev, Vestnik Leningrad Univ., 12 (1957) (in Russian).
  * [37] W.Muniz Oliva, Geometric mechanics, Springer–Verlag, 2002.
  * [38] L.A. Pars, A treatise on analytical dynamics, Heinemannn, London, 1968.
  * [39] H. Poincaré, Hertz’s ideas in mechanics, in addition to H. Hertz, Die Prizipien der Mechanik in neum Zusammemhauge dargestellt, 1894.
  * [40] L.S. Polak, Variation principle of mechanic, Ed. Fisico-matematicheskoi literature, 1960 (in Russian).
  * [41] R. Ramirez and N. Sadovskaia, On the dynamics of nonholonomic systems, Reports on Math. Physics. 60 (2007), 427–451.
  * [42] R. Ramirez, Dynamics of nonholonomic systems, Publisher VINITI 3878 (1985) (in Russian).
  * [43] V.N. Rubanovskii and V.A. Samsonov, Stability of steady motions, in examples and problems, M.: Nauka 1998 (in Russian).
  * [44] V.V. Rumiansev, O principe Hamiltona dlia niegolonomnix system, J. Appl.Math. Mech. 42 (1978), 407–419 (in Russian).
  * [45] N. Sadovskaia, Inverse problem in theory of ordinary differential equations, Thesis Ph. D., Univ. Politécnica de Cataluña, 2002 (in Spanish).
  * [46] J. L. Synge, On the geometry of dynamics, Phil. Trans. Roy. Soc. London ser. A 226 (1927), 31–106.
  * [47] A.S. Sumbatov, Nonholonomic systems, Regular and chaotic dynamics 7 (2002), 221–238.
  * [48] G.K.Suslov, On a particular variant of d’Alembert principle, Math. Sb. 22 (1901), 687–691 (in Russian).
  * [49] G. Zampieri, Nonholonomic versus vakonomic dynamics, J. Differential Equations 163 (2000), 335–347.
  * [50] A.M. Vershik and L.D. Faddeev, Differential geometry and Lagrangian mechanics with constraints, Soviet Physics–Doklady  17, (1972) (in Russian).
  * [51] A. Vierkandt, Über gleitende und rollende, Bewegung Monatshefte der Mathh. und Phys. III (1982) 31–54.
  * [52] V. Volterra,Sopra una classe di equazione dinamiche, Atti Accad. Sci. Torino 33 (1898), 451–475.
  * [53] P. Voronets, On the equations of motion for nonholonomic systems Math. Sb. 22 (1901), 659–686 (in Russian).


# Numerically Probing the Universal Operator Growth Hypothesis

Robin Heveling<EMAIL_ADDRESS>Department of Physics, University of
Osnabrück, D-49076 Osnabrück, Germany  Jiaozi Wang<EMAIL_ADDRESS>Department
of Physics, University of Osnabrück, D-49076 Osnabrück, Germany  Jochen Gemmer
<EMAIL_ADDRESS>Department of Physics, University of Osnabrück, D-49076
Osnabrück, Germany

###### Abstract

Recently, a hypothesis on the complexity growth of unitarily evolving
operators was presented. This hypothesis states that in generic, non-
integrable many-body systems the so-called Lanczos coefficients associated
with an autocorrelation function grow asymptotically linear, with a
logarithmic correction in one-dimensional systems. In contrast, the growth is
expected to be slower in integrable or free models. In the paper at hand, we
numerically test this hypothesis for a variety of exemplary systems, including
1d and 2d Ising models as well as 1d Heisenberg models. While we find the
hypothesis to be practically fulfilled for all considered Ising models, the
onset of the hypothesized universal behavior could not be observed in the
attainable numerical data for the Heisenberg model. The proposed linear bound
on operator growth eventually stems from geometric arguments involving the
locality of the Hamiltonian as well as the lattice configuration. We
investigate such a geometric bound and find that it is not sharply achieved
for any considered model.

##  I Introduction

The issue of the emergence of irreversible behavior from the unitary time
evolution of quantum mechanics has yet to be answered in a satisfying manner
Gogolin and Eisert (2016). In this context, concepts like the “eigenstate
thermalization hypothesis” Deutsch (1991); Srednicki (1994); D’Alessio  _et
al._ (2016) and “quantum typicality” Lloyd (2013); Goldstein  _et al._ (2006);
Reimann (2007) have been introduced as possible fundamental mechanism behind
an eventual equilibration of isolated quantum systems. The idea of typicality
is that an overwhelming majority of pure states (at some energy) give rise to
corresponding thermal expectation values. Thus, it is quite likely that over
the course of time a pure state eventually ends up in the giant “bubble” of
typical states, signaling an apparent equilibration of the system. In the
Heisenberg picture formulation of quantum mechanics, not the states are time-
dependent, but the observables themselves. Hence, it may be somewhat expected
to find a similar notion of typicality for observables, going from initially
simple, few-particle operators to more complex, generic operators. Recent
works have studied this notion of operator growth from various angles von
Keyserlingk  _et al._ (2018); Nahum  _et al._ (2018); Rakovszky  _et al._
(2018); Khemani  _et al._ (2018); Gopalakrishnan  _et al._ (2018); Chan  _et
al._ (2018).
In this paper, we refer to the particular work presented in Ref. Parker  _et
al._ (2019), in which a hypothesis on the universality of operator growth is
brought forth. Said hypothesis is formulated in the framework of the recursion
method Mattis (1981); Viswanath and Müller (2008) and makes a statement about
the growth of the so-called Lanczos coefficients, real numbers that
characterize the complexity growth of operators over the course of time. In
the following, we numerically test this operator growth hypothesis for various
models and observables. Further, we attempt to put the informal version of the
hypothesis into a more quantitative context with regard to the functional form
of the hypothesized universal growth pattern.

The paper at hand is organized as follows: we briefly restate the universal
operator growth hypothesis and introduce related quantities in Sec. II.
Following, in Sec. III, we derive an upper bound on the complexity growth of
operators based on geometric arguments. We present our numerical results and
relate them to the operator growth hypothesis in Sec. IV. We summarize our
main results and conclude in Sec. V.

##  II Operator growth hypothesis

For self-containedness, in this section, we restate the operator growth
hypothesis brought forward in Ref. Parker  _et al._ (2019). To start, the main
quantities that eventually play a role in the hypothesis are introduced. We
consider a system in the thermodynamic limit described by a local Hamiltonian
$\mathcal{H}$ [here, local means short-range, few-body interactions]. An
observable of interest represented by a Hermitian operator $\mathcal{O}$ gives
rise to a corresponding autocorrelation function

$C(t)=\text{Tr}[\mathcal{O}(t)\mathcal{O}]\,,$ (1)

where
$\mathcal{O}(t)=e^{\mathrm{i}\mathcal{H}t}\mathcal{O}e^{-\mathrm{i}\mathcal{H}t}$
is the time-dependent operator in the Heisenberg picture ($\hbar=1$). In the
following, it is convenient to work directly in the Hilbert space of operators
and denote its elements $\mathcal{O}$ as states $|\mathcal{O})$. This Hilbert
space of operators is equipped with an inner product
$(\mathcal{O}_{1}|\mathcal{O}_{2})=\text{Tr}[\mathcal{O}_{1}^{\dagger}\mathcal{O}_{2}]$,
which induces a norm via $||\mathcal{O}||=\sqrt{(\mathcal{O}|\mathcal{O})}$.
The Liouvillian superoperator is defined by
$\mathcal{L}|\mathcal{O})=[\mathcal{H},\mathcal{O}]$ and propagates a state
$|\mathcal{O})$ in time such that the autocorrelation function may be written
as $C(t)=(\mathcal{O}|e^{\mathrm{i}\mathcal{L}t}|\mathcal{O})$.
The Lanczos algorithm can be employed to calculate a tridiagonal
representation of the Liouvillian $\mathcal{L}$ in a (finite) subspace
determined by some “seed” observable $\mathcal{O}$. To start the iterative
scheme detailed below, we take the normalized initial state
$|\mathcal{O}_{0})=|\mathcal{O})$, i.e., $(\mathcal{O}|\mathcal{O})=1$, and
set $b_{1}=||\mathcal{L}\mathcal{O}_{0}||$ as well as
$|\mathcal{O}_{1})=\mathcal{L}|\mathcal{O}_{0})/b_{1}$. Then we iteratively
compute

$\displaystyle|\mathcal{Q}_{n})$
$\displaystyle=\mathcal{L}|\mathcal{O}_{n-1})-b_{n-1}|\mathcal{O}_{n-2})\,,$
(2) $\displaystyle b_{n}$ $\displaystyle=||\mathcal{Q}_{n}||\,,$
$\displaystyle|\mathcal{O}_{n})$ $\displaystyle=|\mathcal{Q}_{n})/b_{n}\,.$

The tridiagonal representation of the Liouvillian in the Krylov basis
$\\{|\mathcal{O}_{n})\\}$ is then given by

$L_{mn}=(\mathcal{O}_{m}|\mathcal{L}|\mathcal{O}_{n})=\begin{pmatrix}0&b_{1}\vphantom{\vdots}&0&...\\\
\vphantom{...}b_{1}&0\vphantom{\vdots}&b_{2}&\\\ 0&b_{2}&0&\ddots\\\
\vdots&&\ddots&\ddots\end{pmatrix}\,,$ (3)

where the Lanczos coefficients $b_{n}$ are real, positive numbers output by
the algorithm. They can be interpreted as hopping amplitudes in a tight-
binding model and their iterative computation is an elementary part of the
recursion method Mattis (1981); Viswanath and Müller (2008).
Before the hypothesis itself is stated, we will briefly present the relation
between the Lanczos coefficients $b_{n}$ and the autocorrelation function
$C(t)$ or, respectively, its Fourier transform the spectral function

$\Phi(\omega)=\int_{-\infty}^{\infty}e^{-\mathrm{i}\omega
t}\,C(t)\,\text{d}t\,.$ (4)

There exists a (non-linear) one-to-one map between the Lanczos coefficients
$b_{n}$ and the spectral function $\Phi(\omega)$, thus, a set of $b_{n}$’s
uniquely determines $\Phi(\omega)$ and vice versa. It can be shown that the
Lanczos coefficients $b_{n}$ appear in the continued fraction expansion of
$\Phi(\omega)$, i.e.,

$\Phi(\omega)=\text{Re
}\dfrac{1}{\mathrm{i}\omega+\dfrac{b_{1}^{2}}{\mathrm{i}\omega+\dfrac{b_{2}^{2}}{\mathrm{i}\omega+...}}}\,.$
(5)

The universal operator growth hypothesis brought forward in Ref. Parker  _et
al._ (2019) concerns the asymptotic behavior of the Lanczos coefficients
$b_{n}$. The hypothesis can informally be stated as follows: The Lanczos
coefficients $b_{n}$ should “grow as fast as possible” in non-integrable
systems. It turns out that [as detailed below] the fastest possible growth
rate is (asymptotically) linear, i.e.,

$b_{n}\sim\alpha n+\gamma+o(1)$ (6)

for some real constants $\alpha>0$ and $\gamma$. In the special case of a one-
dimensional system, the fastest possible growth is sub-linear due to an
additional logarithmic correction, i.e.,

$b_{n}\sim A\,\dfrac{n}{\ln n}+o(n/\ln n)\,,$ (7)

where $A>0$ is a real constant and $o(g_{n})$ denotes some real sequence
$f_{n}$ with $\lim_{n\rightarrow\infty}|f_{n}/g_{n}|=0$.
These bounds on fastest possible (asymptotic) growth eventually originate from
a powerful statement on the behavior of the spectral function $\Phi(\omega)$
for large $\omega$. The spectral function usually features non-vanishing high-
frequency tails for generic many-body systems. By means of geometric arguments
these tails can be rigorously bounded by an exponential function such that

$\Phi(\omega)\leq Ke^{-\kappa|\omega|}\,,$ (8)

for some adequately chosen constant $K>0$ and decay constant $\kappa>0$, which
is related to the geometry of the system Abanin  _et al._ (2015). It can be
shown that spectral functions actually featuring exponentially decaying tails
give rise to asymptotically linear growth in the Lanczos coefficients Lubinsky
(1993); Magnus (1987). Therefore, the operator growth hypothesis is equivalent
to an exponentially decaying spectral function and basically states that the
Lanczos coefficients should grow as fast as “permitted by the geometry”. A few
examples are known for which linear growth is analytically shown to be
achieved Parker  _et al._ (2019); Cao (2021).

##  III Bound on growth via moments

The asymptotically linear bound on the growth of the $b_{n}$ or, respectively,
the exponential bound on the decay of the spectral function $\Phi(\omega)$ are
ultimately a consequence of geometric arguments concerning the locality of the
Hamiltonian and the observable as well as the specific lattice geometry of the
system Abanin  _et al._ (2015). A straightforward way to apply these arguments
is by considering the moments $\mu_{2n}$ of the autocorrelation function and
determining an upper bound on these moments by taking the respective geometry
of the system into account. The moments of the autocorrelation function $C(t)$
are defined by

$\mu_{2n}=\dfrac{\text{d}^{2n}}{\text{d}t^{2n}}C(t)\big{|}_{t=0}$ (9)

or, respectively, in terms of the spectral function

$\mu_{2n}=\int\omega^{2n}\,\Phi(\omega)\,\text{d}\omega\,.$ (10)

Since $C(t)$ is an even function, all odd moments vanish. The information
contained in the moments $\mu_{2n}$ is identical to the information conveyed
by the Lanczos coefficients $b_{n}$. It is detailed in App. A how to translate
between the two quantities.
Employing the Heisenberg equation of motion for time-dependent operators, Eq.
(9) can be written as

$\mu_{2n}=||\mathcal{L}^{n}\mathcal{O}||^{2}\,.$ (11)

This quantity will be bounded from above in the following. We again consider a
local operator $\mathcal{O}$ with $||\mathcal{O}||=1$ and a local Hamiltonian
$\mathcal{H}=\sum_{\ell}h_{\ell}$ with local terms $h_{\ell}$ or,
respectively, a local Liouvillian $\mathcal{L}=\sum_{k}\ell_{k}$ with local
terms $\ell_{k}=[h_{k},\cdot\,]$.
The norm of a local Liovillian $\ell$ [we assume periodicity such that all
$\ell_{k}$ are of the same type] applied to some operator $\mathcal{A}$ can be
bounded by

$||\ell\mathcal{A}||\leq\mathcal{E}\,||\mathcal{A}||\,,$ (12)

where $\mathcal{E}=E_{\text{max}}-E_{\text{min}}$ denotes the maximum
eigenvalue of $\ell$ and $E_{\text{max}}$ ($E_{\text{min}}$) is the maximum
(minimum) eigenvalue of the local Hamiltonian $h$. Equality holds if the
operator $\mathcal{A}$ is an eigenoperator of the local Liouvillian
corresponding to the largest eigenvalue. Utilizing the triangle inequality and
iteratively applying Eq. (12) yields

$\displaystyle||\mathcal{L}^{n}\mathcal{O}||$
$\displaystyle=||\sum_{k_{1},...,k_{n}}\ell_{k_{n}}...\,\ell_{k_{1}}\mathcal{O}||$
$\displaystyle\leq\sum_{k_{1},...,k_{n}}||\ell_{k_{n}}...\,\ell_{k_{1}}\mathcal{O}||$
$\displaystyle\leq\sum_{k_{1},...,k_{n}}\mathcal{E}^{n}=\mathcal{E}^{n}N_{\text{sum}}(n)\,,$
(13)

where $N_{\text{sum}}$ denotes the number of terms in the sum [it is specified
below which terms are actually counted]. This number typically grows quite
fast with $n$ and can be exactly determined for simple geometries, as is
presented at the end of this section. Consequently, the moments can be bounded
as

$\mu_{2n}=||\mathcal{L}^{n}\mathcal{O}||^{2}\leq\mathcal{E}^{2n}N_{\text{sum}}^{2}(n)\,.$
(14)

This bound on the moments is sharp, meaning first, no sequence of moments
$\mu_{2n}$ calculated via Eq. (9) can possibly grow faster and, importantly,
second, that this bound can in principle be actually achieved with an equality
sign. This is the case if and only if two conditions are met. Firstly, to get
an equality in the triangle inequality, all operators occurring in the sum
must be collinear. Secondly, the largest eigenvalue $\mathcal{E}$ must be
realized at each application of a local Liouvillian. Indeed, it would be quite
surprising if this bound were to be achieved tightly for physical systems. On
the other hand, it may not be as far off as one might initially guess, since
multiple applications of the same operator [here $\mathcal{L}$] to a state
increase the overlap of the resulting state with the states at the edges of
the spectrum [given that the initial state has some overlap with eigenstates
corresponding to extremal eigenvalues].
The fastest possible growing moments, i.e., moments equal to the r.h.s. of Eq.
(14), can be translated to corresponding Lanczos coefficients as described in
App. A. For later reference, we denote the coefficients obtained this way by
$\mathcal{B}_{n}$. They depend on the energy scale $\mathcal{E}$ as well as on
the lattice geometry expressed through $N_{\text{sum}}$, both quantities are
exactly known for the models studied below. Even though Eq. (14) gives a
rigorous upper bound on the moments, the $\mathcal{B}_{n}$ obtained from this
bound do not necessarily constitute a pointwise upper bound on the $b_{n}$. In
Sec. IV.1 we give more details on the interpretation of the coefficients
$\mathcal{B}_{n}$.
To end this section, we determine the number $N_{\text{sum}}$ for the simple
case of a one-dimensional chain with nearest-neighbor interactions, i.e.,
$\mathcal{H}=\sum_{\ell}h_{\ell,\ell+1}$. We start with a local operator
$\mathcal{O}$ whose support is only on site zero [the support of an operator
contains all sites on which the operator is not equal to the identity, e.g.,
here $\mathcal{O}=...\otimes I\otimes\mathcal{O}_{0}\otimes I\otimes...$,
where $I$ denotes the identity on a given site]. The operator
$\mathcal{L}\mathcal{O}$ consists of operators $\ell_{-1,0}\mathcal{O}_{0}$
with support on sites $(-1,0)$ as well as $\ell_{0,1}\mathcal{O}_{0}$ with
support on sites $(0,1)$. Next, $\mathcal{L}^{2}\mathcal{O}$ contains six
operators with support on sites $(-2,-1,0)$, $(-1,0)$, $(-1,0,1)$, $(-1,0,1)$,
$(0,1)$, $(0,1,2)$, respectively. In these lists we include [and count]
trivially non-vanishing operators, i.e., operators that vanish due to a lack
of overlap between respective supports are not counted [for example the
operator $\ell_{8,9}\mathcal{O}_{0}$ trivially vanishes], however, operators
with respective overlap between supports [like $\ell_{0,1}\mathcal{O}_{0}$]
are always counted, even though the operator may vanish due to the specific
choice of the local Hamiltonian and initial observable. In this manner, we
iteratively apply the Liouvillian to the initial operator, grow the supports
accordingly and keep track of the number of potentially non-vanishing
operators. For the case at hand, this procedure gives rise to a sequence
$N_{\text{sum}}(n)=1,2,6,22,94,454,2430,14214,...$ of numbers of terms in the
sum in Eq. (14) [for this exact sequence there happens to exist a closed form:
$N_{\text{sum}}(n)=\sum_{k=0}^{n}2^{k}S_{n}^{(k)}$, where $S_{n}^{(k)}$
denotes the Stirling number of the second kind].
Corresponding iteratively computed sequences can be found for more complicated
local Hamiltonians, e.g., next-nearest neighbor interactions, quasi one-
dimensional systems and two-dimensional systems.

##  IV Numerical Analysis

In this section, we numerically check the proposed operator growth hypothesis
by explicitly calculating the Lanczos coefficients $b_{n}$ for various
exemplary setups. These setups include one-dimensional and two-dimensional
Ising models as well as one-dimensional Heisenberg models, all paired with a
variety of different observables. We compare the calculated $b_{n}$ with the
coefficients $\mathcal{B}_{n}$ [obtained from the r.h.s. of Eq. (14)] by
explicitly calculating $\mathcal{E}$ and $N_{\text{sum}}$ for each model.
In practice, it is only possible [for the considered systems at least] to
obtain a finite number $N$ of Lanczos coefficients $b_{n}$, since the
dimension of the operator space grows exponentially. The achievable number of
coefficients $N$ is around $30$ for the $1$d Ising model and around $15$ for
the Heisenberg model. The difference in obtainable $b_{n}$ is due to the fact
that the operator space grows much faster for the Heisenberg model than for
the Ising model. For all considered systems the Hamiltonian $\mathcal{H}$
consist of two terms, an integrable part $\mathcal{H}_{0}$ and an
integrability-breaking part $\mathcal{V}$, i.e.,

$\mathcal{H}=\mathcal{H}_{0}+\lambda\mathcal{V}\,$ (15)

FIG. 1. Lanczos coefficients $b_{n}$ of the transverse Ising model for the
$2$-local observables $a$) $\mathcal{O}^{(1)}$, $b$) slow mode of
$\mathcal{O}^{(2)}_{q}$ with $q=\pi/13$, and $c$) fast mode of
$\mathcal{O}^{(2)}_{q}$ with $q=\pi$. The integrability-breaking magnetic
field attains values from $B_{x}=0.0$ to $B_{x}=0.5$. For all observables, the
transition from a free model to a non-integrable model is evident. The
coefficients $\mathcal{B}_{n}$ are explicitly depicted in $a$) for the case
$B_{x}=0.5$ as yellow dots. Dashed lines indicate the “lower branches” of the
corresponding $\mathcal{B}_{n}$. To avoid clutter, only the dashed lines are
depicted as a guide to the eye for other values of $B_{x}$ and in all
following figures. The coefficients $\mathcal{B}_{n}$ are larger than the
physical $b_{n}$ by a factor of about two for all observables.

[except for the 2d Ising model, where $\mathcal{H}_{0}$ is already non-
integrable]. The parameter $\lambda$ tunes the non-integrability of the model.
We suppose that the total Hamiltonian $\mathcal{H}$ (as well as
$\mathcal{H}_{0}$ and $\mathcal{V}$ individually) can be written in terms of
local Hamiltonians, i.e., $\mathcal{H}=\sum_{\ell}h_{\ell}$. Again, the local
terms usually describe short-range, few-body interactions. For each model, we
consider a number of observables $\mathcal{O}$. Importantly, all observables
should have zero overlap with any conserved quantity Mazur (1969), for example
$(\mathcal{O}|\mathcal{H})=0$.

###  IV.1 Transverse Ising model

The first model under consideration is a transverse Ising model with a tilted
field. Respective unperturbed and total Hamiltonians are given by

$\displaystyle\mathcal{H}_{0}=$
$\displaystyle\sum_{\ell}J_{xx}\sigma^{x}_{\ell}\sigma^{x}_{\ell+1}+B_{z}\sigma^{z}_{\ell}$
(16) $\displaystyle\mathcal{H}$
$\displaystyle=\mathcal{H}_{0}+B_{x}\sum_{\ell}\sigma^{x}_{\ell}\,,$

where $\sigma_{\ell}^{x,y,z}$ denote Pauli operators on site $\ell$. The
magnetic field $B_{x}$ in $x$-direction plays the role of the integrability-
breaking parameter $\lambda$ in Eq. (15), i.e., the system is non-integrable
for $B_{x}\neq 0$ and integrable for $B_{x}=0$. We set $J_{xx}=1.0$ and
$B_{z}=-1.05$ and calculate the Lanczos coefficients for various observables
as detailed in Sec. II. In practice, it is convenient to adopt the set of
Pauli strings as a working basis of the Hilbert space of operators Parker  _et
al._ (2019); Dehaene and De Moor (2003). As a first example, we consider the
$2$-local observable

$\mathcal{O}^{(1)}\propto\sum_{\ell}1.05\,\sigma^{x}_{\ell}\sigma^{x}_{\ell+1}+\sigma^{z}_{\ell}\,,$
(17)

where $2$-local means that the local terms are supported on two sites
respectively. The choice of the parameter $1.05$ in front of the $xx$-coupling
term ensures that $(\mathcal{O}^{(1)}|\mathcal{H})=0$. This exact setup was
also studied in Ref. Parker  _et al._ (2019). The corresponding Lanczos
coefficients are depicted in Fig. 1a for various values of the integrability-
breaking parameter $B_{x}$.
The Hamiltonian of the bare transverse Ising model [with $B_{x}=0$] can be
mapped onto free fermions via a Jordan-Wigner transformation. Further, the
observable in Eq. (17) is local in the fermionic picture. In this non-
interacting case, the Lanczos coefficients seem to be more or less constant.
As soon as a small perturbation that breaks the integrability is introduced,
e.g., $B_{x}=0.01$, the Lanczos coefficients begin to grow. The distinction
between the free case and non-integrable cases is clearly visible in Fig. 1a.
The growth of the Lanczos coefficients for larger values of $B_{x}$ already
looks quite linear, a possible logarithmic correction due to the system’s one-
dimensionality is not visible in the data.
Before we continue, we want to make some remarks on the interpretation of the
coefficients $\mathcal{B}_{n}$. As mentioned, the $\mathcal{B}_{n}$ are
computed by assuming that the moments grow as fast as possible, i.e., an
equality sign in Eq. (14), and translating these maximum moments to Lanczos
coefficients as detailed in App. A. However, even though Eq. (14) gives a
rigorous upper bound on the moments, the resulting $\mathcal{B}_{n}$ do not
necessarily constitute a strict pointwise upper bound on the $b_{n}$ [this is
simply due to the way the Lanczos coefficients are calculated from the
moments]. Thus, the $\mathcal{B}_{n}$ should not be thought of as such.
Rather, the $\mathcal{B}_{n}$ represent a sort of “global uniform” upper
bound, meaning that it is impossible to further increase the value of one
specific coefficient “by hand”, without simultaneously decreasing the value of
another one [or several other ones]. In this sense the $\mathcal{B}_{n}$ are
the “maximum” coefficients. If the “physical” coefficients $b_{n}$ would
follow the $\mathcal{B}_{n}$ tightly, then one could indeed conclude that the
upper bound in Eq. (14) was sharply achieved and that the Lanczos coefficients
would indeed grow “as fast as possible”, given the constraints in Eqs. (12),
(III).
The values of the $\mathcal{B}_{n}$ are explicitly depicted in Fig. 1a for
$B_{x}=0.5$. The coefficients clearly exhibit some even-odd effects. These
even-odd effects also occur for all other considered parameters and models. To
avoid clutter, we only show the “lower branches” of the $\mathcal{B}_{n}$ as
dashed lines for smaller values of $B_{x}$. Therefore, the dashed lines in
Fig. 1a [and all following figures] serve as a guide to the eye for the
“maximal possible” coefficients $\mathcal{B}_{n}$. The energy scale is
$\mathcal{E}=3.2$ [for $B_{x}=0.5$] and $N_{\text{sum}}$ is obtained as
detailed at the end of Sec. III, with the important difference that the
initial operator is now supported on two sites. In Fig. 1a it is evident that
the coefficients $\mathcal{B}_{n}$ are larger than the physical coefficients
$b_{n}$ by about a factor of $2$. Thus, the upper bound on the moments is not
sharply achieved, i.e., Eq. (14) is an overestimate not . This general
behavior of the $b_{n}$ as well as the $\mathcal{B}_{n}$ is reproduced for the
next observable, which constitutes an energy density wave with momentum $q$,
i.e.,

$\mathcal{O}_{q}^{(2)}\propto\sum_{\ell}\cos(q\ell)h_{\ell}\,.$ (18)

Figure 2: Lanczos coefficients $b_{n}$ of the transverse Ising model for a
local observable $\mathcal{O}^{(3)}=\sigma_{0}^{x}$ for various integrability-
breaking parameters $B_{x}$. The distinction between the free and non-
nonintegrable curves is not as striking a before. Green dash-dotted line
indicates a fit $\propto\sqrt{n}$ to the data of the integrable model. Dashed
lines serve as a guide to the eye for the coefficients $\mathcal{B}_{n}$,
which are larger by a factor of about two.

We study a relatively slow dynamic with $q=\pi/13$ and a relatively fast
dynamic with $q=\pi$. Both observables are local in the fermionic picture. The
Lanczos coefficients $b_{n}$ and coefficients $\mathcal{B}_{n}$ are depicted
in Fig. 1b for $q=\pi/13$ and in Fig. 1c for $q=\pi$. In both cases, the
qualitative behavior is quite similar to the first observable. Again, the
Lanczos coefficients of the free model with $B_{x}=0$ seem to be more or less
constant. Once the additional magnetic field is added, the model becomes non-
integrable and at some point the $b_{n}$ grow approximately linearly. Just as
for the first observable, the derived bounds are not tight and the
$\mathcal{B}_{n}$ are larger by a factor of about two.
The final considered observable for the one-dimensional Ising model is a local
operator whose support only contains a single site, i.e.,

$\mathcal{O}^{(3)}=\sigma^{x}_{0}\,.$ (19)

The corresponding Lanczos coefficients are depicted in Fig. 2. There is a
clear qualitative difference compared to the observables investigated thus
far. For the free case with $B_{x}=0$ the Lanczos coefficients seem to no
longer be bounded by a constant. Rather, the growth is quite accurately
described by a square-root $\propto\sqrt{n}$ [see fit]. For this particular
model and observable the square-root-like growth can be understood
analytically Cao (2021). Further, this kind of growth has been observed for a
variety of other integrable models Viswanath and Müller (2008); Parker  _et
al._ (2019); Lee (2001). We suspect that this qualitatively different behavior
compared to the previous cases is due to the specific choice of the
observable, which is, in contrast to all previously considered observables,
non-local in the fermionic picture. Consequently, one could be inclined to
formulate two sufficient conditions, which both have to be met in order for
the $b_{n}$ to be bounded by a constant. First, the Hamiltonian has to
describe a free model and, second, the observable has to be local. For the
observable at hand, which is non-local in the fermionic picture, the second
condition is violated. Therefore, the Lanczos coefficients are not bounded by
a constant, rather they grow as a square-root.
As the Hamiltonian departs from the integrable/free point once $B_{x}\neq 0$,
the $b_{n}$ grow faster [which is not too surprising, since there are simply
more terms in the Hamiltonian]. From the computed data it is not immediately
obvious whether the growth becomes linear (with a logarithmic correction) or
remains more or less square-root-like. Possibly, the data for larger $B_{x}$
in Fig. 2 hints at an onset of linear growth for larger $n$. However, the
distinction is certainly not as clear as for the local observables. Without
previous knowledge of which coefficients belong to which $B_{x}$, it would be
quite difficult to say if a set of $b_{n}$ belongs to an integrable/free or
non-integrable model [barring the absolute values]. Therefore, calculating the
$b_{n}$ as a potential method to determine or define (non-)integrability does
not seem feasible, since the number at which universal behavior sets in may be
larger than the practically computable number of Lanczos coefficients. Again,
the $\mathcal{B}_{n}$ are off by a factor of about $2$.
It is interesting to note that in the non-integrable models all considered
observables seem to more or less lead to similar growth patterns and attain
similar values for larger $n$. For comparison, Fig. 3 depicts the Lanczos
coefficients for all four observables considered thus far for $B_{x}=0.5$.
This figure certainly supports the hypothesis of a universality of operator
growth brought forth in Ref. Parker  _et al._ (2019). Particularly striking is
the relation between the observable $\mathcal{O}^{(1)}$ and the slow Fourier
mode $\mathcal{O}_{q=\pi/13}^{(2)}$ since for $n\gtrsim 10$ the coefficients
practically coincide.

Figure 3: Comparison between Lanczos coefficients for the Ising model with
$B_{x}=0.5$ for all four obervables considered thus far. The growth is quite
similar for larger $n$ hinting at a universality of operator growth.

Before leaving the Ising model and continuing with the Heisenberg model, we
want to briefly present data on the 2d Ising model. As the two-dimensional
Ising model is non-integrable, the hypothesis predicts a strict asymptotically
linear growth [without logarithmic correction] of the Lanczos coefficients
$b_{n}$. Respective Hamiltonians of the two-dimensional Ising model are given
by

$\displaystyle\mathcal{H}_{0}=\sum_{\ell,\ell^{\prime}}J_{xx}$
$\displaystyle\sigma^{x}_{\ell,\ell^{\prime}}\sigma^{x}_{\ell+1,\ell^{\prime}}+J_{xx}^{\prime}\sigma^{x}_{\ell,\ell^{\prime}}\sigma^{x}_{\ell,\ell^{\prime}+1}+B_{z}\sigma^{z}_{\ell,\ell^{\prime}}$
$\displaystyle\mathcal{H}=\mathcal{H}_{0}+B_{x}\sum_{\ell,\ell^{\prime}}\sigma^{x}_{\ell,\ell^{\prime}}\,,$
(20)

where primed indices number the vertical direction and unprimed indices the
horizontal direction. As in the one-dimensional case, we set $J_{xx}=1.0$,
$B_{z}=-1.05$ and vary $B_{x}$. The coupling strength in vertical direction is
set to $J^{\prime}_{xx}=0.5$. Note that this model in non-integrable for all
values of $B_{x}$. The energy scale is given by $\mathcal{E}=3.9$ [for
$B_{x}=0.5$]. We again consider a local observable whose support is only on
one site, i.e.,

$\mathcal{O}=\sigma_{0,0}^{x}\,.$ (21)

The calculated Lanczos coefficients $b_{n}$ are depicted in Fig. 4. Since the
space of operators grows extremely fast, it is only practically possible to
calculate about $13$ coefficients. These coefficients grow in a nicely linear
fashion for all values of $B_{x}$, which is in accord with the operator growth
hypothesis. The coefficients $\mathcal{B}_{n}$ from the derived upper bound on
the moments are very far off. So far off, in fact, that there is an additional
vertical axis for the coefficients $\mathcal{B}_{n}$ on the right side in Fig.
4, which includes a factor of four. As mentioned, the number $N_{\text{sum}}$
is computed iteratively. For the two-dimensional Ising model the operator
space grows so fast that the $\mathcal{B}_{n}$ are only attainable up to
$n=5$.

Summarizing the results from this section, the operator growth hypothesis is
supported by (most of) the numerical data. The Lanczos coefficients of the
one-dimensional non-integrable Ising models seem to eventually attain
approximate linear growth for observables that are local in the fermionic
picture. In these cases, the transition from free to non-integrable is clearly
visible. A possible logarithmic correction is not noticeable in the presented
data. Such a correction would probably only be visible in data that spans
multiple orders of magnitude, which is not the case here. Further, the data
for the two-dimensional Ising model supports the hypothesis for all considered
values of $B_{x}$. Only the data for the third observable
$\mathcal{O}^{(3)}=\sigma_{0}^{x}$ is somewhat inconclusive. The transition is
not as distinct as for the other observables, however, the onset of the
hypothesized universal behavior may be suspected for larger $n$. All these
numerical results support the in Ref. Parker  _et al._ (2019) proposed
operator growth hypothesis in the sense that the $b_{n}$ grow asymptotically
linear. This is the first main result of the paper at hand.
In Ref. Parker  _et al._ (2019), the operator growth hypothesis is several
times informally stated as: the Lanczos coefficients should “grow as fast as
possible” and therefore eventually attain said linear growth. In fact, the
data presented in this section corroborates the notion of asymptotically
linear growth. However, the Lanczos coefficients could, in principle, grow
much faster, as indicated by the coefficients $\mathcal{B}_{n}$, which are not
tightly achieved in any of the considered models. Thus, we want to clarify the
“fastest possible growth” is not to be understood with respect to the absolute
numerical values of the coefficients, but rather in regard to the functional
form of their growth. Indeed, the coefficients $b_{n}$ seem to grow with quite
a similar functional form as the [much to large] coefficients
$\mathcal{B}_{n}$. This is the second main result of the paper at hand.

Figure 4: Lanczos coefficients $b_{n}$ of the two-dimensional transverse Ising
model for an observable $\mathcal{O}=\sigma_{0,0}^{x}$ for various $B_{x}$.
For all values of $B_{x}$ the growth is nicely linear. Dashed lines serve as a
guide to eye for the coefficients $\mathcal{B}_{n}$, which are much larger
(note the additional vertical axis on the right).

FIG. 5. Lanczos coefficients $b_{n}$ of the Heisenberg model for spin density
waves. Depicted are various combinations of anisotropies and momenta: $a$)
$\Delta=0.5$, $q=\pi/13$; $b$) $\Delta=0.5$, $q=\pi$; $c$) $\Delta=1.5$,
$q=\pi/13$; $d$) $\Delta=1.5$, $q=\pi$. The integrability-breaking parameter
$\Delta^{\prime}$ attains values from $0.0$ to $0.5$. Dashed lines serve as a
guide to the eye for the coefficients $\mathcal{B}_{n}$, which are much larger
(note the additional vertical axes on the right).

###  IV.2 Heisenberg model

The second model of interest is a Heisenberg model with an additional next-
nearest-neighbor interaction. Respective Hamiltonians are given by

$\displaystyle\mathcal{H}_{0}=\sum_{\ell}$
$\displaystyle\sigma^{x}_{\ell}\sigma^{x}_{\ell+1}+\sigma^{y}_{\ell}\sigma^{y}_{\ell+1}+\Delta\sigma^{z}_{\ell}\sigma^{z}_{\ell+1}$
(22) $\displaystyle\mathcal{H}$
$\displaystyle=\mathcal{H}_{0}+\Delta^{\prime}\sum_{\ell,\ell^{\prime}}\sigma^{z}_{\ell}\sigma^{z}_{\ell+2}\,.$

The anisotropy of the nearest-neighbor interaction is denoted by $\Delta$ and
the integrability-breaking next-nearest-neighbor interaction is tuned by the
parameter $\Delta^{\prime}$ [which plays the role of $\lambda$ in Eq. (15)].
The bare Heisenberg chain [with $\Delta^{\prime}=0$] is integrable for any
anisotropy $\Delta$.
The model is gapless and exhibits ballistic transport behavior [of spin and
energy] for $|\Delta|<1$, whereas for $|\Delta|>1$ the transport of spin is
diffusive, while energy transport is still ballistic Bertini  _et al._ (2021).
In the following numerics we cover both cases by choosing values
$\Delta=0.5,1.5$. The parameter $\Delta^{\prime}$ that breaks the
integrability is varied in the same fashion as $B_{x}$ in the Ising model.
Note that the full Hamiltonian in Eq. (22) conserves the total magnetization
in $z$-direction. Therefore, it is natural to consider a spin density wave
with momentum $q$, i.e.,

$\mathcal{O}_{q}\propto\sum_{\ell}\cos(q\ell)\sigma^{z}_{\ell}\,.$ (23)

Similar as for the energy density wave in the Ising model, we study a
relatively slow dynamic with $q=\pi/13$ [depicted in Figs. 5a, 5c for
$\Delta=0.5,1.5$ respectively] and a relatively fast dynamic with $q=\pi$
[depicted in Figs. 5b, 5d for $\Delta=0.5,1.5$ respectively].
Since the Heisenberg Hamiltonian contains more coupling terms than the Ising
Hamiltonian, the dimension of the operator space [with respect to the Pauli
basis] grows faster and we are restricted to a smaller number of coefficients,
only about $15$.
Comparing the variance (relative deviations) of the Lanczos coefficients for
$\Delta=0.5$ and $\Delta=1.5$, it is striking that the coefficients for
$\Delta=1.5$ vary much less for different values of $\Delta^{\prime}$. This is
most likely due to the relative strength of the perturbation. Let
$||\mathcal{H}_{0}^{\Delta}||$ denote the norm of the unperturbed Hamiltonian
with anisotropy $\Delta$ and $||\lambda\mathcal{V}||$ the strength of the
perturbation [where $\lambda$ conforms to $\Delta^{\prime}$]. Then, for
example $||0.5\mathcal{V}||/||\mathcal{H}_{0}^{1.5}||=0.34$ and
$||0.5\mathcal{V}||/||\mathcal{H}_{0}^{0.5}||=0.47$ [for comparison, in the
one-dimensional Ising model $||0.5\mathcal{V}||/||\mathcal{H}_{0}||=0.28$].
Again, the stronger the perturbation, the faster the coefficients grow, which
is simply due to additional terms in the Hamiltonian [in order to not obscure
the main points, we refrain from rescaling the Hamiltonian accordingly]. Other
than that, the growth is more irregular than in the Ising model, at least in
the regime where data is available. For both values of $\Delta$ the transition
occurs between an integrable [$\Delta^{\prime}=0$] and non-integrable
[$\Delta^{\prime}\neq 0$] model. However, for $\Delta=0.5$ there is neither
square-root-like growth in the integrable case nor linear growth in the non-
integrable case visible [for both modes with $q=\pi/13$ in Fig. 5a and $q=\pi$
in Fig. 5b]. For $\Delta=1.5$, $\Delta^{\prime}=0.0$ and $q=\pi/13$ the growth
of the coefficients is more similar to a square-root [see fit] and only
visibly deviates for $n\gtrsim 12$, see Fig. 5c. For the faster mode with
$q=\pi$ the growth seems more linear with relatively small deviations, see
Fig. 5d.
The coefficients $\mathcal{B}_{n}$ are much larger than any of the $b_{n}$
such that the additional vertical axes contain a factor of five in all cases.
The number of terms in the sum $N_{\text{sum}}$ grows quite a lot faster than
in the Ising model due to the next-nearest-neighbor interaction. These terms
must in principle be included as soon as $\Delta^{\prime}$ attains an
arbitrarily small value strictly greater than zero. This leads to a somewhat
gross overestimation, since the energy scale $\mathcal{E}$ remains basically
unaltered for small $\Delta^{\prime}$. This is visible in all figures for the
Heisenberg model, where the coefficients $\mathcal{B}_{n}$ for
$\Delta^{\prime}=0$ are calculated with the smaller numbers $N_{\text{sum}}$
of nearest-neighbor interaction. In principle, one could improve the upper
bound on the moments by introducing a second energy scale
$\mathcal{E}_{\text{nnn}}$ of the next-nearest-neighbor interaction and count
terms according to the appearance of nearest-neighbor terms $\ell_{k,k+1}$ and
next-nearest-neighbor terms $\ell_{k,k+2}$. This is, however, more complicated
and not in the spirit of the derivations presented in Refs. Parker  _et al._
(2019); Abanin  _et al._ (2015).
Summarizing, the numerical data for all considered values of $\Delta$ and $q$
can neither really reject nor support the operator growth hypothesis [not
least because data for larger $n$ is not available]. The transitions between
integrable and non-integrable models are certainly less striking than for the
Ising model. Again, only looking at the $b_{n}$ it would be impossible to say
which coefficients belong to an integrable or non-integrable model. This
possibly raises the question whether the distinction between integrability and
non-integrability regarding the growth of the $b_{n}$ is an adequate
distinction to make. As mentioned in Sec. II, the Lanczos coefficients are
uniquely determined by the autocorrelation function $C(t)$. The Heisenberg
model for $\Delta=1.5$ and $\Delta^{\prime}=0$ is integrable and exhibits
diffusive transport behavior Bertini  _et al._ (2021), which is usually
attributed to chaotic, non-integrable systems. In view of this, it may be not
too surprising that the operator growth hypothesis is not supported by the
(limited) numerical data presented in this section. This is the third main
result of the paper at hand.

##  V Conclusion

The first main message of the paper at hand concerns the explicitly calculated
numerical data on the Lanczos coefficients. We numerically probed the
universal operator growth hypothesis proposed in Ref. Parker  _et al._ (2019),
which states that in generic, non-integrable systems the Lanczos coefficients
grow asymptotically linear [with a logarithmic correction in 1d]. We
explicitly calculated Lanczos coefficients $b_{n}$ for various combinations of
models [including 1d and 2d Ising models as well as ballistic and diffusive
Heisenberg models] and observables [including energy and spin density waves as
well as local observables]. We found that the Ising model data generally
supports the operator growth hypothesis. In particular, as soon as an
integrability-breaking perturbation is added to the Hamiltonian, the
coefficients eventually attain a linear growth [this transition is more
pronounced in the case of a free Hamiltonian with a local observable (in the
fermionic picture) than in the case of a free Hamiltonian with a non-local
observable (in the fermionic picture)]. Further, the two-dimensional Ising
model exhibits clear linear growth. Inconclusive, however, remains the data
for the Heisenberg model. For none of the combinations of considered
parameters there is a clear distinction between the integrable and non-
integrable cases. Of course, it may be possible that the hypothesized
universal behavior only sets in at some larger $n$, which is not accessible by
our numerical tools.
The second main message of the paper at hand concerns the coefficients
$\mathcal{B}_{n}$, which are obtained by considering the fastest growing
moments and converting them into Lanczos coefficients. The informal version of
the operator growth hypothesis is stated several times in Ref. Parker  _et
al._ (2019), namely that “the Lanczos coefficients should grow as fast as
possible” and a corresponding bound on the moments leading to linear growth
[with a logarithmic correction in one dimension] is given not . As seen in the
available numerical data, even the “optimized” bound in Eq. (14) is not
remotely tight and the “physical” Lanczos coefficients $b_{n}$ increase much
slower than the “fastest possible growing” coefficients $\mathcal{B}_{n}$ in
all considered models. Therefore, technically, the Lanczos coefficients do not
grow as fast as possible. Nevertheless, the physical Lanczos coefficients
$b_{n}$ seem to grow in a manner that is compatible with the “functional form”
of the maximal growth, i.e., we observe more or less linear growth for the
Ising models, only with a flatter slope than would be induced by the bound on
the moments.

## Acknowledgments

We thank D. Parker for interesting discussion on this topic and for a comment
on an earlier version of this paper. This work was supported by the Deutsche
Forschungsgemeinschaft (DFG) within the Research Unit FOR 2692 under Grant No.
397107022 (GE 1657/3-2).

## References

  * Gogolin and Eisert (2016) C. Gogolin and J. Eisert, “Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems,” Reports on Progress in Physics 79, 056001 (2016).
  * Deutsch (1991) J. M. Deutsch, “Quantum statistical mechanics in a closed system,” Phys. Rev. A 43, 2046–2049 (1991).
  * Srednicki (1994) M. Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E 50, 888–901 (1994).
  * D’Alessio  _et al._ (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics 65, 239–362 (2016).
  * Lloyd (2013) S. Lloyd, “Pure state quantum statistical mechanics and black holes,” (2013), arXiv:1307.0378 .
  * Goldstein  _et al._ (2006) S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghì, “Canonical Typicality,” Phys. Rev. Lett. 96, 050403 (2006).
  * Reimann (2007) P. Reimann, “Typicality for Generalized Microcanonical Ensembles,” Phys. Rev. Lett. 99, 160404 (2007).
  * von Keyserlingk  _et al._ (2018) C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, and S. L. Sondhi, “Operator Hydrodynamics, OTOCs, and Entanglement Growth in Systems without Conservation Laws,” Phys. Rev. X 8, 021013 (2018).
  * Nahum  _et al._ (2018) A. Nahum, S. Vijay, and J. Haah, “Operator Spreading in Random Unitary Circuits,” Phys. Rev. X 8, 021014 (2018).
  * Rakovszky  _et al._ (2018) T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk, “Diffusive Hydrodynamics of Out-of-Time-Ordered Correlators with Charge Conservation,” Phys. Rev. X 8, 031058 (2018).
  * Khemani  _et al._ (2018) V. Khemani, A. Vishwanath, and D. A. Huse, “Operator Spreading and the Emergence of Dissipative Hydrodynamics under Unitary Evolution with Conservation Laws,” Phys. Rev. X 8, 031057 (2018).
  * Gopalakrishnan  _et al._ (2018) S. Gopalakrishnan, D. A. Huse, V. Khemani, and R. Vasseur, “Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems,” Phys. Rev. B 98, 220303 (2018).
  * Chan  _et al._ (2018) A. Chan, A. De Luca, and J. T. Chalker, “Solution of a Minimal Model for Many-Body Quantum Chaos,” Phys. Rev. X 8, 041019 (2018).
  * Parker  _et al._ (2019) D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, “A Universal Operator Growth Hypothesis,” Phys. Rev. X 9, 041017 (2019).
  * Mattis (1981) Daniel C. Mattis, “How to reduce practically any problem to one dimension,” in  _Physics in One Dimension_ (Springer Berlin Heidelberg, 1981) pp. 3–10.
  * Viswanath and Müller (2008) V.S. Viswanath and G. Müller,  _The Recursion Method: Applications to Many-Body Dynamics_ (Springer, New York, 2008).
  * Abanin  _et al._ (2015) D. A. Abanin, W. De Roeck, and F. Huveneers, “Exponentially Slow Heating in Periodically Driven Many-Body Systems,” Phys. Rev. Lett. 115, 256803 (2015).
  * Lubinsky (1993) D. S. Lubinsky, “An update on orthogonal polynomials and weighted approximation on the real line,” Acta Applicandae Mathematica 33, 121–164 (1993).
  * Magnus (1987) A. Magnus, “Asymptotic behaviour of continued fraction coefficients related to singularities of the weight function,” in  _The Recursion Method and Its Applications_ (Springer Berlin Heidelberg, 1987) pp. 22–45.
  * Cao (2021) X. Cao, “A statistical mechanism for operator growth,” Journal of Physics A: Mathematical and Theoretical 54, 144001 (2021).
  * Mazur (1969) P. Mazur, “Non-ergodicity of phase functions in certain systems,” Physica 43, 533–545 (1969).
  * Dehaene and De Moor (2003) J. Dehaene and B. De Moor, “Clifford group, stabilizer states, and linear and quadratic operations over GF(2),” Phys. Rev. A 68, 042318 (2003).
  * (23) In Ref. Parker  _et al._ (2019), App. F a different upper bound on the moments [similar to the bound in Eq. (14)] is derived, i.e., $\mu_{2n}\leq{P}_{n}^{2}(2\mathcal{E}^{\prime})^{2n}$, where ${P}_{n}\leq n4^{n}\max_{\ell\in[0,n]}\ell^{n-\ell}$ and $\mathcal{E}^{\prime}=||h||$. It is shown that this bound corresponds to asymptotically linear growth with a logarithmic correction. Since the spectral function ${\Phi}(\omega)$ attains negative values for moments equal to the r.h.s. of the above bound, it is not possible to obtain corresponding $\mathcal{B}_{n}$. However, the numbers $P_{n}$ grow rapidly, much quicker than the sequence of numbers found at the end of Sec. III, further, $2\mathcal{E}^{\prime}$ is larger than $\mathcal{E}$. Thus, the above bound on the moments is looser than the derived bound in Eq. (14) .
  * Lee (2001) M. H. Lee, “Ergodic Theory, Infinite products, and Long Time Behavior in Hermitian Models,” Phys. Rev. Lett. 87, 250601 (2001).
  * Bertini  _et al._ (2021) B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, “Finite-temperature transport in one-dimensional quantum lattice models,” Rev. Mod. Phys. 93, 025003 (2021).
  * Joslin and Gray (1989) C.G. Joslin and C.G. Gray, “Calculation of transport coefficients using a modified Mori formalism,” Molecular Physics 66, 757–765 (1989).

##  Appendix A Translating between $b_{n}$ and $\mu_{2n}$

For completeness, we briefly present the relation between the Lanczos
coefficients $b_{n}$ and the moments $\mu_{2n}$ Joslin and Gray (1989).
i. From moments to Lanczos coefficients:
To calculate the Lanczos coefficients $b_{n}$ from a given set of moments
$\mu_{2n}$ we proceed as follows: we define $c_{n}=\mu_{2n}/\mu_{0}$ and
compute determinants of certain matrices constructed from the normalized
moments $c_{n}$, i.e.,

$B_{n}=\det(c_{i+j})_{0\leq i,j\leq n-1}$ (24)

where $n\geq 2$ and $B_{0}=B_{1}=1$ as well as

$C_{n}=\det(c_{i+j+1})_{0\leq i,j\leq n-1}$ (25)

where $n\geq 1$ and $C_{0}=1$. Then the Lanczos coefficients are obtained as
fractions of determinants via

$b_{2n}^{2}=\dfrac{B_{n+1}C_{n-1}}{B_{n}C_{n}}\,\,,\quad
b_{2n-1}^{2}=\dfrac{B_{n-1}C_{n}}{B_{n}C_{n-1}}\,.$ (26)

ii. From Lanczos coefficients to moments:
We take the representation $L$ of the Liouvillian $\mathcal{L}$ in the Krylov
subspace spanned by the vectors generated by the Lanczos algorithm, cf. Eq.
(3). The moments $\mu_{2n}$ can be easily read off as the upper-left element
of even powers $2n$ of the matrix $L$, i.e.,

$\mu_{2n}=(L^{2n})_{00}\,.$ (27)