# Krylov complexity of many-body localization: Operator localization in Krylov
basis

Fabian Ballar Trigueros  Perimeter Institute for Theoretical Physics,
Waterloo, Ontario, Canada N2L 2Y5  Cheng-Ju Lin  Perimeter Institute for
Theoretical Physics, Waterloo, Ontario, Canada N2L 2Y5

###### Abstract

Motivated by the recent developments of quantum many-body chaos,
characterizing how an operator grows and its complexity in the Heisenberg
picture has attracted a lot of attentions. In this work, we study the operator
growth problem in a many-body localization (MBL) system from a recently
proposed Lanczos algorithm perspective. Using the Krylov basis, the operator
growth problem can be viewed as a single particle hopping problem on a semi-
infinite chain with the hopping amplitudes given by the Lanczos coefficients.
We find that, in the MBL phase, the Lanczos coefficients scales $\sim
n/\ln(n)$ asymptotically, same as in the ergodic phase, but with an additional
even-odd alteration and effective randomness. Extrapolating the Lanczos
coefficients to the thermodynamic limit, we study the spectral function and
also find that the corresponding single-particle problem is localized for both
unextrapolated and extrapolated Lanczos coefficients, resulting in a bounded
“Krylov complexity” in time. For the MBL phenomenological model, the Lanczos
coefficients also have an even-odd alteration, but approaching to constants
asymptotically. We also find that the Krylov complexity grows linearly in time
for the MBL phenomenological model.

##  I Introduction

Motivated by the recent developments of quantum many-body chaos,
characterizing how a operator grows and its complexity under an unitary many-
body dynamics in the Heisenberg picture has been a research focus across
different subfields of physics Alexei Kitaev (2015); Maldacena  _et al._
(2016); Kitaev and Suh (2018); Roberts  _et al._ (2018); von Keyserlingk  _et
al._ (2018); Nahum  _et al._ (2018). A popular characterization of the growth
of an operator is by the out-of-time-ordered commutators (correlators), which
emphasizes the operator growth in the space-time picture Aleiner  _et al._
(2016); Chowdhury and Swingle (2017); Gärttner  _et al._ (2017); Li  _et al._
(2017); Lin and Motrunich (2018a, b); Gopalakrishnan (2018); Khemani  _et al._
(2018); Xu and Swingle (2020); Sünderhauf  _et al._ (2019); Yan  _et al._
(2020). Another characterization is via the operator entanglement, which has
also been studied in various systems Prosen and Pižorn (2007); Pižorn and
Prosen (2009); Alba  _et al._ (2019); Bertini  _et al._ (2020); MacCormack
_et al._ ; Mascot  _et al._ ; Alba (2021).

Recently, Parker et al. Parker  _et al._ (2019) proposes a different
characterization of operator growth from the perspective of recursive method
(or Lanczos algorithm). The Lanczos algorithm generates a set of so-called
Krylov basis which is used to represent the operator, and the dynamics is
therefore encoded in the time-dependent coefficients of the basis. In a sense
which we will make precise later, the operator growth problem in this basis
can be viewed as the evolution of a single particle wavefunction, hopping on a
semi-infinite chain with the hopping amplitudes given by the Lanczos
coefficients generated by the Lanczos algorithm, initialized on the first
site. One can therefore define various complexity measures from this point of
view. For example, Ref. Parker  _et al._ (2019) defines the Krylov complexity
to be the mean position the corresponding single particle. Ref. Rabinovici
_et al._ (2021) also proposes Krylov entropy, which measures the degree of
spreading of the corresponding single particle wavefunction on the semi-
infinite chain.

Ref. Parker  _et al._ (2019) hypothesizes that, for a chaotic system, the
Lanczos coefficients grow as fast as possible asymptotically. For a local
Hamiltonian, the Lanczos coefficients are upper bounded by a function linearly
in the Krylov order asymptotically (with a logarithmic correction in 1D),
which implies that the Krylov complexity grows exponentially in time (or grows
as a stretched exponential in time in 1D). In contrast, in the integrable
models, the Lanczos coefficients grow as a square root of the Krylov order or
as a constant Parker  _et al._ (2019); Noh (2021). Krylov complexity was also
studied with some general considerations Dymarsky and Gorsky (2020); Caputa
_et al._ and studied in conformal field theories Dymarsky and Smolkin , in
models with holographic correspondence Rabinovici  _et al._ (2021) and in
models with strong and almost strong edge modes Yates  _et al._ (2020); Yates
and Mitra (2021).

In this paper, we study the operator growth and the Krylov complexity in a
many-body localized (MBL) system Pal and Huse (2010); Basko  _et al._ (2006);
Oganesyan and Huse (2007); Nandkishore and Huse (2015) to contrast and
complement with the previous results in the chaotic and integrable models.
Operator growth in MBL systems have been studied using the characterization of
out-of-time-ordered correlators and found that the operators grow with a
logarithmic “operator cone” in space-time He and Lu (2017); Huang  _et al._
(2017); Chen  _et al._ (2017); Slagle  _et al._ (2017); Swingle and Chowdhury
(2017); Fan  _et al._ (2017); Deng  _et al._ (2017). Refs. MacCormack  _et
al._ ; Mascot  _et al._ studied the operator entanglement growth in a MBL
system, finding a logarithmic entanglement growth in time. Moreover, in Ref.
Cao (2021), a lower bound on the Lanczos coefficients is derived for a many-
body localized spin chain. Surprisingly, such a lower bound is also shared by
a generic chaotic spin chain. Accordingly, the operator growth problem in the
Krylov basis and its Krylov complexity certainly warrant studies.

Particularly, we study the operator growth problem in a quantum Ising model
with a longitudinal field and random transverse fields. In the strong
disordered transverse field limit, the system exhibits many-body localization
Imbrie (2016a, b). We find that, while the Lanczos coefficients in the MBL
phase has the same asymptotic behavior as in the ergodic phase, there is an
additional even-odd alternation and apparent effective randomness. By
numerically calculating the evolution of the corresponding singe particle
wavefunction hopping on a semi-infinite chain, we found that the Krylov
complexity is bounded in time in the MBL phase. Furthermore, we also find that
the corresponding single particle hopping problem is localized when
initialized on the first site, and therefore the operator is localized in the
Krylov subspace.

On the other hand, we also study the operator growth problem in the MBL
phenomenological model Huse  _et al._ (2014), where the initial operator is a
“$l$-bit flipping” operator. We find the Lanczos coefficients in this case
again show an even-odd alteration, but approaching to constants
asymptotically. We find that the Krylov complexity grows linearly in this
case, and the time evolution of the corresponding single particle wavefunction
has a traveling wavefront propagating linearly in time.

This paper is organized as follows. In Sec. II, we lay the framework used to
study the operator growth using Krylov basis and define different functions of
interests including the Krylov complexity. We then discuss the results of
operator growth in a MBL system realized in the tilted quantum Ising model in
Sec. III, where we show the numerical results and the conjectured behavior of
the Lanczos coefficients. Using the Lanczos coefficients, we then calculate
the corresponding single particle hopping problem, showing its localization
and boundedness of Krylov complexity in time. In Sec. IV, we study the
operator growth in the MBL phenomenological model, by calculating the Lanczos
coefficients and the Krylov complexity. We conclude in Sec. V with some
discussion of our results, open questions and future directions of research.

##  II Setup

We begin by setting up the framework to study the operator growth using Krylov
basis. Consider a quantum many-body system described by a Hamiltonian $H$ and
an initial local Hermitian operator ${\mathcal{O}}$. For the later
convenience, we use an operator-state correspondence description: For any
operator ${\mathcal{O}}=\sum_{i,j}O_{ij}|i\rangle\langle j|$, where
$|i\rangle$ and $|j\rangle$ are from an orthonormal basis of the Hilbert space
of consideration, we can define the corresponding operator state as
$|\mathcal{O})\equiv\sum_{i,j}O_{ij}|i\rangle|j\rangle$. Since we will
consider quantities at an infinite temperature, we define the inner product
between two operator states $|\mathcal{A})$ and $|\mathcal{B})$ as

$(\mathcal{A}\rvert\mathcal{B})\equiv\frac{\text{Tr}[\mathcal{A}^{\dagger}\mathcal{B}]}{\text{Tr}[I]}~{},$
(1)

where $I$ is the identity matrix.

The Heisenberg evolution of an operator is

${\mathcal{O}}(t)=e^{iHt}\mathcal{O}e^{-iHt}=\sum_{n=0}^{\infty}\frac{(it)^{n}}{n!}\mathcal{L}^{n}{\mathcal{O}}~{},$
(2)

where the Liouvillian superoperator is defined as
$\mathcal{L}{\mathcal{O}}\equiv[H,{\mathcal{O}}]$. Somewhat abusing the
notation, we also define
$\mathcal{L}|{\mathcal{O}})\equiv|[H,{\mathcal{O}}])$. We therefore have,
equivalently,

$|{\mathcal{O}}(t))=\sum_{n=0}^{\infty}\frac{(it)^{n}}{n!}\mathcal{L}^{n}{|\mathcal{O}})~{}.$
(3)

To motivate the Lanczos algorithm, note that we can approximate the Heisenberg
evolution Eq. (3) by summing over $n$ up to $n_{max}$. The resulting
approximated operator is therefore in the so-called Krylov space

$\mathcal{H}_{\mathcal{O}}=\text{span}\left\\{\mathcal{L}^{n}|\mathcal{O}),n=0\dots
n_{\text{max}}\right\\}~{}.$ (4)

Lanczos algorithm generates an orthonormal basis $\\{|\mathcal{O}_{n})\\}$,
$n=0\dots n_{max}$ of the Krylov space recursively which tridiagonalizes
$\mathcal{L}$. Starting with the initial operator
$\left|\mathcal{O}_{0}\right)\equiv|\mathcal{O})$, we have
$|\mathcal{O}_{1})=b_{1}^{-1}\mathcal{L}|\mathcal{O}_{0})$, where
$b_{1}^{2}=(\mathcal{L}\mathcal{O}_{0}|\mathcal{L}\mathcal{O}_{0})$. For
$n\geq 2$, we recursively define

$\displaystyle|\mathcal{A}_{n})$
$\displaystyle=\mathcal{L}|\mathcal{O}_{n-1})-b_{n-1}|\mathcal{O}_{n-2})~{},$
$\displaystyle b_{n}^{2}$
$\displaystyle\equiv(\mathcal{A}_{n}\rvert\mathcal{A}_{n})~{},$
$\displaystyle|\mathcal{O}_{n})$
$\displaystyle=\frac{1}{b_{n}}\rvert\mathcal{A}_{n})~{},$ (5)

where $b_{n}$ is the Lanczos coefficient. In this basis, $\mathcal{L}$ becomes
tridiagonal with the matrix representation

$\mathcal{L}=\begin{pmatrix}0&b_{1}&0&0&...\\\ b_{1}&0&b_{2}&0&...\\\
0&b_{2}&0&b_{3}&...\\\ 0&0&b_{3}&0&...\\\
\vdots&\vdots&\vdots&\vdots\end{pmatrix}~{}.$ (6)

By expanding the operator in the Krylov basis

$|{\mathcal{O}}(t))=\sum_{n=0}^{\infty}\varphi_{n}(t)|\mathcal{O}_{n})~{},$
(7)

the Heisenberg equation of motion becomes

$-i\partial_{t}\varphi_{n}(t)=b_{n}\varphi_{n-1}(t)+b_{n+1}\varphi_{n+1}(t)~{},$
(8)

for $n=0...n_{\text{max}}$ with the initial condition
$\varphi_{n}(0)=\delta_{n0}$, where we also define $b_{0}\equiv 0$. We
therefore see that, the equation of motion governing the coefficients in the
Krylov basis can be viewed as a single particle hopping problem on a semi-
infinite chain, with the hopping amplitudes $b_{n}$’s given by the Lanczos
coefficients. Denoting
$\vec{\varphi}(t)=(\varphi_{0}(t),\varphi_{1}(t),\dots)^{T}$, we have
$\vec{\varphi}(t)=e^{i\mathcal{L}t}\vec{\varphi}(0)$, where
$e^{i\mathcal{L}t}$ is the matrix exponential from Eq. (6).

This perspective indeed motivates a natural consideration of complexity. One
intuition of the complexity of $|\mathcal{O}(t))$ comes from the complexity of
$|\mathcal{O}_{n})$. In particular, $|\mathcal{O}_{n})$ involves an $n$-nested
commutator with the Hamiltonian, and the operator is therefore more complex
and more nonlocal with the increased order $n$. The order $n$ therefore can be
served as a measure of the operator complexity, which motivates the definition
of Krylov complexity in Ref. Parker  _et al._ (2019),

$C_{K}(t)\equiv\sum_{n=0}^{\infty}n|(\mathcal{O}(t)|\mathcal{O}_{n})|^{2}~{}.$
(9)

It can also be interpreted as the mean position of the particle in the
corresponding single particle hopping problem. On the other hand,
$\vec{\varphi}(t)$ also gives us a measure of how much resource one has to use
to have a good approximation of $|\mathcal{O}(t))$. If $\vec{\varphi}(t)$ is
concentrated or localized at small orders of $n$, then one does not need too
high of a truncation order $n_{\text{max}}$ to obtain a good approximation of
$|\mathcal{O}(t))$ in the Lanczos algorithm. As we will see, this is indeed
the case for the operators in the MBL systems.

It is worth to mention that, to reconstruct the full information of the
operator $|\mathcal{O}(t))$, one would need both $\vec{\varphi}(t)$ (or
equivalently $b_{n}$) and $|\mathcal{O}_{n})$. However, there are some
physical quantities which can be obtained from $\vec{\varphi}(t)$ or $b_{n}$
solely. The notable ones are the auto-correlation function

$F(t)\equiv\frac{\text{Tr}[\mathcal{O}(t)^{\dagger}\mathcal{O}(0)]}{\text{Tr}[I]}=(\mathcal{O}_{0}|\mathcal{O}(t))=\varphi_{0}(t)~{}.$
(10)

The spectral function, defined by its Fourier transformation
$\Phi(\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(t)e^{-i\omega t}dt$, can
then be expressed as

$\Phi(\omega)=\sum_{m}|(\mathcal{O}_{0}|E_{m})|^{2}\delta(\omega-E_{m})~{},$
(11)

where $|E_{m})$ is an eigenstate of $\mathcal{L}$ in Eq. (6) with the
eigenvalue $E_{m}$. We also note that the analytical continuation of the
spectral function can also be formally expressed through the continuous
fraction of $b_{n}$, which is obtained by solving Eq. (8) with Laplace
transformation Viswanath and Müller (1994); Parker  _et al._ (2019).

##  III Operator growth in a MBL system

In this section, we study the operator growth problem from the perspective of
Lanczos algorithm in a MBL system. We consider a 1d spin chain with $L$ sites
with an open boundary condition

$H=-J\sum_{\langle ij\rangle}Z_{i}Z_{j}-g\sum_{j}X_{j}+\sum_{j}h_{j}Z_{j}~{},$
(12)

where $\langle ij\rangle$ denotes the nearest neighbors, $X_{j}$, $Y_{j}$,
$Z_{j}$ are the Pauli matrices on site $j$, $J=1$ is the energy unit,
$g=-1.05$ and $h_{i}$ is the random transverse field drawing from a uniform
distribution $[-h,h]$. (For an odd chain, we consider sites
$j=-(L\\!-\\!1)/2\dots(L\\!-\\!1)/2$; for an even chain, we consider sites
$j=-L/2\\!+\\!1\dots L/2$.) Such a spin chain exhibits MBL when the disorder
strength $h\gtrsim 3$ as shown in Appendix A. For all of our following
calculations, we generate $10^{3}$ disorder realizations, and all the
disorder-averaged quantities are denoted with an overline.

###  III.1 Lanczos coefficients

Figure 1: Disorder-averaged Lanczos coefficients $\overline{b_{n}}$ for
several disorder strengths $h$ in the tilted quantum Ising model Eq. (12),
starting with the $Z_{0}$ initial operator and system size $L=13$. In the MBL
phase $h\gtrsim 3$, $\overline{b_{n}}$ has the same asymptotic behavior as in
the ergodic phase, but with an additional even-odd alteration.  Figure 2:
Disorder-averaged Lanczos coefficients $\overline{b_{n}}$ for several system
sizes $L$ and $h=7$, with the initial operator $Z_{0}$. The onset point of the
plateau of $\overline{b_{n}}$ increases with the system size $L$. We therefore
conclude that the saturation of $\overline{b_{n}}$ is due to the finite size.
Inset: the disorder averaged Lanczos coefficients starting with the the
$X_{0}$, $Y_{0}$ and $Z_{0}$ operators. The observed even-odd alteration of
the Lanczos coefficients is independent of the choice of the initial
operators.

First we generate the Lanczos coefficients and examine their behaviors for
various disorder strengths $h$. In Fig. 1, we show the disorder averaged
Lanczos coefficients in a system size $L=13$ with the initial operator
$\mathcal{O}_{0}\\!=\\!Z_{j=0}$. We note that before the saturation or the
plateau at large $n$, $\overline{b_{n}}$ appears to have an asymptotic
behavior $\sim n/\ln(n)$, consistent with Refs. Parker  _et al._ (2019); Cao
(2021). In the linear-linear plot (Fig. 2), the traces indeed appear to be
linear and the logarithmic correction is difficult to notice. However, as we
see in the log-log plot in Fig. 1, the logarithmic correction has the effect
of skewing the apparent exponent $b_{n}\sim n^{\delta}$ of $\delta$ from
$\delta=1$ to $\delta\approx 0.8$. We also plot our data with $\sim n/\ln(n)$
to show that the logarithmic correction is indeed present, and such a behavior
is expected in both the MBL and ergodic phases. On the other hand, we observe
that in the MBL phase $h\gtrsim 3$, the averaged Lanczos coefficients $b_{n}$
starts to show an even-odd alteration, with the amplitude of the alteration
increases with an increased disorder strength $h$.

In Fig. 2, we plotted $\overline{b_{n}}$ generated from different system
sizes. We see that the starting $n$ of the saturation or plateau increases as
the system size $L$ increases. We therefore conclude that the plateau behavior
of $\overline{b_{n}}$ is due to the finite system size $L$. In fact, this
appears to be a common finite-size feature of the Lanczos coefficients in
various models Viswanath and Müller (1994); Yates  _et al._ (2020).

To further confirm that the data before the plateau is representative of the
thermodynamic limit, we push the calculation up to $L=31$ using a different
algorithm which uses the product of Pauli matrices as the basis and only
tracks the coefficients. This enables a large $L$ computation, but can only
allow calculation to reach smaller $n$ since it requires exponentially more
resources to represent $|\mathcal{O}_{n})$ with increasing $n$. We note that
the Lanczos coefficients for $L=31$ show the same behavior for small $n$ as
the smaller system sizes. This again supports our conclusion that the
asymptotic scaling and the even-odd alteration of Lanczos coefficients before
the saturation is representative of the thermodynamic limit. In the inset of
Fig. 2, we also show that the aforementioned observed behaviors of $b_{n}$ is
independent of the initial local operator.

With the above observations, the mathematical results in Ref. Cao (2021)
(despite the different choice of the disorder distribution), and looking at
several disorder realizations of $b_{n}$ such as the ones shown in Fig. 3, we
conjecture the Lanczos coefficients in the MBL system will have the following
asymptotic behavior:

$b_{n}=\begin{cases}a_{e}\ \frac{n}{\ln n}+c_{e}+\Gamma_{n},&\text{if }n\
\textrm{even}\\\ a_{o}\ \frac{n}{\ln n}+c_{o}+\Gamma_{n},&\text{if }n\
\textrm{odd}~{}.\end{cases}$ (13)

More precisely, by the above formula, we expect the parameters $a_{e(o)}$ and
$c_{e(o)}$ to be different for each disorder realization, drawing from some
probability distributions. On the other hand, we also expect an additional
effective randomness $\Gamma_{n}$ for each $n$, drawing from a probability
distribution independent of $n$.

Figure 3: (a)-(d) Disorder realizations of the Lanczos coefficients for Eq.
(12), with the initial operator $Z_{0}$, $h=7$ and $L=13$. We use the data
points in the range of $n=10$ to $n=40$ to fit the conjectured asymptotic
behavior of the Lanczos coefficients Eq. (13), and then use the fitted results
to extrapolate the Lanczos coefficients for $n>40$ for each realization.

To further verify our conjectured formula, for each disorder realization, we
take the Lanczos coefficients in the range $n=10$ to $n=40$ and fit them with
$y_{\text{fit}}=a_{e(o)}n/\ln n+c_{e(o)}$ using least-square fit for the even
and odd branch respectively. The resulting fits are shown in Fig. 3 for some
disorder realizations and the histogram of the fitted parameters are shown in
Fig. 4 (a)-(d). After fitting the data, we then examine the difference between
the data and the fit $y-y_{\text{fit}}$ and plot the histogram of the
differences collected from each $n$ in $n=10$ to $n=40$ from both even and odd
branch, as shown in Fig. 4(e). We see that the histograms in Fig. 4 all look
very close to Gaussian distributions, which supports our conjecture that these
parameters can be viewed as drawing from some probability distributions
independently, and the means and variances of the probability distributions
can be estimated from the Gaussian distributions shown in Fig. 4, according to
the central limit theorem.

Figure 4: (a)(b) The histograms of the fitted parameters $a_{o}$ and $c_{o}$
from Eq. (13) on the odd branch. (c)(d) The histograms of the fitted data
$a_{e}$ and $c_{e}$ from the even branch. (e) The histogram of $\Gamma_{n}$ in
Eq. (13), obtained by taking the difference between the data and the fitted
results $y-y_{\text{fit}}$, where $y_{\text{fit}}=a_{o(e)}n/\ln(n)+c_{o(e)}$
with the fitted parameter for each realization. The Gaussian distributions of
the histogram supports the conjectured effective randomness of the parameters.

An advantage of using recursive method or the Lanczos algorithm is that it
allows us to extrapolate the Lanczos coefficients in an attempt to assert some
quantities such as the spectral function in the thermodynamic limit with
various extrapolation schemes Viswanath and Müller (1994); Khait  _et al._
(2016); Parker  _et al._ (2019). With the conjectured asymptotic behavior Eq.
(13), we extrapolate the Lanczos coefficients to large $n$ and to the
thermodynamic limit using the following scheme. For each disorder realization
of $b_{n}$, we obtained the fitted parameters of $a_{e(o)}$ and $c_{e(o)}$
from fitting the data in the range $n=10$ to $n=40$ as mentioned previously.
We then take these $b_{n}$ from the numerical results up to $n=40$, and then
extrapolate $b_{n}$ for $n>40$ to $n_{\max}=1000$ using Eq. (13) with the
fitted $a_{e(o)}$ and $c_{e(o)}$ for this realization and the randomness
$\Gamma_{n}$ drawing from a Gaussian ensemble with the mean and the variance
given by the parameters in Fig. 4(e). In the following, we will study various
quantities using the unextrapolated $b_{n}$ (usually denoted by solid lines in
the figures) and using the extrapolated $b_{n}$ (usually denoted by dashed
lines in the figures).

###  III.2 Spectral Function

First let us examine the behavior of the spectral function using Eq. (11). In
terms of the corresponding single particle hopping problem, the spectral
function is just the density of states weighted by the weight function which
is the probability of the eigenfunction on the first site
$|(\mathcal{O}_{0}|E_{j})|^{2}$ and therefore can be obtained by diagonalizing
Eq. (6). Our attempt is to extract the analytical behavior at low and high
frequencies in the thermodynamic limit. Since numerically, one can only obtain
discrete eigenvalues from diagonalizing Eq. (6), the spectral function will be
a sum of delta functions with some weights. It is therefore more customary to
consider the cumulative spectral function

$\chi(\omega)\equiv\int_{\Omega=-\infty}^{\omega}\Phi(\Omega)d\Omega=\sum_{E_{j}<\omega}|(\mathcal{O}_{0}|E_{j})|^{2}~{},$
(14)

and extract the smooth part of the cumulative spectral function. This also
gives us a convenient way to average over disorder realizations, where one
would just sum over all the weights below $\omega$ in Eq. (14) from all the
disorder realizations and then divide it by the number of realizations.

We examine the high frequency behavior of $\overline{\chi}(\omega)$ first,
using the unextrapolated and extrapolated $b_{n}$. In Ref. Parker  _et al._
(2019), it was proposed that if the Lanczos coefficients grows linearly in
$n$, then the spectral function at the high frequencies decays exponentially
$\Phi(\omega)\sim\exp(-\frac{\pi|\omega|}{2\alpha})$, where $\alpha$ is the
coefficient of the linear growth $b_{n}\sim\alpha n$. Since
$\overline{\chi}(\omega\\!\rightarrow\\!\infty)=1$, we plot
$1-\overline{\chi}(\omega\\!\rightarrow\\!\infty)$ as shown in Fig. 5. From
the figure, we see that at high frequency, $\overline{\chi}(\omega)$ decays
exponentially, and therefore so does $\overline{\Phi}(\omega)$ for both the
unextrapolated and extrapolated $b_{n}$. Furthermore, we extract the
exponential decay as $\chi(\omega)\sim\exp(-\frac{\pi\omega}{2\alpha})$, where
$\alpha\approx 4.36$, which is close to the mean values of $a_{e(o)}$ as shown
in Fig. 4. We see that the high frequency behavior of the spectral function is
indeed blind to the MBL or ergodic phase Cao (2021).

On the other hand, at the zero frequency, we observe the presence of a delta
function in the spectral function $\Phi(\omega)\sim A\delta(\omega)+\dots$
(manifested as a step function for the cumulative spectral function
$\chi(\omega)$) for the unextrapolated data $n_{\text{max}}=40$ and also for
the extrapolated data. The existence of the delta function is indeed a
hallmark of MBL, since it implies that the auto-correlation function $F(t)$
will decay to $A\neq 0$ at infinite time, assuming the rest part of the
spectral function is analytical. These results also suggest that our
extrapolation scheme is reasonable, since it does not change the expected
qualitative behaviors of the spectral function in MBL.

It is also interesting to see if the spectral function shows any singular
behavior at low frequencies in addition to the zero-frequency delta function.
The low-frequency singular behavior of the spectral function would manifest in
the long-time behavior of the auto-correlation function. While not shown in
the figure, we attempt to fit the low frequency part of the cumulative
spectral function using the $n_{\text{max}}=1000$ result with
$\omega^{\gamma+1}$ and found that $\gamma\approx 0.39$ (so that
$\Phi(\omega)\sim\omega^{0.39}$). However, we note that it is still fairly
challenging to extract the low-frequency singular behavior from our present
numerical data. It would be ideal if one can obtain an analytical low-
frequency result from the physics of MBL or from our conjectured asymptotic
formula Eq. (13) and to compare with our disorder-averaged spectral function
Nandkishore  _et al._ (2014); Johri  _et al._ (2015); Pekker  _et al._ (2017);
Varma  _et al._ (2019); Gopalakrishnan and Parameswaran (2020); Herbrych  _et
al._ . On the other hand, if one has an expected low-frequency behavior of the
spectral function, one can also use this information to modify the conjectured
behavior of $b_{n}$ Viswanath and Müller (1994); Khait  _et al._ (2016). In
any case, the low-frequency behavior of the spectral function in MBL warrants
a further investigation in the future using the recursive method Khait  _et
al._ (2016).

Figure 5: The disorder-averaged cumulative spectral function
$\overline{\chi(\omega)}$. The high frequency behavior appears to decay
exponentially, irregardless of calculated from the unextrapolated $b_{n}$ with
$n_{\max}=40$ or from the extrapolated $b_{n}$ with $n_{\max}>40$. Note that
the cumulative spectral function appears to be converged when $n_{\max}=200$.
Inset: the overall disorder-averaged cumulative spectral function
$\overline{\chi(\omega)}$. Noe that the existence of $\delta(\omega)$ in the
spectral function, a hallmark of many-body localization, is unaffected by our
extrapolation scheme of the Lanczos coefficients.

###  III.3 Krylov complexity

As discussed in the Sec. II, we can view the operator growth problem as a
single particle hopping problem on a semi-infinite chain. The time evolution
of the single particle wavefunction $\vec{\varphi}(t)$ therefore reflects the
growth of the operator in the Krylov basis
$|\mathcal{O}(t))=\sum_{n}\varphi_{n}(t)|\mathcal{O}_{n})$.

First, we examine the mean position of the corresponding single particle
wavefunction, namely the Krylov complexity, from both unextrapolated and
extrapolated $b_{n}$. In Fig. 6, we show the disorder-averaged Krylov
complexity in the MBL phase ($h=7$) in time. The long-time averaged Krylov
complexity appears to be bounded in time and saturates quickly to a value much
smaller than the truncation $n_{\text{max}}$ from both extrapolated and
unextrapolated $b_{n}$. Note that $\overline{C_{K}}(t)$ shows a fluctuation
around the long-time averaged value at long time, and the magnitude of the
fluctuation appears to grow logarithmically in time. However, by examining
$\overline{C_{K}}(t)$ for several $n_{\text{max}}$, we notice the onset time
of such a fluctuation increases with $n_{\text{max}}$, suggesting it as a
result of finite $n_{\text{max}}$. We therefore conclude that the Krylov
complexity is bounded in time in MBL systems.

Figure 6: Disorder-averaged Krylov complexity $C_{K}(t)$ in MBL ($h=7$) from
the unextrapolated Lanczos coefficents (solid lines) and from the extrapolated
Lanczos coefficents (dashed lines). The Krylov complexity $C_{K}(t)$ appears
to be bounded in time and we expect the fluctuations in time at the long times
are due to finite $n_{max}$. Note that in any case, the truncation $n_{max}$
is much greater than $C_{K}(t)$.

Figure 7: Disorder-averaged probability distribution
$\overline{|\varphi_{n}(t)|^{2}}$ in the MBL phase $h=7$ at different times
from (a) unextrapolated $b_{n}$ with $n_{max}=40$ and (b) extrapolated $b_{n}$
with $n_{max}=1000$. The wavefunctions appear to be localized at small $n$
even at large times.

The growth of the operator can be more easily visualized by examining the
probability distribution of the single particle wavefunction at various times.
In Fig. 7, we show $\overline{|\varphi_{n}(t)|^{2}}$ from both unextrapolated
and extrapolated Lanczos coefficients. Remarkably, we see the wavefunctions
are localized around the first site even for a exponentially long time. We
therefore conclude that the corresponding single particle hopping problem in
the MBL phase is localized if initialized on the first site . This is in stark
contrast with the case in the ergodic phase, where the single particle
wavefunction propagates to large $n$ superpolynomially fast Parker  _et al._
(2019).

Our result also implies that the Krylov method can be an efficient method to
simulate the operator dynamics in MBL. In practice, we have to choose a
truncation $n_{\text{max}}$ as the cut-off order of the Lanczos algorithm. If
$\varphi_{n}(t)$ is localized at small $n$, then the truncation error
$\epsilon(t)=\sum_{n=n_{\text{max}}+1}^{\infty}|\varphi_{n}(t)|^{2}$ will be
small. In particular, from our results of $|\varphi_{n}(t)|^{2}$, we expect
this error will be small for MBL systems even at exponentially long times.
Note that while $|\varphi_{n}(t)|^{2}$ seems to be localized and not changing
too much with times, the phases are still evolving. This can still cause the
operator entanglement entropy or out-of-time-ordered correlator to change at
long times in MBL systems.

##  IV Operator growth in the MBL phenomenological model

In this section, we study the operator growth and the Krylov complexity in the
MBL phenomenological model Huse  _et al._ (2014) to compare with the results
in the previous section. A defining feature of the MBL system is the existence
of the local integrals of motion Serbyn  _et al._ (2013); Chandran  _et al._
(2015), namely the “$\ell$-bits”, resulting in the following phenomenological
model

$\displaystyle H_{\text{ph}}$
$\displaystyle=\sum_{j}J_{j}^{(0)}\hat{\tau}^{z}_{j}+\sum_{i<j}J^{(1)}_{ij}\hat{\tau}^{z}_{i}\hat{\tau}^{z}_{j}$
$\displaystyle+\sum_{n=2}^{j-i}\sum_{i<j,\\{k\\}}J^{(n)}_{i\\{k\\}j}\hat{\tau}^{z}_{i}\hat{\tau}^{z}_{k_{1}}\dots\hat{\tau}^{z}_{k_{n-1}}\hat{\tau}^{z}_{j}~{},$
(15)

where $\\{k\\}$ denotes the set of sites $i<k_{1},\dots,k_{n-1}<j$,
$\hat{\tau}^{z}_{j}$ are the $\ell$-bits which can be related to the physical
bits $Z_{j}$ by the unitary which diagonalizes the Hamiltonian, and
$J^{(2)}_{ij}$ and $J^{(n)}_{i\\{k\\}j}$ are the interactions among the
$\ell$-bits which decays exponentially in the number of $\ell$-bits involved
and the distances between the $\ell$-bit. To simplify the problem, we assume
$J^{(n)}_{i\\{k\\}j}=K^{(n)}_{i\\{k\\}j}\exp(-nr/\xi)$ (including $n=0$ and
$1$), where the range $r=j-i$, $\xi$ being the localization length, and the
parameters $K^{(n)}_{i\\{k\\}j}$ are drawing from an uniform distribution
$[-W,W]$ Swingle and Chowdhury (2017); MacCormack  _et al._ . In particular,
we consider the parameters $L=13$ as the system size, $\xi=0.6$ and $W=5.0$,
with $10^{3}$ disorder realizations.

###  IV.1 Lanczos coefficients

Here we generate the Lanczos coefficients from the Hamiltonian Eq. (IV) with
the initial operator $\mathcal{O}_{0}=\hat{\tau}^{x}_{0}$. To calculate the
Lanczos coefficients, it is easy to split the Hamiltonian to a part that
involves $\tau_{0}^{z}$ and the rest:
$H_{\text{ph}}=J_{\text{eff}}\tau^{z}_{0}+H_{\text{other}}$, where

$\displaystyle J_{\text{eff}}=J^{(0)}_{0}+\sum_{j\neq
0}J^{(1)}_{0j}\tau^{z}_{j}+\sum_{n=2}^{L-1}\sum_{\\{j\\}\neq
0}J^{(n)}_{0\\{j\\}}\tau^{z}_{j_{1}}\dots\tau^{z}_{j_{n}}~{}.$ (16)

It is easy to see that $J_{\text{eff}}\tau^{z}$ is indeed the only part that
generates the Heisenberg evolution of $\tau^{x}_{0}$. In fact, one can obtain
$\mathcal{L}^{2n-1}\tau_{0}^{x}=i\tau_{0}^{y}(J_{\text{eff}})^{2n-1}$ and
$\mathcal{L}^{2n}\tau_{0}^{x}=\tau_{0}^{x}(J_{\text{eff}})^{2n}$. By
expressing $J_{\text{eff}}$ as a diagonal matrix with a dimension
$2^{L-1}\\!\times\\!2^{L-1}$ where the entries corresponding to different
$\tau^{z}$ configurations, we obtain the following recursive method to
generate the Lanczos coefficents: Define $B_{0}=I$ as a $2^{L-1}\times
2^{L-1}$ identity matrix, $B_{1}=J_{\text{eff}}$, $b_{0}=1$ and
$b_{1}^{2}=\frac{1}{2^{L}}\text{Tr}[B_{1}^{2}]$, then

$\displaystyle B_{n+1}$
$\displaystyle=\frac{1}{b_{n}}J_{\text{eff}}B_{n}-\frac{b_{n}}{b_{n-1}}B_{n-1}~{},$
$\displaystyle b_{n+1}^{2}$
$\displaystyle=\frac{1}{2^{L}}\text{Tr}[B_{n+1}^{2}]~{},$ (17)

for $n>1$. This gives us $O_{2n-1}=b_{2n-1}^{-1}i\tau_{0}^{y}B_{2n-1}$ and
$O_{2n}=b_{2n}^{-1}\tau_{0}^{x}B_{2n}$ from the Lanczos algorithm.

Figure 8: (a) The disorder-averaged Lanczos coefficients $\overline{b_{n}}$
from the MBL phenomenological model Eq. (IV) with the initial operator
$\tau^{x}_{0}$. The Lanczos coefficients have an even-odd alteration, both
approaching to a constant in large $n$. Inset: The disorder averaged Lanczos
coefficients $\overline{b_{n}}$ at small $n$. (b) The Lanczos coefficients
$b_{n}$ from the MBL phenomenological model for several disorder realizations.

In Fig. 8(a), we show the disorder averaged Lanczos coefficients
$\overline{b_{n}}$ from the phenomenological model. Note that
$\overline{b_{n}}$ appears to approach a constant, but also with an even-odd
alteration and the odd branch usually has a higher value than the even branch.
We also show several disorder realizations of $b_{n}$ in Fig. 8(b).

###  IV.2 Krylov Complexity

For each disorder realization of $b_{n}$, we can solve for the corresponding
single particle dynamics, representing its operator growth. In Fig. 9, we
plotted the disorder averaged Krylov complexity $\overline{C_{K}(t)}$ of MBL
phenomenological model. We notice that, after $t\gtrsim 1$, the Krylov
complexity appears to grow linearly in time and saturate. However, by plotting
$\overline{C_{K}(t)}$ with several different $n_{max}$, we conclude that the
saturation is due to the finite-$n_{max}$, and therefore we expect the
linearly growing behavior of $\overline{C_{K}(t)}$ will continue in the true
$n_{max}\rightarrow\infty$ limit.

Figure 9: The disorder averaged Krylov complexity $\overline{C_{K}(t)}$ in the
MBL phenomenological model. After a short time, $\overline{C_{K}(t)}$ appears
to grow linearly in time before it reaches the saturation. By comparing the
results with several $n_{\text{max}}$, we conclude that the saturation is a
finite-$n_{\text{max}}$ effect.

The operator growth in the Krylov basis can be even better visualized by
looking at $\overline{|\varphi_{n}(t)|^{2}}$ at different times, which we
plotted in Fig. 10. As we can see, the wavefunction has a feature of a
propagating “wavefront” towards large $n$, which appears to propagate linearly
in time. The time scale where the “wavefront” hits $n_{max}$ is also the time
scale $\overline{C_{K}(t)}$ starts to saturate. It is therefore clear that the
saturation of $\overline{C_{K}(t)}$ is indeed a finite-$n_{max}$ effect, and
the single particle wavefunction $\vec{\varphi}(t)$ will propagate
indefinitely in the $n_{max}\\!\rightarrow\\!\infty$ limit. This is in stark
contrast with the growth of the “p-bit” shown in Sec. III, where the operator
appears to be localized.

Figure 10: The disorder averaged probability distribution
$\overline{|\varphi_{n}(t)|^{2}}$ at different times. The wavefunction has a
propagating wavefront, where the wavefront appears to propagate towards higher
$n$ linearly in time.

##  V Conclusions

In this paper, we study the operator growth problem in a MBL system using
Krylov basis. From this point of view, the operator growth problem can be
mapped to a single-particle hopping problem on a semi-infinite chain with the
hopping amplitudes given by the Lanczos coefficients. In particular, we
consider the problem with the initial operators as a “p-bit” ($X_{0}$) and as
an “$\ell$-bit” ($\tau^{x}_{0}$). We find that, in the “p-bit” case, the
asymptotic behavior of the Lanczos coefficients behaves as $b_{n}\sim
n/\log(n)$, same as in the ergodic case. However, the presence of the even-odd
alteration and relatively strong effective randomness are the features that
are distinguishable from the ergodic phase. These features could potentially
be used as an alternative diagnostic for MBL. It is worth to mention that, we
also examine the Lanczos coefficients in the MBL system realized by a
Heisenberg chain (shown elsewhere Trigueros (2021)). We find that the
aforementioned features are still present, and therefore we conclude that
these features are generic for MBL and independent of the choices of the
model. On the other hand, for the “$\ell$-bit”, the Lanczos coefficients also
show an even-odd alteration, but approach to constants at large $n$, in
contrast with the results for the “p-bit”.

With our conjectured asymptotic behavior of the Lanczos coefficients in MBL
for the “p-bit”, we extrapolate the Lanczos coefficients to the thermodynamic
limit, and study the spectral function from the resulting extrapolated Lanczos
coefficients. We found that the existence of the delta function at zero
frequency, which is a hallmark of MBL phase, is unaffected by our
extrapolation. Moreover, we also show that the spectral function decays
exponentially at high frequencies. The low-frequency behavior of the spectral
function is more challenging to extract. It is therefore interesting to see if
the conjectured asymptotic Lanczos coefficients can provide us some analytical
understanding of the low-frequency asymptotic behavior of the spectral
function.

Perhaps somewhat surprisingly, even though the asymptotic behaviour of the
Lanczos coefficients are the same in both MBL and ergodic phases, the
corresponding single particle hopping problem is localized in the MBL phase if
initialized on the first site, resulting in a bounded Krylov complexity in
time. Note that we show this for both extrapolated and unextrpolated Lanczos
coefficients. In contrast, the Krylov complexity in the ergodic phase grows
superpolynomially in time Parker  _et al._ (2019), corresponding to the
delocalization of the single particle. On the other hand, the Krylov
complexity for “$\ell$-bit” grows linearly in time, due to the different
choices of the Krylov basis.

We comment on some questions and future directions motivated from our results.
As our results have shown, the single particle hopping problem is localized
when initialized on the first site in MBL phase but propagates
superpolynomially fast and therefore is delocalized in the ergodic phase. From
this perspective, can we study or understand the MBL-ergodic transition from
this corresponding single particle localization-delocalization transition? It
would also be interesting to understand the properties of the spectral
function or auto-correlation function at the MBL transition from this point of
view, and to deepen our understanding of possible experimental signatures of
MBL-ergodic transition.

Since in the MBL phase, the operator is exponentially localized in the Krylov
space even at long times, Lanczos algorithm is a good method to obtain the
operator dynamics. In other words, one does not need too high of a truncation
$n_{\max}$ to obtain a good approximation of $|\mathcal{O}(t))$. However, we
note that even in the MBL case, it can still require exponential resources to
represent $|\mathcal{O}_{n})$. Can we design a Krylov-based numerical method
to simulate the MBL dynamics more efficiently? Some Krylov-based hybrid
numerical methods might be possibilities to alleviate the obstacle, enabling
us to simulate MBL dynamics efficiently in larger system sizes and longer
times on a classical computer.

Finally, our result suggests a correspondence of the “dynamical phases” and
“computational complexity”. That is to say, the complexity of using Krylov
method to calculate the operator dynamics is low for systems in MBL phase,
independent of the microscopic model. It is therefore interesting to see how
much of this correspondence holds for different types of the dynamical phase.
Analogous to the recent developments of using computational complexity to
understand and characterize different phases of matters, our results could
suggest a new fruitful direction of using complexity to understand and
characterize dynamical phases.

###### Acknowledgements.

We thank Timothy Hsieh and Liujun Zou for the valuable discussions. We
acknowledge supports from Perimeter Institute for Theoretical Physics. Fabian
Ballar Trigueros would like to thank the PSI program for facilitating this
research. Research at Perimeter Institute is supported in part by the
Government of Canada through the Department of Innovation, Science and
Economic Development Canada and by the Province of Ontario through the
Ministry of Colleges and Universities.

##  Appendix A Gap ratio statistics of the random tilted Ising model

Figure 11: The gap ratio statistics $r$ for various disorder strength $h$ and
system size $L$. Note that the dashed lines corresponding to the mean values
from a Gaussian orthogonal ensemble ($\overline{r}\approx 0.53$) and a Poisson
ensemble ($\overline{r}\approx 0.39$). We therefore see that the system is in
MBL when $h\gtrsim 3$.

In this appendix, we show the gap ratio statistics Oganesyan and Huse (2007)
for the tilted quantum Ising model Eq. (12). In particular, we define
$r=\min(\delta_{n},\delta_{n-1})/\max(\delta_{n},\delta_{n-1})$, where
$\delta_{n}=E_{n+1}-E_{n}$ is the energy difference between two consecutive
eigenstates. We take the middle half of the spectrum to calculate the mean
value of $r$, and then average over $10^{3}$ disorder realizations to obtain
$\overline{r}$ for various disorder strength $h$ and system size $L$, plotted
in Fig. (11). Note that for the Gaussian orthogonal ensemble, one expects
$\overline{r}\approx 0.53$ while $\overline{r}\approx 0.39$ for the Poisson
distribution. We therefore see that, for the Hamiltonian Eq. (12), $h\gtrsim
3$ is the regime where the system is in the many-body localization phase.

## References

  * Alexei Kitaev (2015) Alexei Kitaev, “A Simple Model of Quantum Holography,”  (2015).
  * Maldacena  _et al._ (2016) Juan Maldacena, Stephen H. Shenker, and Douglas Stanford, “A Bound on Chaos,” J. High Energ. Phys. 2016, 106 (2016).
  * Kitaev and Suh (2018) Alexei Kitaev and S. Josephine Suh, “The Soft Mode in the Sachdev-Ye-Kitaev Model and Its Gravity Dual,” J. High Energ. Phys. 2018, 183 (2018).
  * Roberts  _et al._ (2018) Daniel A. Roberts, Douglas Stanford, and Alexandre Streicher, “Operator growth in the SYK model,” J. High Energy Phys. 2018, 122 (2018).
  * von Keyserlingk  _et al._ (2018) C. W. von Keyserlingk, Tibor Rakovszky, Frank 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) Adam Nahum, Sagar Vijay, and Jeongwan Haah, “Operator spreading in random unitary circuits,” Phys. Rev. X 8, 021014 (2018).
  * Aleiner  _et al._ (2016) Igor L. Aleiner, Lara Faoro, and Lev B. Ioffe, “Microscopic Model of Quantum Butterfly Effect: Out-of-Time-Order Correlators and Traveling Combustion Waves,” Ann. Phys. 375, 378–406 (2016).
  * Chowdhury and Swingle (2017) Debanjan Chowdhury and Brian Swingle, “Onset of many-body chaos in the O(N) model,” Phys. Rev. D 96, 065005 (2017).
  * Gärttner  _et al._ (2017) Martin Gärttner, Justin G. Bohnet, Arghavan Safavi-Naini, Michael L. Wall, John J. Bollinger, and Ana Maria Rey, “Measuring Out-of-Time-Order Correlations and Multiple Quantum Spectra in a Trapped-Ion Quantum Magnet,” Nat. Phys. 13, 781–786 (2017).
  * Li  _et al._ (2017) Jun Li, Ruihua Fan, Hengyan Wang, Bingtian Ye, Bei Zeng, Hui Zhai, Xinhua Peng, and Jiangfeng Du, “Measuring Out-of-Time-Order Correlators on a Nuclear Magnetic Resonance Quantum Simulator,” Phys. Rev. X 7, 031011 (2017).
  * Lin and Motrunich (2018a) Cheng-Ju Lin and Olexei I. Motrunich, “Out-of-time-ordered correlators in a quantum Ising chain,” Phys. Rev. B 97, 144304 (2018a).
  * Lin and Motrunich (2018b) Cheng-Ju Lin and Olexei I. Motrunich, “Out-of-time-ordered correlators in short-range and long-range hard-core boson models and in the Luttinger-liquid model,” Phys. Rev. B 98, 134305 (2018b).
  * Gopalakrishnan (2018) Sarang Gopalakrishnan, “Operator Growth and Eigenstate Entanglement in an Interacting Integrable Floquet System,” Phys. Rev. B 98, 060302 (2018).
  * Khemani  _et al._ (2018) Vedika Khemani, David A. Huse, and Adam Nahum, “Velocity-dependent Lyapunov exponents in many-body quantum, semiclassical, and classical chaos,” Phys. Rev. B 98, 144304 (2018).
  * Xu and Swingle (2020) Shenglong Xu and Brian Swingle, “Accessing scrambling using matrix product operators,” Nature Physics 16, 199–204 (2020).
  * Sünderhauf  _et al._ (2019) Christoph Sünderhauf, Lorenzo Piroli, Xiao-Liang Qi, Norbert Schuch, and J. Ignacio Cirac, “Quantum chaos in the Brownian SYK model with large finite N : OTOCs and tripartite information,” J. High Energy Phys. 2019, 38 (2019).
  * Yan  _et al._ (2020) Bin Yan, Lukasz Cincio, and Wojciech H. Zurek, “Information scrambling and loschmidt echo,” Phys. Rev. Lett. 124, 160603 (2020).
  * Prosen and Pižorn (2007) Tomaž Prosen and Iztok Pižorn, “Operator space entanglement entropy in a transverse Ising chain,” Phys. Rev. A 76, 032316 (2007).
  * Pižorn and Prosen (2009) Iztok Pižorn and Tomaž Prosen, “Operator space entanglement entropy in XY spin chains,” Phys. Rev. B 79, 184416 (2009).
  * Alba  _et al._ (2019) V. Alba, J. Dubail, and M. Medenjak, “Operator entanglement in interacting integrable quantum systems: The case of the rule 54 chain,” Phys. Rev. Lett. 122, 250603 (2019).
  * Bertini  _et al._ (2020) Bruno Bertini, Pavel Kos, and Tomaz Prosen, “Operator entanglement in local quantum circuits I: Chaotic dual-unitary circuits,” SciPost Phys. 8, 67 (2020).
  * (22) Ian MacCormack, Mao Tian Tan, Jonah Kudler-Flam, and Shinsei Ryu, “Operator and entanglement growth in non-thermalizing systems: Many-body localization and the random singlet phase,”  arXiv:2001.08222 .
  * (23) Eric Mascot, Masahiro Nozaki, and Masaki Tezuka, “Local Operator Entanglement in Spin Chains,”  arXiv:2012.14609 .
  * Alba (2021) Vincenzo Alba, “Diffusion and operator entanglement spreading,” Phys. Rev. B 104, 094410 (2021).
  * Parker  _et al._ (2019) Daniel E. Parker, Xiangyu Cao, Alexander Avdoshkin, Thomas Scaffidi, and Ehud Altman, “A Universal Operator Growth Hypothesis,” Phys. Rev. X 9, 041017 (2019).
  * Rabinovici  _et al._ (2021) E. Rabinovici, A. Sánchez-Garrido, R. Shir, and J. Sonner, “Operator complexity: A journey to the edge of Krylov space,” Journal of High Energy Physics 2021, 62 (2021).
  * Noh (2021) Jae Dong Noh, “Operator growth in the transverse-field Ising spin chain with integrability-breaking longitudinal field,” Phys. Rev. E 104, 034112 (2021).
  * Dymarsky and Gorsky (2020) Anatoly Dymarsky and Alexander Gorsky, “Quantum chaos as delocalization in Krylov space,” Phys. Rev. B 102, 085137 (2020).
  * (29) Pawel Caputa, Javier M. Magan, and Dimitrios Patramanis, “Geometry of Krylov Complexity,”  arXiv:2109.03824 .
  * (30) Anatoly Dymarsky and Michael Smolkin, “Krylov complexity in conformal field theory,”  arXiv:2104.09514 .
  * Yates  _et al._ (2020) Daniel J. Yates, Alexander G. Abanov, and Aditi Mitra, “Dynamics of almost strong edge modes in spin chains away from integrability,” Phys. Rev. B 102, 195419 (2020).
  * Yates and Mitra (2021) Daniel J. Yates and Aditi Mitra, “Strong and almost strong modes of Floquet spin chains in Krylov subspaces,” Phys. Rev. B 104, 195121 (2021).
  * Pal and Huse (2010) Arijeet Pal and David A. Huse, “Many-body localization phase transition,” Phys. Rev. B 82, 174411 (2010).
  * Basko  _et al._ (2006) D.M. Basko, I.L. Aleiner, and B.L. Altshuler, “Metal–Insulator transition in a weakly interacting many-electron system with localized single-particle states,” Ann. Phys. 321, 1126–1205 (2006).
  * Oganesyan and Huse (2007) Vadim Oganesyan and David A. Huse, “Localization of interacting fermions at high temperature,” Phys. Rev. B 75, 155111 (2007).
  * Nandkishore and Huse (2015) Rahul Nandkishore and David A. Huse, “Many-body localization and thermalization in quantum statistical mechanics,” Annu. Rev. Condens. Matter Phys. 6, 15–38 (2015).
  * He and Lu (2017) Rong-Qiang He and Zhong-Yi Lu, “Characterizing many-body localization by out-of-Time-Ordered correlation,” Phys. Rev. B 95, 054201 (2017).
  * Huang  _et al._ (2017) Yichen Huang, Yong-Liang Zhang, and Xie Chen, “Out-of-Time-Ordered Correlators in Many-Body Localized Systems,” Ann. Phys. 529, 1600318–n/a (2017).
  * Chen  _et al._ (2017) Xiao Chen, Tianci Zhou, David A. Huse, and Eduardo Fradkin, “Out-of-time-order correlations in many-body localized and thermal phases,” Ann. Phys. 529, 1600332 (2017).
  * Slagle  _et al._ (2017) Kevin Slagle, Zhen Bi, Yi-Zhuang You, and Cenke Xu, “Out-of-time-order correlation in marginal many-body localized systems,” Phys. Rev. B 95, 165136 (2017).
  * Swingle and Chowdhury (2017) Brian Swingle and Debanjan Chowdhury, “Slow scrambling in disordered quantum systems,” Phys. Rev. B 95, 060201 (2017).
  * Fan  _et al._ (2017) Ruihua Fan, Pengfei Zhang, Huitao Shen, and Hui Zhai, “Out-of-time-order correlation for many-body localization,” Sci. Bull. 62, 707–711 (2017).
  * Deng  _et al._ (2017) Dong-Ling Deng, Xiaopeng Li, J. H. Pixley, Yang-Le Wu, and S. Das Sarma, “Logarithmic Entanglement Lightcone in Many-Body Localized Systems,” Phys. Rev. B 95, 024202 (2017).
  * Cao (2021) Xiangyu Cao, “A statistical mechanism for operator growth,” J. Phys. Math. Theor. 54, 144001 (2021).
  * Imbrie (2016a) John Z. Imbrie, “Diagonalization and many-body localization for a disordered quantum spin chain,” Phys. Rev. Lett. 117, 027201 (2016a).
  * Imbrie (2016b) John Z. Imbrie, “On Many-Body Localization for Quantum Spin Chains,” J. Stat. Phys. 163, 998–1048 (2016b).
  * Huse  _et al._ (2014) David A. Huse, Rahul Nandkishore, and Vadim Oganesyan, “Phenomenology of fully many-body-localized systems,” Phys. Rev. B 90, 174202 (2014).
  * Viswanath and Müller (1994) V. S. Viswanath and Gerhard Müller,  _The Recursion Method: Applications to Many-Body Dynamics_, Lecture Notes in Physics Monographs (Springer-Verlag, Berlin Heidelberg, 1994).
  * Khait  _et al._ (2016) Ilia Khait, Snir Gazit, Norman Y. Yao, and Assa Auerbach, “Spin transport of weakly disordered Heisenberg chain at infinite temperature,” Phys. Rev. B 93, 224205 (2016).
  * Nandkishore  _et al._ (2014) Rahul Nandkishore, Sarang Gopalakrishnan, and David A. Huse, “Spectral features of a many-body localized system weakly coupled to a heat bath,” Phys. Rev. B 90, 064203 (2014).
  * Johri  _et al._ (2015) Sonika Johri, Rahul Nandkishore, and R. N. Bhatt, “Many-body localization in imperfectly isolated quantum systems,” Phys. Rev. Lett. 114, 117401 (2015).
  * Pekker  _et al._ (2017) David Pekker, Bryan K. Clark, Vadim Oganesyan, and Gil Refael, “Fixed Points of Wegner-Wilson Flows and Many-Body Localization,” Phys. Rev. Lett. 119, 075701 (2017).
  * Varma  _et al._ (2019) Vipin Kerala Varma, Abhishek Raj, Sarang Gopalakrishnan, Vadim Oganesyan, and David Pekker, “Length scales in the many-body localized phase and their spectral signatures,” Phys. Rev. B 100, 115136 (2019).
  * Gopalakrishnan and Parameswaran (2020) Sarang Gopalakrishnan and S. A. Parameswaran, “Dynamics and Transport at the Threshold of Many-Body Localization,” Phys. Rep. 862, 1–62 (2020).
  * (55) J. Herbrych, M. Mierzejewski, and P. Prelovšek, “Relaxation at different length-scales in models of many-body localization,”  arXiv:2110.15635 .
  * Serbyn  _et al._ (2013) Maksym Serbyn, Z. Papić, and Dmitry A. Abanin, “Local conservation laws and the structure of the many-body localized states,” Phys. Rev. Lett. 111, 127201 (2013).
  * Chandran  _et al._ (2015) Anushya Chandran, Isaac H. Kim, Guifre Vidal, and Dmitry A. Abanin, “Constructing local integrals of motion in the many-body localized phase,” Phys. Rev. B 91, 085425 (2015).
  * Trigueros (2021) Fabian Ballar Trigueros,  _Operator Growth and Krylov Complexity in Many-Body Localization_ , Master’s thesis, University of Waterloo, Ontario, Canada (2021).